Modeling quantum bounces of atomic systems

Olivier Rousselle
Laboratoire Kastler Brossel (Sorbonne Université)

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

This research is part of the GBAR collaboration, aiming at timing the free-fall of antihydrogen in the Earth gravitational states. We consider an atomic system (antihydrogen), where we want to produce interference figure at the detection to improve the precision of the experiment. For that, we use a small perfect mirror in which the atoms perform quantum bounces (they are trapped by the combined action of gravity and Casimir-Polder potential). We use Mathematica to model quantum bounces and trajectories, and to plot the intensity signal obtained. Such system can be generalized to other quantum systems (like hydrogen atoms).

Experimental design

We present here the 2D view experimental design used, with the following parameters:
- Initial height h=10 μm;
- Initial horizontal speed of the atoms is v0 = 0.8 m/s;
- Velocity dispersion of the initial wave packet Δv = 0,06 m/s;
- Length of reflecting mirror d = 5 cm;
- Free-fall heigh H ≃ 30 cm; 
- Detection on an horizontal surface.

Atoms bounce several times above the mirror (quantum reflection on Casimir-Polder potential). Quantum paths corresponding to different GQS (Gravitational Quantum States) interfere.

Mathematica code

We model this design by using a semi-classical optical approach. We look for trajectories from the source x0=0, z0=h, t0=0 to the detector (X,T), with n bounces taking place on a perfectly reflecting mirror at zb= 0 at times tk with k = 1,2,...,n.
In the code, we define different functions:

- Classical action, with depends on the initial and detection 	positions.

Scl[z0_,t0_,z_,t_]:=m*(z-z0)^2/(2*(t-t0))-m*g*(t-t0)*(z+z0)/2-m*g^2*(t-t0)^3/24;

- A function which computes the time corresponding to the first bounce (t1) in the reflection mirror. More precisely, we solve a polynomial equation of 4th degree to get t1 and then we eliminate some solutions by applying physical conditions.

findtb[x_,z_,t_,n_]:=(
	td=d*t/x;
	coeffA=n*(n+1)*g^2;
	coeffB=-(2*n+1)*g^2*t;
	coeffC=2*g*(z-zb)+2*(2*n^2-1)*g*(h-zb)+(g*t)^2;
	coeffD=-2*(2*n-1)*(h-zb)*g*t;
	coeffE=4*n*(n-1)*(h-zb)^2;
  t1roots=Reduce[coeffA*y^4+coeffB*y^3+coeffC*y^2+coeffD*y+coeffE==0,y,Reals];
  t1=y/.{ToRules[t1roots]};
	vb=-(zb-h)/t1+g*t1/2;
	tn=t1+2*(n-1)*vb/g;
	t1final={};
For [i=1,i<=Length[t1],i++,
If[t1[[i]]>t0 && tn[[i]]<td &&  tn[[i]]+2*vb[[i]]/g>td,t1final=Append[t1final,t1[[i]]]]];
Return[t1final])

- A function with computes the wave function at the detection. Its arguments are the detection positions (x, z, t) and the number of bounces n on the mirror.

Term[x_,z_,t_,n_]:=(
term=0;
Tb1=findtb[x,z,t,n];
For[k=1,k≤Length[Tb1],k++, 
tb1=Tb1[[k]];
vb=-(zb-h)/tb1+g*tb1/2;
Tb={t0};
For[i=1,i≤n,i++,Tb=Append[Tb,tb1+2*(i-1)*vb/g ]];
Tb=Append[Tb,t];
Z={h};
For[i=1,i≤n,i++,Z=Append[Z,zb]];
Z=Append[Z,z];
ΔS=Sum[Scl[Z[[i]],Tb[[i]],Z[[i+1]],Tb[[i+1]]],{i,1,n+1}];

term=term+Exp[i*ΔS/hbar]];
Return[term])

- The  interference  signal  between  the  semi-classical  trajectories  with 0,1,2,...,n bounces between the source and the detection. We add a factor (-1)^i to take into account π-dephasings at each bounce and a term for the spreading of the wave-function.

intensity[x_,z_,t_,bounces_]:=(
factor=Term0[x,z,t];
factor=factor+Sum[Term[x,z,t,n]*(-1)^n,{n,1,bounces}];Psit=1/Sqrt[2*Pi*(zeta^2+(hbar*t/(2*m*zeta))^2)]*Exp[-((z-h+g*t^2/2)^2)/(2*zeta^2+(hbar*t/(Sqrt[2]*m*zeta))^2)];

Return[Psit*Abs[factor]^2]);

Results

We plot the shape of the intensity on the horizontal detection, as function of x. We clearly see the apparence of interference figure, with fringes.
Plots of the current and number of bounces (in 1D):

bounces=1;
Xmin=d+v0*tau-DeltaX;
Xmax=d+v0*tau+DeltaX;
X=Range[Xmin,Xmax,(Xmax-Xmin)/500];
Plot[intensity[x,-H,x/v0,bounces],{x,Xmin,Xmax},PlotRange->Full,AxesLabel->{x,Current},PlotStyle->Red]

2D plot of the current :

Xmin2D = d+v0*tau-3*DeltaX;
Xmax2D=d+v0*tau+3*DeltaX;
Tmin2D=(d/v0+tau)*0.95;
Tmax2D=(d/v0+tau)*1.05;
DensityPlot[intensity[x,-H,t,bounces],{x,Xmin2D,Xmax2D},{t,Tmin2D,Tmax2D},AxesLabel->Automatic,ColorFunction->”Rainbow”,PlotLegends->Automatic,PlotRange->{0,200}]

(Legended[, (BarLegend)])

Classical case, d=0:
      (Legended[, (BarLegend)])

Thank you for your attention !