Laser Cooling

When we load ions into this trap they are confined along the axis of the trap but, as shown in the previous notebook, they can execute harmonic motion in the xy plane. Our goal is to align ions along the axis and minimize xy motion as much as possible. From the previous 
example we see that we can do this by making (Subscript[v, ox]),  (Subscript[v, oy]) as small as possible. Since velocity can be related to kinetic energy, which in turn,
can be characterized by a temperature, the process of minimizing the kinetic energy of the trapped ion is called cooling. A powerfull
method for accomplishing this is called laser cooling.

Atoms and ions have internal structure, they are called energy levels and in an oversimplistic description we can consider our
ions to have two energy levels as shown in the picture below

Roughly speaking, the ion absorbs a photon with energy  ℏω  that is equal to the difference in the energies of the two 
levels of the ion. Here ℏ is Plank's constant and ω is the angular frequency of the photon which is related to its wavelength λ by
the relation ω = 2 π c/λ where c is the speed of light. In addition to absorbing the energy of the photon the ion can also "gain" a momentum kick from that photon. This momentum is a vector and is given by the expression

Clear[ω]
parameters={n0->10,ω0->3,Γ->0.2,ℏ->1,c->10,mass->1};

nphotons[ω_]= n0 (Γ/2)^2/((ω-ω0)^2+(Γ/2)^2)

Plot[nphotons[ω]/. parameters,{ω,0,6},PlotRange->All]

Note that the function N(ω) is highly peaked near the resonance frequency (Subscript[ω, 0]) and the width of the curve is determined by the
parameter Γ which is also called the full width at half maximum (FWHM). This type of function is called a Lorentzian
and  these parameters depend on which ion one chooses.Though N(ω) determines the fluorescence spectrum, it also determines the efficiency in which a photon of frequency ω can be absorbed by the ion. The closer ω is to the resonance frequency  the more likely the ion will absorb the photon.

Notice also that the integral  ∫N(ω) dω = N((Subscript[ω, 0])) ((π Γ)/2).

Integrate[(Γ/2)^2/((ω-ω0)^2+(Γ/2)^2),{ω,-∞,∞},Assumptions->ω0>0 && Γ>0]

Radiation pressure

(* with the Doppler shifted frequency the radiation force becomes *)
 ℏ ω (Sqrt[1-v/c])/(Sqrt[1+v/c])/c  nphotons[ω(Sqrt[1-v/c])/(Sqrt[1+v/c])]  
ftotal[ω_]=% /. v/c->γ

Normal[Series[ftotal[ω],{γ,0,1}]]

(* f0 is velocity independent *) 
f0[ω_,v_]=FullSimplify[((n0 (Γ^2) ω ℏ)/(c ((Γ^2)+4 (ω^2)-8 ω ω0+4 ((ω0)^2))))]
(* f1 depends linearly on the velocity *)
f1[ω_,v_]=FullSimplify[-((n0 γ (Γ^2) ω ((Γ^2)-4 (ω^2)+4 ((ω0)^2)) ℏ)/(c ((((Γ^2)+4 (ω^2)-8 ω ω0+4 ((ω0)^2)))^2)))] /. γ->v/c

Plot[f0[ω,0] /. parameters,{ω,0,4},PlotRange->All]
Plot3D[(f0[ω,v] /. parameters),{ω,0,4},{v,-1,1},PlotRange->All,PlotPoints->50]

Plot3D[f1[ω,v] /. parameters,{ω,0,4},{v,-1,1},PlotRange->All,PlotPoints->50]

frad[ω_,t_]:= ℏ ω/c (Sqrt[1-x'[t]/c])/(Sqrt[1+x'[t]/c])  nphotons[ω(Sqrt[1-x'[t]/c])/(Sqrt[1+x'[t]/c])] /mass  /. parameters

Omega=3;
omega=ω0 -Γ/. parameters
x0=1;
v0=1;

sols=Flatten[NDSolve[ {x''[t] == - Omega^2 x[t] +frad[omega,t],x[0]==x0,x'[0]==v0},x[t],{t,0,100}]]
fun[t_] = x[t] /. sols

Plot[fun[t],{t,0,100},PlotRange->All]

F0=ℏ omega/c nphotons[omega] /. parameters;
Print[" F0=",F0]
beta=ℏ^2 omega^2/c^2 nphotons'[omega] /. parameters;
Print["β=",beta]
xs=F0/mass/Omega^2 /. parameters;
Print["xs=",xs]
Omegatilde=Sqrt[Omega^2-(beta/2)^2];
Print["Omegatilde=",Omegatilde]

fun2[t_]=xs +  x0 Exp[-beta t/2] Cos[Omegatilde t]+ v0/Omega Exp[-beta t/2] Sin[Omegatilde t];

Plot[fun2[t],{t,0,100},PlotRange->All]