(Subscript[ψ, 0])[x_]=(((((m ω)/(π ℏ))))^(1/4))(E^(-m ω (x^2)/(2ℏ)))
(Subscript[ψ, 1])[x_]=(((((m ω)/(π ℏ))))^(1/4))(Sqrt[((2m ω)/ℏ)])x (E^(-m ω (x^2)/(2ℏ)))

(* Let us check that the normalizations are correct *)

Simplify[(∫_(-∞)^∞)((Abs[(Subscript[ψ, 0])[x]])^2)\[DifferentialD]x, Assumptions->((m ω)/ℏ)>0] 
Simplify[(∫_(-∞)^∞)((Abs[(Subscript[ψ, 1])[x]])^2)\[DifferentialD]x, Assumptions->((m ω)/ℏ)>0]

(* Wigner distributions — using the fact that the wave functions are real *)

(Subscript[ρ, 0])[x_,p_]=Simplify[(1/(2π ℏ))(∫_(-∞)^∞)(E^(-I p r/ℏ))(Subscript[ψ, 0])[x+(r/2)](Subscript[ψ, 0])[x-(r/2)]\[DifferentialD]r, Assumptions->((m ω)/ℏ)>0] 
(Subscript[ρ, 1])[x_,p_]=Simplify[(1/(2π ℏ))(∫_(-∞)^∞)(E^(-I p r/ℏ))(Subscript[ψ, 1])[x+(r/2)](Subscript[ψ, 1])[x-(r/2)]\[DifferentialD]r, Assumptions->((m ω)/ℏ)>0]

(* One can plot these two distributions *)

Plot3D[(Subscript[ρ, 0])[x,p]/.{ℏ-> 1, m-> 1, ω->1}, {x,-2.2,2.2}, {p,-2.2,2.2},PlotRange->All,AxesLabel->Automatic,PlotPoints->50,PlotStyle->Green,ImageSize->600,BaseStyle->{FontFamily->"Times",FontSize->18}]
Plot3D[(Subscript[ρ, 1])[x,p]/.{ℏ-> 1, m-> 1, ω->1}, {x,-2.2,2.2}, {p,-2.2,2.2},PlotRange->All,AxesLabel->Automatic,PlotPoints->50,PlotStyle->Directive[Opacity[0.8],Green],ImageSize->600,BaseStyle->{FontFamily->"Times",FontSize->18}]

(* Here are the density plots of these Wigner distributions *)

DensityPlot[(Subscript[ρ, 0])[x,p]/.{ℏ-> 1, m-> 1, ω->1}, {x,-2.2,2.2}, {p,-2.2,2.2},PlotRange->All,AxesLabel->Automatic,PlotPoints->50,ImageSize->400,BaseStyle->{FontFamily->"Times",FontSize->18},PlotLegends->Automatic]
DensityPlot[(Subscript[ρ, 1])[x,p]/.{ℏ-> 1, m-> 1, ω->1}, {x,-2.2,2.2}, {p,-2.2,2.2},PlotRange->All,AxesLabel->Automatic,PlotPoints->50,ImageSize->400,BaseStyle->{FontFamily->"Times",FontSize->18},PlotLegends->Automatic]

(* Expectation value of the energy / the Hamiltonian *)

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

(* Expectation value of the square of the Hamiltonian *)

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

(* Consider a superposition of the two lower states. One can first compute and plot its Wigner distribution *)

Ψ[x_]=(1/(Sqrt[2]))((Subscript[ψ, 0])[x]+(Subscript[ψ, 1])[x]);
ρ[x_,p_]=Simplify[(1/(2π ℏ))(∫_(-∞)^∞)(E^(-I p r/ℏ))Ψ[x+(r/2)]Ψ[x-(r/2)]\[DifferentialD]r, Assumptions->((m ω)/ℏ)>0]

Plot3D[ρ[x,p]/.{ℏ-> 1, m-> 1, ω->1}, {x,-2.2,2.2}, {p,-2.2,2.2},PlotRange->All,AxesLabel->Automatic,PlotPoints->50,PlotStyle->Green,ImageSize->600,BaseStyle->{FontFamily->"Times",FontSize->18}]

(* … and now compute the expectation values of the Hamiltonian… *)

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

(* … and of the Hamiltonian squared *)

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

(* Thus <(H^2)> ≠ <H((>)^2): for this superposition, there is a spread in energy *)