Wigner.html Wigner.mws

Wigner distribution

The Wigner function is the basic quantity to describe ensembles in a phase-space formulation of quantum mechanics. It is defined as a paricular Fourier transform of the density matrix.

f(q,p) = int(exp(-I*p*s)*psi(q+s/2)*psi^*(q-s/2),s=-infinity..infinity)

(the integral corresponds to an inverse Fourier transform)

Let us start with the simplest examples, such as harmonic oscillator eigenstates:

>    restart; with(inttrans):

>    fx:=x*exp(-x^2/2); #chose the unnormalized state here

fx := x*exp(-1/2*x^2)

>    No:=1/sqrt(int((fx)^2,x=-infinity..infinity));

No := 2^(1/2)/Pi^(1/4)

>    psi:=unapply(No*fx,x);

psi := proc (x) options operator, arrow; 2^(1/2)/Pi^(1/4)*x*exp(-1/2*x^2) end proc

>    int(psi(x)^2,x=-infinity..infinity);

1

>    rho:=psi(x+s/2)*psi(x-s/2);

rho := 2/Pi^(1/2)*(x+1/2*s)*exp(-1/2*(x+1/2*s)^2)*(x-1/2*s)*exp(-1/2*(x-1/2*s)^2)

>    fW:=invfourier(rho,s,p);

fW := (-1+2*p^2+2*x^2)/Pi*exp(-x^2-p^2)

>    Rp:=int(fW,x=-infinity..infinity);

Rp := exp(-p^2)*(-1+2*p^2)/Pi^(1/2)+exp(-p^2)/Pi^(1/2)

>    Rx:=int(fW,p=-infinity..infinity);

Rx := exp(-x^2)*(-1+2*x^2)/Pi^(1/2)+exp(-x^2)/Pi^(1/2)

>    plot(Rx,x=-3..3);

[Maple Plot]

>    int(Rp,p=-infinity..infinity);

1

>    int(Rx,x=-infinity..infinity);

1

>    plot3d(fW,x=-3..3,p=-3..3,axes=boxed,shading=zhue);

[Maple Plot]

Nothing particularly exciting happens here for the ground state. For the first excited state we have something non-classical: the Wigner function is not positive definite!

Exercise 1:

Choose different excited states of the harmonic oscillator and calculate the Wigner function. How is the number of nodes in the wavefunction related to the structure in phase space?

What do we obtain for a wavepacket?

We choose a Gaussian wavepacket centered at the origin with an average momentum of p0 . We pick a simple width (which controls the width of the momentum distribution for all times). We compute the answer at time zero.

>    p0:=1;

p0 := 1

>    fx:=exp(-(x)^2/2)*exp(I*p0*x); #chose the unnormalized state here

fx := exp(-1/2*x^2)*exp(x*I)

>    #fx:=1/(1+x^2)*exp(I*p0*x); # the Lorentzian wavepacket doesn't work.

>    No:=1/sqrt(int(abs(fx)^2,x=-infinity..infinity));

No := 1/(Pi^(1/4))

>    psi:=unapply(No*fx,x);

psi := proc (x) options operator, arrow; 1/Pi^(1/4)*exp(-1/2*x^2)*exp(x*I) end proc

>    int(abs(psi(x))^2,x=-infinity..infinity);

1

>    assume(x,real,s,real);

>    rho:=simplify(expand(psi(x+s/2)*conjugate(psi(x-s/2))));

rho := (-1)^(s/Pi)/Pi^(1/2)*exp(-x^2-1/4*s^2)

>    fW:=invfourier(rho,s,p);

fW := exp(-x^2-2*I*p+1-p^2)/Pi

>    Rp:=int(fW,x=-infinity..infinity);

Rp := exp(-2*I*p+1-p^2)/Pi^(1/2)

>    Rx:=int(fW,p=-infinity..infinity);

Rx := exp(-x^2)/Pi^(1/2)

>    plot([Re(Rp),Im(Rp)],p=-3..3,color=[red,blue]);

[Maple Plot]

Something went wrong. The Wigner function should be real-valued.

>    int(Rx,x=-infinity..infinity);

1

>    fW:=simplify(expand(int(exp(-I*p*s)*rho,s=-infinity..infinity)));

fW := 2*exp(-x^2-1+2*p-p^2)

>    subs(p=1,x=1,%);

2*exp(-1)

>    plot3d(fW,x=-3..3,p=-3..3,axes=boxed,shading=zhue);

[Maple Plot]

The Wigner function for a Gaussian wavepacket state is quite straightforward, and is easily interpreted.

Exercise 2:

Displace the Gaussian wavepacket and re-compute the Wigner function.

Exercise 3:

Broaden the Gaussian wavepacket in coordinate space and re-compute the Wigner function. Make your observation about the momentum distribution.

We explore now the time-dependent Gaussian wavepacket. The width parameter a  determines the spatial extent and momentum spread about the average.

>    a:=1; p0:='p0';

a := 1

p0 := 'p0'

>    psi := unapply(sqrt(a)*1/Pi^(1/4)*exp(-1/2*(x-I*a^2*p0)^2/(a^2+t*I)-1/2*a^2*p0^2)/(sqrt(-a^2-I*t)),x);

psi := proc (x) options operator, arrow; 1/Pi^(1/4)*exp(-1/2*(x-I*p0)^2/(1+t*I)-1/2*p0^2)/(-1-I*t)^(1/2) end proc

>    assume(t>0);

>    int(abs(psi(x))^2,x=-infinity..infinity);

int(1/Pi^(1/2)*exp(1/2*Re(-(x-I*p0)^2/(1+t*I)-p0^2))^2/(1+t^2)^(1/2),x = -infinity .. infinity)

>    int(evalc(simplify(abs(psi(x))^2)),x=-infinity..infinity);

1

>    rho:=simplify(expand(psi(x+s/2)*conjugate(psi(x-s/2))));

rho := 1/Pi^(1/2)*exp(1/4*(-2*I*s*abs(p0)^2*p0+2*s*abs(p0)^2*p0*t+4*I*x*abs(p0)^2*p0-4*x*abs(p0)^2*p0*t-4*I*x*p0^3-4*x*p0^3*t-2*I*s*p0^3-2*s*p0^3*t-4*I*x*s*p0^2*t+4*x^2*p0^2+s^2*p0^2+2*I*p0^4*t+2*p0^4*...
rho := 1/Pi^(1/2)*exp(1/4*(-2*I*s*abs(p0)^2*p0+2*s*abs(p0)^2*p0*t+4*I*x*abs(p0)^2*p0-4*x*abs(p0)^2*p0*t-4*I*x*p0^3-4*x*p0^3*t-2*I*s*p0^3-2*s*p0^3*t-4*I*x*s*p0^2*t+4*x^2*p0^2+s^2*p0^2+2*I*p0^4*t+2*p0^4*...

>    fW:=simplify(expand(int(exp(-I*p*s)*rho,s=-infinity..infinity)));

fW := 2*exp(-1/4*(4*p^2*p0^2*t^2+4*I*p*p0^3*t-8*p*p0^2*x*t-4*I*abs(p0)^2*t*p*p0-4*p*p0*abs(p0)^2+4*I*x*abs(p0)^2*p0+abs(p0)^4+4*x^2*p0^2-4*p*p0^3+2*abs(p0)^2*p0^2+4*p^2*p0^2-4*I*x*p0^3+p0^4)/p0^2)

>    fW:=evalc(fW);

fW := 2*exp((-p^2*p0^2*t^2+2*p*p0^2*x*t+2*p*p0^3-p0^4-x^2*p0^2-p^2*p0^2)/p0^2)

>    p0:=1;

p0 := 1

>    fW1:=subs(t=1,fW);

fW1 := 2*exp(-2*p^2+2*p*x+2*p-1-x^2)

>    fW2:=subs(t=2,fW);

fW2 := 2*exp(-5*p^2+4*p*x+2*p-1-x^2)

>    fW3:=subs(t=3,fW);

fW3 := 2*exp(-10*p^2+6*p*x+2*p-1-x^2)

>    Rx1:=int(fW1,p=-infinity..infinity);

Rx1 := exp(-1/2-1/2*x^2+x)*2^(1/2)*Pi^(1/2)

>    Rx2:=int(fW2,p=-infinity..infinity);

Rx2 := 2/5*exp(-4/5-1/5*x^2+4/5*x)*5^(1/2)*Pi^(1/2)

>    Rx3:=int(fW3,p=-infinity..infinity);

Rx3 := 1/5*exp(-9/10-1/10*x^2+3/5*x)*10^(1/2)*Pi^(1/2)

>    plot([Rx1,Rx2,Rx3],x=-5..10,color=[red,blue,green]);

[Maple Plot]

>    Rp1:=int(fW1,x=-infinity..infinity);

Rp1 := 2*exp(-p^2+2*p-1)*Pi^(1/2)

>    Rp2:=int(fW2,x=-infinity..infinity);

Rp2 := 2*exp(-p^2+2*p-1)*Pi^(1/2)

>    Rp3:=int(fW3,x=-infinity..infinity);

Rp3 := 2*exp(-p^2+2*p-1)*Pi^(1/2)

The momentum distribution remains constant! This is required, since the wavepacket undergoes free dispersion.

>    plots[animate3d](fW,x=-5..15,p=-4..4,t=0..4,axes=boxed,grid=[50,50],shading=zhue,style=patchcontour,frames=20);

[Maple Plot]

The apparent ripples at later times are an artefact of the graphing algorithm.

Exercise 4:

Observe what happens when the spatial width parameter a  is increased or decreased.

>   

>