(Subscript[ρ, 0])[x_,p_]=(1/(π ℏ))(E^(-ξ (p^2)/(ℏ m ω)))(E^(-m ω (((x-x0))^2)/(ℏ ξ)));
ρ[t_,x_,p_]=(Subscript[ρ, 0])[x Cos[ω t]-(p/(m ω))Sin[ω t],p Cos[ω t]+m ω x Sin[ω t]];

(* Let us check that the normalization is correct *)

Simplify[(∫_(-∞)^∞)(∫_(-∞)^∞)(Subscript[ρ, 0])[x,p]\[DifferentialD]x \[DifferentialD]p, Assumptions->{ℏ>0, ξ>0, m>0, ω>0}]

(* Time evolution of the distribution for ξ=1 *)

Plot3D[{ρ[0,x,p],ρ[(π/(2ω)),x,p],ρ[(π/ω),x,p],ρ[((3π)/(2ω)),x,p]}/.{ξ-> 1, x0-> 5,ℏ-> 1, m-> 1, ω->1}, {x,-8,8}, {p,-8,8},PlotRange->All,AxesLabel->Automatic,PlotPoints->50,PlotStyle->Green,ImageSize->600,BaseStyle->{FontFamily->"Times",FontSize->18}]

ListAnimate[Table[Plot3D[ρ[t,x,p]/.{ξ->1,x0->5,ℏ->1,m->1,ω->1}, {x,-8,8}, {p,-8,8},PlotRange->{0,0.32},AxesLabel->Automatic,PlotPoints->50,PlotStyle->Green,ImageSize->600,BaseStyle->{FontFamily->"Times",FontSize->18}],{t,0,((2π)/ω),(π/(10ω))}]]

(* Time evolution of the distribution for ξ=3 *)

Plot3D[{ρ[0,x,p],ρ[(π/(2ω)),x,p],ρ[(π/ω),x,p],ρ[((3π)/(2ω)),x,p]}/.{ξ->3,x0->5,ℏ->1,m->1,ω->1}, {x,-8,8}, {p,-8,8},PlotRange->All,AxesLabel->Automatic,PlotPoints->50,PlotStyle->Green,ImageSize->600,BaseStyle->{FontFamily->"Times",FontSize->18}]

ListAnimate[Table[Plot3D[ρ[t,x,p]/.{ξ->3,x0->5,ℏ->1,m->1,ω->1}, {x,-8,8}, {p,-8,8},PlotRange->{0,0.32},AxesLabel->Automatic,PlotPoints->50,PlotStyle->Green,ImageSize->600,BaseStyle->{FontFamily->"Times",FontSize->18}],{t,0,((2π)/ω),(π/(10ω))}]]

(* Position probability distribution in the case ξ=1 *)

(Subscript[ρ, sym])[t_,x_,p_]=Simplify[ρ[t,x,p]/.ξ-> 1];
 (Subscript[𝒫, sym])[t_,x_]=Simplify[(∫_(-∞)^∞)(Subscript[ρ, sym])[t,x,p]\[DifferentialD]p ,Assumptions-> {ℏ>0, ξ>0, m>0,ω>0}]

ListAnimate[Table[Plot[(Subscript[𝒫, sym])[t,x]/.{x0->5,ℏ->1,m->1,ω->1}, {x,-8,8},PlotRange-> {0,0.58},Filling->Bottom,AxesLabel->Automatic,PlotPoints->50,PlotStyle->Blue,ImageSize->600,BaseStyle->{FontFamily->"Times",FontSize->18}],{t,0,((2π)/ω),(π/(10ω))}]]

(Subscript[ρ, asym])[t_,x_,p_]=Simplify[ρ[t,x,p]/.ξ-> 3];
 (Subscript[𝒫, asym])[t_,x_]=Simplify[(∫_(-∞)^∞)(Subscript[ρ, asym])[t,x,p]\[DifferentialD]p ,Assumptions-> {ℏ>0, ξ>0, m>0,ω>0,t∈Reals}]

ListAnimate[Table[Plot[(Subscript[𝒫, asym])[t,x]/.{x0->5,ℏ->1,m->1,ω->1}, {x,-8,8},PlotRange-> {0,1},Filling->Bottom,AxesLabel->Automatic,PlotPoints->50,PlotStyle->Blue,ImageSize->600,BaseStyle->{FontFamily->"Times",FontSize->18}],{t,0,((2π)/ω),(π/(10ω))}]]