Orca memory equations



In [1]:

    
from fast import *
from fast.config import fast_path
from matplotlib import pyplot
from sympy import sin,cos,exp,sqrt,pi,zeros,I, Integral, oo, integrate
from numpy import array
%matplotlib inline
init_printing()
print_ascii=True; print_ascii=False



In [2]:

    
path="/home/oscar/oxford/inhomogeneous_broadening/complete_model/" 
name='orca'

We will be deriving the optical Bloch equations for a three level system in a ladder configuration as that in the figure.



In [3]:

    
fig=pyplot.figure(); ax=fig.add_subplot(111,aspect="equal")

p1=[0.5,1]; p2=[1.5,3]; p3=[2.5,5]
draw_state(ax,p1,text=r"$|1\rangle$",l=1.0,alignment='right',label_displacement=0.05,fontsize=25,linewidth=4.0)
draw_state(ax,p2,text=r"$|2\rangle$",l=1.0,alignment='right',label_displacement=0.05,fontsize=25,linewidth=4.0)
draw_state(ax,p3,text=r"$|3\rangle$",l=1.0,alignment='right',label_displacement=0.05,fontsize=25,linewidth=4.0)

excitation(ax,[p1[0]+0.25,p1[1]],[p2[0]+0.25,p2[1]], fc="r", ec="r",width=0.01, head_width=0.2, head_length=0.2)
excitation(ax,[p2[0]+0.25,p2[1]],[p3[0]+0.25,p3[1]], fc="b", ec="b",width=0.01, head_width=0.2, head_length=0.2)

decay(     ax,[p1[0]-0.25,p1[1]],[p2[0]-0.25,p2[1]], 0.05,10.0,color="r",linewidth=1.0)
decay(     ax,[p2[0]-0.25,p2[1]],[p3[0]-0.25,p3[1]], 0.05,10.0,color="b",linewidth=1.0)

pyplot.axis('off')
pyplot.savefig(path+name+'_diagram.png',bbox_inches="tight")

We define the number of states and of radiation fields.



In [4]:

    
Ne=3
Nl=2

We define a few important symbols.



In [5]:

    
c,hbar,e,mu0, epsilon0=symbols("c hbar e mu0 varepsilon0",positive=True)
# fprint([c,hbar,e,mu0, epsilon0],print_ascii=print_ascii)



In [6]:

    
t,X,Y,Z,R,Phi=symbols("t X Y Z R Phi",real=True)
RR=Matrix([R*cos(Phi),R*sin(Phi),Z])
# RR=Matrix([X,Y,Z])
# fprint([t,RR],print_ascii=print_ascii)

We define the variables related to the laser field.



In [7]:

    
E0,omega_laser=define_laser_variables(Nl, variables=[t,R,Z])
# fprint(E0,print_ascii=print_ascii)



In [8]:

    
# fprint(omega_laser,print_ascii=print_ascii)

We write two electric fields propagating trough the $\hat{x}$ direction polarized in the $\hat{z}$ direction. First the wave vectors:



In [9]:

    
phi1=0 ; theta1=0; alpha1=pi/2; beta1=pi/8
phi2=0; theta2=pi; alpha2=pi/2; beta2=pi/8

k1=Matrix([cos(phi1)*sin(theta1),sin(phi1)*sin(theta1),cos(theta1)])
k2=Matrix([cos(phi2)*sin(theta2),sin(phi2)*sin(theta2),cos(theta2)])

k=[k1,k2]

# fprint(k,print_ascii=print_ascii)

The polarization vectors.



In [10]:

    
ep1=polarization_vector(phi1,theta1,alpha1,beta1, 1)
ep2=polarization_vector(phi2,theta2,alpha2,beta2, 1)

em1=ep1.conjugate()
em2=ep2.conjugate()

ep=[ep1,ep2]
em=[em1,em2]

# fprint([ep,em],print_ascii=print_ascii)

The electric field (evaluated in $\vec{R}=0$).



In [11]:

    
zero_vect=Matrix([0,0,0])
E_cartesian = [(+E0[l]            *ep[l]*exp(-I*omega_laser[l]*(t-k[l].dot(RR)/c)) 
                +E0[l].conjugate()*em[l]*exp(+I*omega_laser[l]*(t-k[l].dot(RR)/c)))/2 
                    for l in range(Nl)]

# fprint(E_cartesian,print_ascii=print_ascii)



In [12]:

    
l1=PlaneWave(phi1,theta1,alpha1,beta1,color="red")
l2=PlaneWave(phi2,theta2,alpha2,beta2,color="blue")

laseres=[l1,l2]
Nl=len(laseres)

fig = pyplot.figure(); ax = fig.gca(projection='3d')
draw_lasers_3d(ax,laseres,path+name+'_lasers.png')

We write the electric fields in the helicity basis (see notebook "Vectors in the helicity basis and the electric field").



In [13]:

    
E=[cartesian_to_helicity(E_cartesian[l]).expand() for l in range(Nl)]
# fprint(E,print_ascii=print_ascii)

We define the position operator.



In [14]:

    
r=define_r_components(Ne,helicity=True,explicitly_hermitian=True)
#Ladder means that r_{p;31}=0
r=[ri.subs({r[0][2,0]:0,r[1][2,0]:0,r[2][2,0]:0}) for ri in r]
# fprint(r,print_ascii=print_ascii)

The frequencies of the energy levels, the resonant frequencies, and the decay frequencies.



In [15]:

    
omega_level,omega,gamma=define_frequencies(Ne,explicitly_antisymmetric=True)
#Ladder means gamma31=0
gamma=gamma.subs({gamma[2,0]:0})

# fprint(omega_level,print_ascii=print_ascii)



In [16]:

    
# fprint(omega, print_ascii=print_ascii)



In [17]:

    
# fprint(gamma, print_ascii=print_ascii)

The atomic hamiltonian is



In [18]:

    
H0=Matrix([[hbar*omega_level[i]*KroneckerDelta(i,j) for j in range(Ne)] for i in range(Ne)])
# fprint(H0, print_ascii=print_ascii)

The interaction hamiltonian is



In [19]:

    
zero_matrix=zeros(Ne,Ne)
H1=sum([ e*helicity_dot_product(E[l],r) for l in range(Nl)],zero_matrix)
# fprint(H1,print_ascii=print_ascii)

and the complete hamiltonian is



In [20]:

    
H=H0+H1

Rotating wave approximation

Notice that the electric field can be separated by terms with positive and negative frequency:



In [21]:

    
E_cartesian_p=[E0[l]/2            *ep[l]*exp(-I*omega_laser[l]*(t-k[l].dot(RR)/c)) for l in range(Nl)]
E_cartesian_m=[E0[l].conjugate()/2*em[l]*exp(I*omega_laser[l]*(t-k[l].dot(RR)/c)) for l in range(Nl)]

E_p=[cartesian_to_helicity(E_cartesian_p[l]) for l in range(Nl)]
E_m=[cartesian_to_helicity(E_cartesian_m[l]) for l in range(Nl)]

# fprint([E_p,E_m], print_ascii=print_ascii)



In [22]:

    
# fprint( simplify(sum([E[l] for l in range(Nl)],zero_vect)-(sum([E_p[l]+E_m[l] for l in range(Nl)],zero_vect) )), print_ascii=print_ascii)

The position operator can also be separated in this way. We go to the interaction picture (with $\hat{H}_0$ as the undisturbed hamiltonian)



In [23]:

    
r_I=[ Matrix([[exp(I*omega[i,j]*t)*r[p][i,j] for j in range(Ne)] for i in range(Ne)]) for p in range(3)]
# fprint(r_I[0], print_ascii=print_ascii)



In [24]:

    
# fprint(r_I[1], print_ascii=print_ascii)



In [25]:

    
# fprint(r_I[2], print_ascii=print_ascii)

Which can be decomposed in positive and negative frequencies as



In [26]:

    
r_I_p=[ Matrix([[ delta_greater(j,i)*exp(-I*omega[j,i]*t)*r[p][i,j] for j in range(Ne)]for i in range(Ne)]) for p in range(3)]
# fprint(r_I_p[0], print_ascii=print_ascii)



In [27]:

    
# fprint(r_I_p[1], print_ascii=print_ascii)



In [28]:

    
# fprint(r_I_p[2], print_ascii=print_ascii)



In [29]:

    
r_I_m=[ Matrix([[ delta_lesser( j,i)*exp( I*omega[i,j]*t)*r[p][i,j] for j in range(Ne)]for i in range(Ne)]) for p in range(3)]
# fprint(r_I_m[0],print_ascii=print_ascii)



In [30]:

    
# fprint(r_I_m[1], print_ascii=print_ascii)



In [31]:

    
# fprint(r_I_m[2], print_ascii=print_ascii)

that summed equal $\vec{\hat{r}}_I$



In [32]:

    
# fprint( [r_I[p]-(r_I_p[p]+r_I_m[p]) for p in range(3)] , print_ascii=print_ascii)

Thus the interaction hamiltonian in the interaciton picture is \begin{equation} \hat{H}_{1I}=e\vec{E}\cdot \vec{\hat{r}}_I= e(\vec{E}^{(+)}\cdot \vec{\hat{r}}^{(+)}_I + \vec{E}^{(+)}\cdot \vec{\hat{r}}^{(-)}_I + \vec{E}^{(-)}\cdot \vec{\hat{r}}^{(+)}_I + \vec{E}^{(-)}\cdot \vec{\hat{r}}^{(-)}_I) \end{equation}



In [33]:

    
H1I=sum([ e*helicity_dot_product(E[l],r_I) for l in range(Nl)],zero_matrix)
# fprint(H1I,print_ascii=print_ascii)

Since both $\omega^l$ and $\omega_{ij}$ are in the order of THz, the terms that have frequencies with the same sign are summed, and thus also of the order of THz. The frequencies in the terms with oposite signs however, are detunings of the order of MHz. Since we are only interested in the coarse-grained evolution of the density matrix, we may omit the fast terms and approximate

\begin{equation} \hat{H}_{1I} \simeq \hat{H}_{1I,RWA}= e( \vec{E}^{(+)}\cdot \vec{\hat{r}}^{(-)}_I + \vec{E}^{(-)}\cdot \vec{\hat{r}}^{(+)}_I ) \end{equation}

That is known as the rotating wave approximation (RWA).



In [34]:

    
H1IRWA=sum( [ (e*(helicity_dot_product(E_p[l],r_I_m)+helicity_dot_product(E_m[l],r_I_p))) for l in range(Nl)],zero_matrix)
# fprint(H1IRWA,print_ascii=print_ascii)

The matrix element $(\hat{H}_{1I,RWA})_{21}$ element is



In [35]:

    
# fprint(H1IRWA[1,0].expand(),print_ascii=print_ascii)

But if the detuning $\omega_{21}-\omega^1 \ll \omega_{21}-\omega^2$ (the second field is far detuned from the $1 \rightarrow 2$ transition), then $\omega_{21}-\omega^2$ may be also considered too high a frequency to be relevant to coarse-grained evolution. So we might neclect that term in $(\hat{H}_{1I,RWA})_{21}$ and similarly neglect the $\omega_{32}-\omega^1$ for term in $(\hat{H}_{1I,RWA})_{32}$:



In [36]:

    
# fprint(H1IRWA[2,1].expand(),print_ascii=print_ascii)

In other words, if the detunings in our experiments allow the approximmation, we might choose which frequency components $\omega^l$ excite which transitions. Let us say that $L_{ij}$ is the set of $l$ such that $\omega^l$ excites the transition $i\rightarrow j$



In [37]:

    
Lij=[[1,2,[1]],[2,3,[2]]]
Lij=formatLij(Lij,Ne)
print array(Lij)









    



[[[] [1] []]
 [[1] [] [2]]
 [[] [2] []]]

Thus the interacion hamiltonian in the interaction picture can be approximated as



In [38]:

    
H1IRWA =sum([ e*( helicity_dot_product( E_p[l],vector_element(r_I_m,i,j)) ) * ket(i+1,Ne)*bra(j+1,Ne) 
            for l in range(Nl) for j in range(Ne) for i in range(Ne) if l+1 in Lij[i][j] ],zero_matrix)
H1IRWA+=sum([ e*( helicity_dot_product( E_m[l],vector_element(r_I_p,i,j)) ) * ket(i+1,Ne)*bra(j+1,Ne) 
            for l in range(Nl) for j in range(Ne) for i in range(Ne) if l+1 in Lij[i][j] ],zero_matrix)

# fprint(H1IRWA, print_ascii=print_ascii)

Returning to the Schrödinger picture we have.



In [39]:

    
r_p=[ Matrix([[ delta_greater(j,i)*r[p][i,j] for j in range(Ne)]for i in range(Ne)]) for p in range(3)]
# fprint(r_p, print_ascii=print_ascii)



In [40]:

    
r_m=[ Matrix([[ delta_lesser( j,i)*r[p][i,j] for j in range(Ne)]for i in range(Ne)]) for p in range(3)]
# fprint(r_m, print_ascii=print_ascii)



In [41]:

    
# fprint( [r[p]-(r_p[p]+r_m[p]) for p in range(3)] , print_ascii=print_ascii)

Thus the interaction hamiltonian in the Schrödinger picture in the rotating wave approximation is



In [42]:

    
H1RWA =sum([ e*( helicity_dot_product( E_p[l],vector_element(r_m,i,j)) ) * ket(i+1,Ne)*bra(j+1,Ne) 
            for l in range(Nl) for j in range(Ne) for i in range(Ne) if l+1 in Lij[i][j] ],zero_matrix)
H1RWA+=sum([ e*( helicity_dot_product( E_m[l],vector_element(r_p,i,j)) ) * ket(i+1,Ne)*bra(j+1,Ne) 
            for l in range(Nl) for j in range(Ne) for i in range(Ne) if l+1 in Lij[i][j] ],zero_matrix)

# fprint(H1RWA, print_ascii=print_ascii)

And the complete hamiltonian in the Schrödinger picture in the rotating wave approximation is



In [43]:

    
HRWA=H0+H1RWA
# fprint(HRWA, print_ascii=print_ascii)

Rotating Frame

Next we will make a phase transformation in order to eliminate the explicit time dependance of the equations.



In [44]:

    
cc,cctilde,phase=define_psi_coefficients(Ne)
# fprint([cc,cctilde,phase], print_ascii=print_ascii)



In [45]:

    
phase=Matrix([ Function("theta_"+str(i+1),real=True)(t,Z) for i in range(Ne)])
# phase



In [46]:

    
psi=Matrix([ exp(I*phase[i])*cctilde[i] for i in range(Ne)])
# fprint(psi, print_ascii=print_ascii)

The Schrödinger equation $i\hbar \partial_t |\psi\rangle=\hat{H}_{RWA}$ is



In [47]:

    
lhs=Matrix([(I*hbar*Derivative(psi[i],t).doit()).expand() for i in range(Ne)])
# fprint(lhs, print_ascii=print_ascii)



In [48]:

    
rhs=HRWA*psi

We multiply each of these equations by $e^{-i \theta_i}$ and substracting $i \theta_i \tilde{c}_i$



In [49]:

    
lhs_new=Matrix([simplify(  lhs[i]*exp(-I*phase[i]) +hbar*Derivative(phase[i],t)*cctilde[i] ) for i in range(Ne)])
# fprint(lhs_new, print_ascii=print_ascii)



In [50]:

    
rhs_new=Matrix([simplify(  rhs[i]*exp(-I*phase[i])
                         +hbar*Derivative(phase[i],t)*cctilde[i] ).expand() for i in range(Ne)])
# fprint(rhs_new, print_ascii=print_ascii)

It can be seen that the equations loose their explicit time dependance only if $\omega^{1}(-Z/c+t) - \theta_{1} + \theta_{2}=0$ and $\omega^{2}(Z/c+t) - \theta_{2} + \theta_{3}=0$. Which is satisfied if



In [51]:

    
eq1=omega_laser[0]*(-Z/c+t)+phase[1]-phase[0]
eq2=omega_laser[1]*( Z/c+t)+phase[2]-phase[1]
pt=solve([eq1,eq2],[phase[1],phase[2]])
# pt

Thus the equations become



In [52]:

    
# fprint(lhs_new, print_ascii=print_ascii)



In [53]:

    
rhs_new=simplify(rhs_new.subs(pt)).expand().simplify()
# fprint(rhs_new, print_ascii=print_ascii)

It can be seen that this is the Schrödinger equation derived from an effective hamiltonian $\tilde{H}$



In [54]:

    
Htilde=Matrix([ [Derivative(rhs_new[i],cctilde[j]).doit().simplify() for j in range(Ne)] for i in range(Ne)])
# fprint(Htilde, print_ascii=print_ascii)

We can see that it is convenient to choose $\theta_1=-\omega_1$ to simplify the hamiltonian. Also, we can recognize $\omega^1-\omega_2+\omega_1=\delta^1$ as the detuning of the first field relative to the atomic transition $\omega_{21}=\omega_2-\omega_1$, and the same for $\omega^2-\omega_3+\omega_2=\delta^2$. And choosing $\theta_1=\omega_1 t$



In [55]:

    
delta1,delta2=symbols("delta1 delta2",real=True)
Htilde=Htilde.subs({phase[0]:-omega_level[0]*t}).doit()
Htilde=Htilde.subs({omega_laser[0]:delta1+omega_level[1]-omega_level[0]})
Htilde=Htilde.subs({omega_laser[1]:delta2+omega_level[2]-omega_level[1]})

Htilde=Htilde.expand()

# fprint(Htilde, print_ascii=print_ascii)

If we define the Rabi frequencies $\Omega_1 =e E_0^1 r_{0;21}/\hbar$ and $\Omega_2 =e E_0^2 r_{0;32}/\hbar$



In [56]:

    
Omega1,Omega2=symbols("Omega1 Omega2",real=True)
Omega1,Omega2=symbols("Omega1 Omega2")
Omega1=Function("Omega1")(t,R,Z)
Omega2=Function("Omega2")(t,R,Z)

Htilde=Htilde.subs({E0[0]:Omega1*hbar/r[2][1,0]/e})
Htilde=Htilde.subs({E0[1]:Omega2*hbar/r[0][2,1]/e})

fprint(Htilde, print_ascii=print_ascii)









    Out[56]:





$$\left[\begin{matrix}0 & \frac{\hbar}{2} \overline{\Omega_{1}{\left (t,R,Z \right )}} & 0\\\frac{\hbar}{2} \Omega_{1}{\left (t,R,Z \right )} & - \delta_{1} \hbar & \frac{\hbar}{2} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\0 & \frac{\hbar}{2} \Omega_{2}{\left (t,R,Z \right )} & - \delta_{1} \hbar - \delta_{2} \hbar\end{matrix}\right]$$

Optical Bloch Equations

We define the density matrix.



In [57]:

    
rho = define_density_matrix(Ne, variables=[t,R,Z])
fprint(rho , print_ascii=print_ascii)









    Out[57]:





$$\left[\begin{matrix}\rho_{11}{\left (t,R,Z \right )} & \rho_{12}{\left (t,R,Z \right )} & \rho_{13}{\left (t,R,Z \right )}\\\rho_{21}{\left (t,R,Z \right )} & \rho_{22}{\left (t,R,Z \right )} & \rho_{23}{\left (t,R,Z \right )}\\\rho_{31}{\left (t,R,Z \right )} & \rho_{32}{\left (t,R,Z \right )} & \rho_{33}{\left (t,R,Z \right )}\end{matrix}\right]$$

The hamiltonian part of the equations is \begin{equation} \dot{\hat{\rho}}=\frac{i}{\hbar}[\hat{\rho}, \hat{\tilde{H}}] \end{equation}



In [58]:

    
hamiltonian_terms=(I/hbar*(rho*Htilde-Htilde*rho)).expand()
# fprint(hamiltonian_terms, print_ascii=print_ascii)

There are two Lindblad operators, since there are two spontaneous decay channels.



In [59]:

    
lindblad_terms =gamma[1,0]*lindblad_operator(ket(1,Ne)*bra(2,Ne),rho)
lindblad_terms+=gamma[2,1]*lindblad_operator(ket(2,Ne)*bra(3,Ne),rho)

# fprint(lindblad_terms, print_ascii=print_ascii)

The Optical Bloch equations are thus.



In [60]:

    
eqs = hamiltonian_terms + lindblad_terms



In [61]:

    
eqsign=symbols("=")
eqs_list=[]
for mu in range(0,Ne**2-1 -(Ne**2 - Ne)/2+1):
    ii,jj,s=IJ(mu,Ne)
    i=ii-1; j=jj-1
    eqs_list+=[[Derivative(rho[i,j],t),eqsign,eqs[i,j]]]
eqs_list=Matrix(eqs_list)
fprint(eqs_list, print_ascii=print_ascii)









    Out[61]:





$$\left[\begin{matrix}\frac{\partial}{\partial t} \rho_{11}{\left (t,R,Z \right )} & = & \gamma_{21} \rho_{22}{\left (t,R,Z \right )} + \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{12}{\left (t,R,Z \right )} - \frac{i}{2} \rho_{21}{\left (t,R,Z \right )} \overline{\Omega_{1}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{22}{\left (t,R,Z \right )} & = & - \gamma_{21} \rho_{22}{\left (t,R,Z \right )} + \gamma_{32} \rho_{33}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{12}{\left (t,R,Z \right )} + \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{23}{\left (t,R,Z \right )} + \frac{i}{2} \rho_{21}{\left (t,R,Z \right )} \overline{\Omega_{1}{\left (t,R,Z \right )}} - \frac{i}{2} \rho_{32}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{33}{\left (t,R,Z \right )} & = & - \gamma_{32} \rho_{33}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{23}{\left (t,R,Z \right )} + \frac{i}{2} \rho_{32}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{21}{\left (t,R,Z \right )} & = & i \delta_{1} \rho_{21}{\left (t,R,Z \right )} - \frac{\gamma_{21}}{2} \rho_{21}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{11}{\left (t,R,Z \right )} + \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{22}{\left (t,R,Z \right )} - \frac{i}{2} \rho_{31}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{31}{\left (t,R,Z \right )} & = & i \delta_{1} \rho_{31}{\left (t,R,Z \right )} + i \delta_{2} \rho_{31}{\left (t,R,Z \right )} - \frac{\gamma_{32}}{2} \rho_{31}{\left (t,R,Z \right )} + \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{32}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{21}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t} \rho_{32}{\left (t,R,Z \right )} & = & i \delta_{2} \rho_{32}{\left (t,R,Z \right )} - \frac{\gamma_{21}}{2} \rho_{32}{\left (t,R,Z \right )} - \frac{\gamma_{32}}{2} \rho_{32}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{22}{\left (t,R,Z \right )} + \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{33}{\left (t,R,Z \right )} + \frac{i}{2} \rho_{31}{\left (t,R,Z \right )} \overline{\Omega_{1}{\left (t,R,Z \right )}}\end{matrix}\right]$$

Wave equation

From Maxwell's equations in a dielectric medium it can be shown that in the abscence of bound charges, and magnetization currents, the electric field and the polazation of the medium follow the inhomogeneous wave equation.

\begin{equation} \nabla^2 \vec{E} - \frac{1}{c^2} \partial^2_t \vec{E}= \mu_0 \partial_t^2 \vec{P} \end{equation}

We have also taken our fields to be of the form

\begin{equation} \vec{E}(t,\vec{R})=\vec{E}^{+}+\vec{E}^{-}=\frac{1}{2} \sum_l \vec{E}^{l(+)}(t,\vec{R}) e^{i(\vec{k}^l \cdot \vec{R} -\omega^l t)} +c.c. \end{equation}

The $(+)$ part of the field is explicitly



In [62]:

    
E_cartesian_p









    Out[62]:





$$\left [ \left[\begin{matrix}- \frac{\sqrt{2}}{4} \operatorname{E_{01}}{\left (t,R,Z \right )} e^{- i \varpi_{1} \left(- \frac{Z}{c} + t\right)}\\- \frac{\sqrt{2} i}{4} \operatorname{E_{01}}{\left (t,R,Z \right )} e^{- i \varpi_{1} \left(- \frac{Z}{c} + t\right)}\\0\end{matrix}\right], \quad \left[\begin{matrix}\frac{\sqrt{2}}{4} \operatorname{E_{02}}{\left (t,R,Z \right )} e^{- i \varpi_{2} \left(\frac{Z}{c} + t\right)}\\- \frac{\sqrt{2} i}{4} \operatorname{E_{02}}{\left (t,R,Z \right )} e^{- i \varpi_{2} \left(\frac{Z}{c} + t\right)}\\0\end{matrix}\right]\right ]$$

introducing this into the wave equation we get



In [63]:

    
def laplacian_cylindric(scalar,coords,full=False):
    R,Phi,Z=coords
    return Derivative(scalar,Z,2).doit()

def laplacian_vec_cylindric(vector,coords,full=False):
    return Matrix([ laplacian_cylindric(vi,coords,full=full) for vi in vector])



In [64]:

    
E_cartesian_p_tot = sum([E_cartesian_p[l] for l in range(Nl)], zero_vect)
E_cartesian_p_tot









    Out[64]:





$$\left[\begin{matrix}- \frac{\sqrt{2}}{4} \operatorname{E_{01}}{\left (t,R,Z \right )} e^{- i \varpi_{1} \left(- \frac{Z}{c} + t\right)} + \frac{\sqrt{2}}{4} \operatorname{E_{02}}{\left (t,R,Z \right )} e^{- i \varpi_{2} \left(\frac{Z}{c} + t\right)}\\- \frac{\sqrt{2} i}{4} \operatorname{E_{01}}{\left (t,R,Z \right )} e^{- i \varpi_{1} \left(- \frac{Z}{c} + t\right)} - \frac{\sqrt{2} i}{4} \operatorname{E_{02}}{\left (t,R,Z \right )} e^{- i \varpi_{2} \left(\frac{Z}{c} + t\right)}\\0\end{matrix}\right]$$



In [65]:

    
term1=laplacian_vec_cylindric(E_cartesian_p_tot,[R,Phi,Z]) 
term2=-1/c**2*Matrix([Derivative(vi,t,2).doit() for vi in E_cartesian_p_tot])
lhs=term1+term2
#pprint(lhs,num_columns=150)

And if we consider the amplitudes as slowly varying envelopes (both in time and space) we can approximate them as



In [66]:

    
svea_subs={Derivative(E0[0],Z,2):0,Derivative(E0[1],Z,2):0,
           Derivative(E0[0],t,2):0,Derivative(E0[1],t,2):0}
svea_subs









    Out[66]:





$$\left \{ \frac{\partial^{2}}{\partial Z^{2}}  \operatorname{E_{01}}{\left (t,R,Z \right )} : 0, \quad \frac{\partial^{2}}{\partial t^{2}}  \operatorname{E_{01}}{\left (t,R,Z \right )} : 0, \quad \frac{\partial^{2}}{\partial Z^{2}}  \operatorname{E_{02}}{\left (t,R,Z \right )} : 0, \quad \frac{\partial^{2}}{\partial t^{2}}  \operatorname{E_{02}}{\left (t,R,Z \right )} : 0\right \}$$



In [67]:

    
lhs=lhs.subs(svea_subs)
lhs.simplify()
# lhs

On the other hand, we may approximate the macroscopic polarization of the atoms as varying only at the frequencies of the electric field components and at the same polarizations:

\begin{equation} \vec{P}=\frac{1}{2}\sum_l \vec{P}^{l(+)} e^{i(\vec{k}^l \cdot \vec{R} -\omega^l t)} +c.c. \end{equation}



In [68]:

    
P0 = [Function("P_0^1")(t,R,Z), Function("P_0^2")(t,R,Z)]
# P0



In [69]:

    
P_cartesian_p=[P0[l]/2            *ep[l]*exp(-I*omega_laser[l]*(t-k[l].dot(RR)/c)) for l in range(Nl)]
P_cartesian_m=[P0[l].conjugate()/2*em[l]*exp(I*omega_laser[l]*(t-k[l].dot(RR)/c)) for l in range(Nl)]

P_cartesian_p_tot=sum(P_cartesian_p, zero_vect)

# P_cartesian_p

The right-hand side of the wave equation



In [70]:

    
rhs=mu0*Matrix([Derivative(vi,t,2).doit() for vi in P_cartesian_p_tot])
#pprint(rhs)

And in another slowly varying approximation, the terms with $(\omega^1)^2$ are much larger than those with derivatives of the amplitudes.



In [71]:

    
svea_subs2={Derivative(P0[0],t,1):0, Derivative(P0[1],t,1):0,
            Derivative(P0[0],t,2):0, Derivative(P0[1],t,2):0}
# svea_subs2



In [72]:

    
rhs=rhs.subs(svea_subs2)
# pprint(rhs)



In [73]:

    
eqs_wave=lhs-rhs
# pprint(eqs_wave)

These are three scalar equations each of which has coupled terms varying at high frequencies $\omega^l$. However, these frequencies can be decoupled if the polarization of the beams are orthogonal (as they are in our case). Taking a dot product with the polarizations we obtain one equation for each frequency component.



In [74]:

    
fact1=2*c**2*exp(I*omega_laser[0]*(t-Z/c))/(2*I*omega_laser[0])
lhs1=(fact1*ep[0].conjugate().dot(lhs)).expand()
rhs1=(fact1*ep[0].conjugate().dot(rhs)).expand()

# lhs1, rhs1



In [75]:

    
fact2=2*c**2*exp(I*omega_laser[1]*(t+Z/c))/(2*I*omega_laser[1])
lhs2=(fact2*ep[1].conjugate().dot(lhs)).expand()
rhs2=(fact2*ep[1].conjugate().dot(rhs)).expand()

# lhs2, rhs2

We take the spacial derivatives to the right-hand side. Just because we can.



In [76]:

    
lhs1=lhs1-c*Derivative(E0[0],Z)
rhs1=rhs1-c*Derivative(E0[0],Z)

lhs2=lhs2+c*Derivative(E0[1],Z)
rhs2=rhs2+c*Derivative(E0[1],Z)

We put these equations in terms of Rabi frequencies.



In [77]:

    
fact1=e*r[2][1,0]/hbar
lhs1=(fact1*lhs1.subs({E0[0]:Omega1*hbar/r[2][1,0]/e}).doit()).expand()
rhs1=(fact1*rhs1.subs({E0[0]:Omega1*hbar/r[2][1,0]/e}).doit()).expand()

Matrix([lhs1,eqsign,rhs1]).transpose()









    Out[77]:





$$\left[\begin{matrix}\frac{\partial}{\partial t} \Omega_{1}{\left (t,R,Z \right )} & = & \frac{i e \mu_{0} r_{+1;21}}{2 \hbar} c^{2} \varpi_{1} \operatorname{P^{1}_{0}}{\left (t,R,Z \right )} - c \frac{\partial}{\partial Z} \Omega_{1}{\left (t,R,Z \right )}\end{matrix}\right]$$



In [78]:

    
fact2=e*r[0][2,1]/hbar
lhs2=(fact2*lhs2.subs({E0[1]:Omega2*hbar/r[0][2,1]/e}).doit()).expand()
rhs2=(fact2*rhs2.subs({E0[1]:Omega2*hbar/r[0][2,1]/e}).doit()).expand()

Matrix([lhs2,eqsign,rhs2]).transpose()









    Out[78]:





$$\left[\begin{matrix}\frac{\partial}{\partial t} \Omega_{2}{\left (t,R,Z \right )} & = & \frac{i e \mu_{0} r_{-1;32}}{2 \hbar} c^{2} \varpi_{2} \operatorname{P^{2}_{0}}{\left (t,R,Z \right )} + c \frac{\partial}{\partial Z} \Omega_{2}{\left (t,R,Z \right )}\end{matrix}\right]$$

We may relate the macroscopic polarization to the density matrix if we identify the quantum mechanics operator that corresponds to the polarization. Since the polarization is nothing but the density of dipole moment, we can see that

\begin{equation} \vec{P}=-n \mathrm{Tr}(e \vec{\hat{r}} \hat{\rho}) \end{equation}

notice that the minus sign comes from the fact that $\vec{\hat{r}}$ points in the direction of the electron relative to the proton, while the electric dipole moment points in the opposite direction. If we further make the asumption that each frequency component $l$ of the polarization is only driven by the transition $|i\rangle \leftrightarrow |j\rangle$ if $l\in L_{ij}$, then we may decompose $\vec{\hat{r}}$ into $l$ components as.

\begin{equation} \vec{\hat{r}} = \sum_l \vec{\hat{r}}^{l(+)} + \vec{\hat{r}}^{l(+)} \end{equation}

Explicitly:



In [79]:

    
r_p_component=[ [ Matrix([ [ r_p[p][i,j] if l+1 in Lij[i][j] else 0
                            for j in range(Ne)  ] for i in range(Ne)])
                 for p in range(3)] for l in range(Nl)]

r_m_component=[ [ Matrix([ [ r_m[p][i,j] if l+1 in Lij[i][j] else 0
                            for j in range(Ne)  ] for i in range(Ne)])
                 for p in range(3)] for l in range(Nl)]



In [80]:

    
rpl1=r_p_component[0]
rpl2=r_p_component[1]



In [81]:

    
n=Function("n")(R,Z)
n









    Out[81]:





$$n{\left (R,Z \right )}$$



In [82]:

    
vh=Matrix(symbols("r_-,r0,r_+",real=True))
# vh,helicity_to_cartesian(vh)



In [83]:

    
rpl1_cartesian=[(rpl1[0]-rpl1[2])/sqrt(2),(rpl1[0]+rpl1[2])*I/sqrt(2),rpl1[1]]
# rpl1_cartesian



In [84]:

    
rpl2_cartesian=[(rpl2[0]-rpl2[2])/sqrt(2),(rpl2[0]+rpl2[2])*I/sqrt(2),rpl2[1]]
# rpl2_cartesian

The following factor of 2 comes from the fact that $P^{l(+)}/2 = -e n Tr(\vec{\hat{r}}^{l(+)}\hat{\rho})$



In [85]:

    
Ppl1=-2*n*e*Matrix([ (rpl1_cartesian[i]*rho).trace() for i in range(3)])
Ppl2=-2*n*e*Matrix([ (rpl2_cartesian[i]*rho).trace() for i in range(3)])

# Ppl1, Ppl2

Taking the dot product with the polarizations we get the polarization amplitudes in terms of density matrix elements.



In [86]:

    
Ppl1=ep[0].conjugate().dot(Ppl1).expand()
Ppl2=ep[1].conjugate().dot(Ppl2).expand()

# Ppl1, Ppl2



In [87]:

    
rhs1=rhs1.subs({P0[0]:Ppl1})
rhs2=rhs2.subs({P0[1]:Ppl2})

Maxwell-Bloch equations

And we add these equations to the Bloch equations



In [88]:

    
eqsign=symbols("=")
eqs_list=[]
for mu in range(0,Ne**2-1 -(Ne**2 - Ne)/2+1):
    ii,jj,s=IJ(mu,Ne)
    i=ii-1; j=jj-1
    eqs_list+=[[Derivative(rho[i,j],t),eqsign,eqs[i,j]]]

eqs_list += [[lhs1, eqsign, rhs1]]
eqs_list += [[lhs2, eqsign, rhs2]]

eqs_list=Matrix(eqs_list)
fprint(eqs_list, print_ascii=print_ascii)









    Out[88]:





$$\left[\begin{matrix}\frac{\partial}{\partial t} \rho_{11}{\left (t,R,Z \right )} & = & \gamma_{21} \rho_{22}{\left (t,R,Z \right )} + \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{12}{\left (t,R,Z \right )} - \frac{i}{2} \rho_{21}{\left (t,R,Z \right )} \overline{\Omega_{1}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{22}{\left (t,R,Z \right )} & = & - \gamma_{21} \rho_{22}{\left (t,R,Z \right )} + \gamma_{32} \rho_{33}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{12}{\left (t,R,Z \right )} + \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{23}{\left (t,R,Z \right )} + \frac{i}{2} \rho_{21}{\left (t,R,Z \right )} \overline{\Omega_{1}{\left (t,R,Z \right )}} - \frac{i}{2} \rho_{32}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{33}{\left (t,R,Z \right )} & = & - \gamma_{32} \rho_{33}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{23}{\left (t,R,Z \right )} + \frac{i}{2} \rho_{32}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{21}{\left (t,R,Z \right )} & = & i \delta_{1} \rho_{21}{\left (t,R,Z \right )} - \frac{\gamma_{21}}{2} \rho_{21}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{11}{\left (t,R,Z \right )} + \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{22}{\left (t,R,Z \right )} - \frac{i}{2} \rho_{31}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{31}{\left (t,R,Z \right )} & = & i \delta_{1} \rho_{31}{\left (t,R,Z \right )} + i \delta_{2} \rho_{31}{\left (t,R,Z \right )} - \frac{\gamma_{32}}{2} \rho_{31}{\left (t,R,Z \right )} + \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{32}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{21}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t} \rho_{32}{\left (t,R,Z \right )} & = & i \delta_{2} \rho_{32}{\left (t,R,Z \right )} - \frac{\gamma_{21}}{2} \rho_{32}{\left (t,R,Z \right )} - \frac{\gamma_{32}}{2} \rho_{32}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{22}{\left (t,R,Z \right )} + \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{33}{\left (t,R,Z \right )} + \frac{i}{2} \rho_{31}{\left (t,R,Z \right )} \overline{\Omega_{1}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \Omega_{1}{\left (t,R,Z \right )} & = & - \frac{i \mu_{0}}{\hbar} c^{2} e^{2} r_{+1;21}^{2} \varpi_{1} n{\left (R,Z \right )} \rho_{21}{\left (t,R,Z \right )} - c \frac{\partial}{\partial Z} \Omega_{1}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t} \Omega_{2}{\left (t,R,Z \right )} & = & - \frac{i \mu_{0}}{\hbar} c^{2} e^{2} r_{-1;32}^{2} \varpi_{2} n{\left (R,Z \right )} \rho_{32}{\left (t,R,Z \right )} + c \frac{\partial}{\partial Z} \Omega_{2}{\left (t,R,Z \right )}\end{matrix}\right]$$

Question: how to use calculations from detunings as a replacement for calculations for velocity classes.

Linear approximation

We can linearize this equations by approximating all the population to be in $\rho_{11}$, and taking field $\Omega_1$ to be weak, and the coherences to be small. We will add an $\epsilon$ to all terms involving either coherences or signal field.



In [89]:

    
epsilon=symbols("varepsilon", real=True)
lin_subs={rho[0,0]:1, rho[1,1]:epsilon**2*rho[1,1], rho[2,2]:epsilon**2*rho[2,2],
          rho[1,0]:epsilon*rho[1,0], rho[0,1]:epsilon*rho[0,1],
          rho[2,0]:epsilon*rho[2,0], rho[0,2]:epsilon*rho[0,2],
          rho[2,1]:epsilon*rho[2,1], rho[1,2]:epsilon*rho[1,2],
          Omega1: epsilon*Omega1}
# lin_subs



In [90]:

    
eqs_lin=eqs_list
eqs_lin=eqs_lin.subs(lin_subs)
fprint(eqs_lin, print_ascii=print_ascii)









    Out[90]:





$$\left[\begin{matrix}\frac{d}{d t} 1 & = & \gamma_{21} \varepsilon^{2} \rho_{22}{\left (t,R,Z \right )} + \frac{i \varepsilon^{2}}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{12}{\left (t,R,Z \right )} - \frac{i \varepsilon^{2}}{2} \rho_{21}{\left (t,R,Z \right )} \overline{\Omega_{1}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t}\left(\varepsilon^{2} \rho_{22}{\left (t,R,Z \right )}\right) & = & - \gamma_{21} \varepsilon^{2} \rho_{22}{\left (t,R,Z \right )} + \gamma_{32} \varepsilon^{2} \rho_{33}{\left (t,R,Z \right )} - \frac{i \varepsilon^{2}}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{12}{\left (t,R,Z \right )} + \frac{i \varepsilon^{2}}{2} \rho_{21}{\left (t,R,Z \right )} \overline{\Omega_{1}{\left (t,R,Z \right )}} + \frac{i \varepsilon}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{23}{\left (t,R,Z \right )} - \frac{i \varepsilon}{2} \rho_{32}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t}\left(\varepsilon^{2} \rho_{33}{\left (t,R,Z \right )}\right) & = & - \gamma_{32} \varepsilon^{2} \rho_{33}{\left (t,R,Z \right )} - \frac{i \varepsilon}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{23}{\left (t,R,Z \right )} + \frac{i \varepsilon}{2} \rho_{32}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t}\left(\varepsilon \rho_{21}{\left (t,R,Z \right )}\right) & = & i \delta_{1} \varepsilon \rho_{21}{\left (t,R,Z \right )} - \frac{\gamma_{21} \varepsilon}{2} \rho_{21}{\left (t,R,Z \right )} + \frac{i \varepsilon^{3}}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{22}{\left (t,R,Z \right )} - \frac{i \varepsilon}{2} \Omega_{1}{\left (t,R,Z \right )} - \frac{i \varepsilon}{2} \rho_{31}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t}\left(\varepsilon \rho_{31}{\left (t,R,Z \right )}\right) & = & i \delta_{1} \varepsilon \rho_{31}{\left (t,R,Z \right )} + i \delta_{2} \varepsilon \rho_{31}{\left (t,R,Z \right )} - \frac{\gamma_{32} \varepsilon}{2} \rho_{31}{\left (t,R,Z \right )} + \frac{i \varepsilon^{2}}{2} \Omega_{1}{\left (t,R,Z \right )} \rho_{32}{\left (t,R,Z \right )} - \frac{i \varepsilon}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{21}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t}\left(\varepsilon \rho_{32}{\left (t,R,Z \right )}\right) & = & i \delta_{2} \varepsilon \rho_{32}{\left (t,R,Z \right )} - \frac{\gamma_{21} \varepsilon}{2} \rho_{32}{\left (t,R,Z \right )} - \frac{\gamma_{32} \varepsilon}{2} \rho_{32}{\left (t,R,Z \right )} - \frac{i \varepsilon^{2}}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{22}{\left (t,R,Z \right )} + \frac{i \varepsilon^{2}}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{33}{\left (t,R,Z \right )} + \frac{i \varepsilon^{2}}{2} \rho_{31}{\left (t,R,Z \right )} \overline{\Omega_{1}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t}\left(\varepsilon \Omega_{1}{\left (t,R,Z \right )}\right) & = & - \frac{i \varepsilon}{\hbar} c^{2} e^{2} \mu_{0} r_{+1;21}^{2} \varpi_{1} n{\left (R,Z \right )} \rho_{21}{\left (t,R,Z \right )} - c \frac{\partial}{\partial Z}\left(\varepsilon \Omega_{1}{\left (t,R,Z \right )}\right)\\\frac{\partial}{\partial t} \Omega_{2}{\left (t,R,Z \right )} & = & - \frac{i \varepsilon}{\hbar} c^{2} e^{2} \mu_{0} r_{-1;32}^{2} \varpi_{2} n{\left (R,Z \right )} \rho_{32}{\left (t,R,Z \right )} + c \frac{\partial}{\partial Z} \Omega_{2}{\left (t,R,Z \right )}\end{matrix}\right]$$

And so, terms which are doubly small can be neglected.



In [91]:

    
eqs_lin=eqs_lin.subs({epsilon**2:0}).subs({epsilon:1})
eqs_lin=eqs_lin.subs({Derivative(1,t):0,Derivative(0,t):0})
fprint(eqs_lin, print_ascii=print_ascii)









    Out[91]:





$$\left[\begin{matrix}0 & = & 0\\0 & = & \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{23}{\left (t,R,Z \right )} - \frac{i}{2} \rho_{32}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\0 & = & - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{23}{\left (t,R,Z \right )} + \frac{i}{2} \rho_{32}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{21}{\left (t,R,Z \right )} & = & i \delta_{1} \rho_{21}{\left (t,R,Z \right )} - \frac{\gamma_{21}}{2} \rho_{21}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} - \frac{i}{2} \rho_{31}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{31}{\left (t,R,Z \right )} & = & i \delta_{1} \rho_{31}{\left (t,R,Z \right )} + i \delta_{2} \rho_{31}{\left (t,R,Z \right )} - \frac{\gamma_{32}}{2} \rho_{31}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{21}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t} \rho_{32}{\left (t,R,Z \right )} & = & i \delta_{2} \rho_{32}{\left (t,R,Z \right )} - \frac{\gamma_{21}}{2} \rho_{32}{\left (t,R,Z \right )} - \frac{\gamma_{32}}{2} \rho_{32}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t} \Omega_{1}{\left (t,R,Z \right )} & = & - \frac{i \mu_{0}}{\hbar} c^{2} e^{2} r_{+1;21}^{2} \varpi_{1} n{\left (R,Z \right )} \rho_{21}{\left (t,R,Z \right )} - c \frac{\partial}{\partial Z} \Omega_{1}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t} \Omega_{2}{\left (t,R,Z \right )} & = & - \frac{i \mu_{0}}{\hbar} c^{2} e^{2} r_{-1;32}^{2} \varpi_{2} n{\left (R,Z \right )} \rho_{32}{\left (t,R,Z \right )} + c \frac{\partial}{\partial Z} \Omega_{2}{\left (t,R,Z \right )}\end{matrix}\right]$$

The equations for $\rho_{22}, \rho_{33}$ make a statement about the imaginary part of $\Omega_2\rho_{23}$ being zero. Let's just ignore that. The equation for $\rho_{32}$ is only coupled to itself, and can be easily solved



In [92]:

    
rho32_sol=exp((-(gamma[1,0]+gamma[2,1])/2+I*delta2)*t)
rho32_sol









    Out[92]:





$$e^{t \left(i \delta_{2} - \frac{\gamma_{21}}{2} - \frac{\gamma_{32}}{2}\right)}$$



In [93]:

    
Derivative(rho32_sol,t).doit().expand()-eqs_lin[5,2].subs({rho[2,1]:rho32_sol}).expand()









    Out[93]:





$$0$$

So the dynamics are captured by the equations for $\rho_{21}, \rho_{31}$ and the wave propagation equations.



In [94]:

    
eqss=[[eqs_lin[3:5,:][0,0],eqsign,eqs_lin[3:5,:][0,2]],
      [eqs_lin[3:5,:][1,0],eqsign,eqs_lin[3:5,:][1,2]]]

eqss=eqss+[[eqs_lin[6: ,:][0,0],eqsign,eqs_lin[6: ,:][0,2]],
      [eqs_lin[6: ,:][1,0],eqsign,eqs_lin[6: ,:][1,2]]]

eqss=Matrix(eqss)
eqs_lin=eqss
eqs_lin









    Out[94]:





$$\left[\begin{matrix}\frac{\partial}{\partial t} \rho_{21}{\left (t,R,Z \right )} & = & i \delta_{1} \rho_{21}{\left (t,R,Z \right )} - \frac{\gamma_{21}}{2} \rho_{21}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} - \frac{i}{2} \rho_{31}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{31}{\left (t,R,Z \right )} & = & i \delta_{1} \rho_{31}{\left (t,R,Z \right )} + i \delta_{2} \rho_{31}{\left (t,R,Z \right )} - \frac{\gamma_{32}}{2} \rho_{31}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{21}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t} \Omega_{1}{\left (t,R,Z \right )} & = & - \frac{i \mu_{0}}{\hbar} c^{2} e^{2} r_{+1;21}^{2} \varpi_{1} n{\left (R,Z \right )} \rho_{21}{\left (t,R,Z \right )} - c \frac{\partial}{\partial Z} \Omega_{1}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t} \Omega_{2}{\left (t,R,Z \right )} & = & - \frac{i \mu_{0}}{\hbar} c^{2} e^{2} r_{-1;32}^{2} \varpi_{2} n{\left (R,Z \right )} \rho_{32}{\left (t,R,Z \right )} + c \frac{\partial}{\partial Z} \Omega_{2}{\left (t,R,Z \right )}\end{matrix}\right]$$

If we further consider $\Omega_2$ to be very strong we can neglect the first term on the right-hand side on the last equation.



In [95]:

    
eqOmega2=eqs_lin[3,:][0]-eqs_lin[3,:][2]
eqOmega2=eqOmega2.subs({n:0})
eqOmega2









    Out[95]:





$$- c \frac{\partial}{\partial Z} \Omega_{2}{\left (t,R,Z \right )} + \frac{\partial}{\partial t} \Omega_{2}{\left (t,R,Z \right )}$$

And thus solve that equation and take $\Omega_2$ as a given quantity. For instance, it might be a gaussian pulse with time duration $\tau_2$:



In [96]:

    
from sympy import log
Omega2_peak, tau2, t0=symbols("Omega2_peak, tau2, t0", positive=True)
Omega2_sol=Omega2_peak*exp(-4*log(2)*(t-t0+Z/c)/tau2**2)
Omega2_sol









    Out[96]:





$$\Omega_{2 peak} e^{- \frac{4}{\tau_{2}^{2}} \left(\frac{Z}{c} + t - t_{0}\right) \log{\left (2 \right )}}$$



In [97]:

    
eqOmega2.subs({Omega2:Omega2_sol}).doit()









    Out[97]:





$$0$$

Then the equations are reduced to



In [98]:

    
eqs_lin=eqs_lin[:-1,:].subs({c**2: 1/mu0/epsilon0})
eqs_lin









    Out[98]:





$$\left[\begin{matrix}\frac{\partial}{\partial t} \rho_{21}{\left (t,R,Z \right )} & = & i \delta_{1} \rho_{21}{\left (t,R,Z \right )} - \frac{\gamma_{21}}{2} \rho_{21}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} - \frac{i}{2} \rho_{31}{\left (t,R,Z \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\\\frac{\partial}{\partial t} \rho_{31}{\left (t,R,Z \right )} & = & i \delta_{1} \rho_{31}{\left (t,R,Z \right )} + i \delta_{2} \rho_{31}{\left (t,R,Z \right )} - \frac{\gamma_{32}}{2} \rho_{31}{\left (t,R,Z \right )} - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{21}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t} \Omega_{1}{\left (t,R,Z \right )} & = & - c \frac{\partial}{\partial Z} \Omega_{1}{\left (t,R,Z \right )} - \frac{i e^{2} r_{+1;21}^{2} \varpi_{1}}{\hbar \varepsilon_{0}} n{\left (R,Z \right )} \rho_{21}{\left (t,R,Z \right )}\end{matrix}\right]$$

The Doppler effect

The detunings of the fields however, will be different for different velocity classes in an atomic vapour.



In [99]:

    
vX, vY, vZ=symbols("v_X, v_Y, v_Z", real=True)
v = Matrix([vX, vY, vZ])
# v

For velocities $v \ll c$ the optical frequencies will be

\begin{equation} \varpi_l' = \varpi_l - \vec{k}_l \cdot \vec{v}_l \end{equation}

And so, in our case \begin{equation} \delta_1' = \delta_1 - \varpi_1 v_Z/c, \hspace{1cm} \delta_2' = \delta_2 + \varpi_2 v_Z/c \end{equation}

And the density matrix will be a mixture of contributions from all velocity classes \begin{equation} \rho = \int dv_Z g(v_Z) \rho{v_Z}' \end{equation}

So we will rewrite our equations in terms of this new $v_Z$



In [100]:

    
rhop=define_density_matrix(Ne, variables=[t,Z,vZ])
g=Function("g")(vZ)
g









    Out[100]:





$$g{\left (v_{Z} \right )}$$



In [101]:

    
delta1p=delta1-omega_laser[0]*vZ/c
delta2p=delta2+omega_laser[1]*vZ/c
delta1p,delta2p









    Out[101]:





$$\left ( \delta_{1} - \frac{v_{Z} \varpi_{1}}{c}, \quad \delta_{2} + \frac{v_{Z} \varpi_{2}}{c}\right )$$

So our final equations are:



In [102]:

    
eqs_linv=eqs_lin.subs({delta1: delta1p, delta2: delta2p, rho[1,0]: rhop[1,0], rho[2,0]: rhop[2,0]})
eqs_linv[2,:]= eqs_lin[2,:].subs({rho[1,0]: Integral(g*rhop[1,0],(vZ,-oo,oo))})
eqs_linv=eqs_linv.expand()
eqs_linv









    Out[102]:





$$\left[\begin{matrix}\frac{\partial}{\partial t} \rho_{21}{\left (t,Z,v_{Z} \right )} & = & i \delta_{1} \rho_{21}{\left (t,Z,v_{Z} \right )} - \frac{\gamma_{21}}{2} \rho_{21}{\left (t,Z,v_{Z} \right )} - \frac{i}{2} \Omega_{1}{\left (t,R,Z \right )} - \frac{i}{2} \rho_{31}{\left (t,Z,v_{Z} \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}} - \frac{i v_{Z}}{c} \varpi_{1} \rho_{21}{\left (t,Z,v_{Z} \right )}\\\frac{\partial}{\partial t} \rho_{31}{\left (t,Z,v_{Z} \right )} & = & i \delta_{1} \rho_{31}{\left (t,Z,v_{Z} \right )} + i \delta_{2} \rho_{31}{\left (t,Z,v_{Z} \right )} - \frac{\gamma_{32}}{2} \rho_{31}{\left (t,Z,v_{Z} \right )} - \frac{i}{2} \Omega_{2}{\left (t,R,Z \right )} \rho_{21}{\left (t,Z,v_{Z} \right )} - \frac{i v_{Z}}{c} \varpi_{1} \rho_{31}{\left (t,Z,v_{Z} \right )} + \frac{i v_{Z}}{c} \varpi_{2} \rho_{31}{\left (t,Z,v_{Z} \right )}\\\frac{\partial}{\partial t} \Omega_{1}{\left (t,R,Z \right )} & = & - c \frac{\partial}{\partial Z} \Omega_{1}{\left (t,R,Z \right )} - \frac{i e^{2} r_{+1;21}^{2} \varpi_{1}}{\hbar \varepsilon_{0}} n{\left (R,Z \right )} \int_{-\infty}^{\infty} g{\left (v_{Z} \right )} \rho_{21}{\left (t,Z,v_{Z} \right )}\, dv_{Z}\end{matrix}\right]$$

Adiabatic approximation



In [103]:

    
rho21_ad = solve(eqs_linv[0, 2], rhop[1, 0])[0]
# rho21_ad



In [104]:

    
Gamma21 = symbols("Gamma_21")
Gamma32 = symbols("Gamma_32")
#Gamma/2 = gamma[1, 0]/2 - I*(delta1-vZ*omega_laser[0]/c)
ss_Gamma21 = {gamma[1, 0]: 2*(I*(delta1-vZ*omega_laser[0]/c)-Gamma21/2)}
ss_Gamma32 = {gamma[2, 1]: 2*(I*(delta2+vZ*omega_laser[1]/c)-Gamma32/2)}

ss_Gamma21_inv = {Gamma21: gamma[1, 0]- 2*I*(delta1-vZ*omega_laser[0]/c)}
ss_Gamma32_inv = {Gamma32: gamma[2, 1]- 2*I*(delta2+vZ*omega_laser[1]/c)}

ss_Gamma21, ss_Gamma21_inv









    Out[104]:





$$\left ( \left \{ \gamma_{21} : - \Gamma_{21} + 2 i \left(\delta_{1} - \frac{v_{Z} \varpi_{1}}{c}\right)\right \}, \quad \left \{ \Gamma_{21} : \gamma_{21} - 2 i \left(\delta_{1} - \frac{v_{Z} \varpi_{1}}{c}\right)\right \}\right )$$



In [105]:

    
rho21_ad = simplify(rho21_ad.subs(ss_Gamma21))
ss_rho21_ad = {rhop[1, 0]: rho21_ad}
ss_rho21_ad









    Out[105]:





$$\left \{ \rho_{21}{\left (t,Z,v_{Z} \right )} : \frac{i}{\Gamma_{21}} \left(\Omega_{1}{\left (t,R,Z \right )} + \rho_{31}{\left (t,Z,v_{Z} \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\right)\right \}$$



In [106]:

    
eq1 = eqs_linv[1, 2].subs(ss_rho21_ad)#.subs(ss_Gamma32).expand()
# eq1

We set this in two-photon resonance



In [107]:

    
eq1 = eq1.subs(delta1, -delta2).expand()
# eq1



In [108]:

    
eq2 = eqs_lin[2,2].subs(rho[1, 0], rhop[1, 0]).subs(ss_rho21_ad).expand()
# eq2



In [109]:

    
def group(expr, variables, return_coeffs=False):
    coefs = [Derivative(expr, variables[i]).doit() for i in range(len(variables))]
    expr_lin = sum([coefs[i]*variables[i] for i in range(len(variables))])
    remainder = (expr - expr_lin).expand()
    if return_coeffs:
        return coefs
    return expr_lin + remainder
variables = [rhop[2, 0], Omega1]
eq1 = group(eq1, variables)
eq2 = group(eq2, variables)

So the equations in the adiabatic approximation are



In [110]:

    
eqs_linva = Matrix([[Derivative(rhop[2, 0], t), eqsign, eq1],
                    [Derivative(Omega1, t), eqsign, eq2]])
eqs_linva









    Out[110]:





$$\left[\begin{matrix}\frac{\partial}{\partial t} \rho_{31}{\left (t,Z,v_{Z} \right )} & = & \left(- \frac{\gamma_{32}}{2} - \frac{i v_{Z}}{c} \varpi_{1} + \frac{i v_{Z}}{c} \varpi_{2} + \frac{\Omega_{2}{\left (t,R,Z \right )}}{2 \Gamma_{21}} \overline{\Omega_{2}{\left (t,R,Z \right )}}\right) \rho_{31}{\left (t,Z,v_{Z} \right )} + \frac{\Omega_{1}{\left (t,R,Z \right )}}{2 \Gamma_{21}} \Omega_{2}{\left (t,R,Z \right )}\\\frac{\partial}{\partial t} \Omega_{1}{\left (t,R,Z \right )} & = & - c \frac{\partial}{\partial Z} \Omega_{1}{\left (t,R,Z \right )} + \frac{e^{2} r_{+1;21}^{2} \varpi_{1}}{\Gamma_{21} \hbar \varepsilon_{0}} \Omega_{1}{\left (t,R,Z \right )} n{\left (R,Z \right )} + \frac{e^{2} r_{+1;21}^{2} \varpi_{1}}{\Gamma_{21} \hbar \varepsilon_{0}} n{\left (R,Z \right )} \rho_{31}{\left (t,Z,v_{Z} \right )} \overline{\Omega_{2}{\left (t,R,Z \right )}}\end{matrix}\right]$$

We transform into the control field reference frame.



In [111]:

    
tau = symbols("tau", real=True)

eqs_linva2 = eqs_linva.subs(t, tau)

eqs_linva2[1, 2] = eqs_linva2[1, 2]/2

Omega2_p = Function("Omega2")(tau)
n_p = symbols("n", positive=True)
eqs_linva2 = eqs_linva2.subs({Omega2.subs(t, tau): Omega2_p})
eqs_linva2 = eqs_linva2.subs(n, n_p)
eqs_linva2









    Out[111]:





$$\left[\begin{matrix}\frac{\partial}{\partial \tau} \rho_{31}{\left (\tau,Z,v_{Z} \right )} & = & \left(- \frac{\gamma_{32}}{2} - \frac{i v_{Z}}{c} \varpi_{1} + \frac{i v_{Z}}{c} \varpi_{2} + \frac{\Omega_{2}{\left (\tau \right )}}{2 \Gamma_{21}} \overline{\Omega_{2}{\left (\tau \right )}}\right) \rho_{31}{\left (\tau,Z,v_{Z} \right )} + \frac{\Omega_{1}{\left (\tau,R,Z \right )}}{2 \Gamma_{21}} \Omega_{2}{\left (\tau \right )}\\\frac{\partial}{\partial \tau} \Omega_{1}{\left (\tau,R,Z \right )} & = & - \frac{c}{2} \frac{\partial}{\partial Z} \Omega_{1}{\left (\tau,R,Z \right )} + \frac{e^{2} n r_{+1;21}^{2} \varpi_{1} \Omega_{1}{\left (\tau,R,Z \right )}}{2 \Gamma_{21} \hbar \varepsilon_{0}} + \frac{e^{2} n r_{+1;21}^{2} \varpi_{1} \overline{\Omega_{2}{\left (\tau \right )}}}{2 \Gamma_{21} \hbar \varepsilon_{0}} \rho_{31}{\left (\tau,Z,v_{Z} \right )}\end{matrix}\right]$$

We can then take a fourier transform of the Z dependence



In [112]:

    
kz = symbols("k_z", real=True)
Omega1_hat = Function(r"\hat{\Omega}_1")(tau, kz)
rho31_hat = Function(r"\hat{\rho}_{31}")(tau, kz, vZ)
# Omega1_hat, rho31_hat



In [113]:

    
eqs_linva_hat = eqs_linva2

eqs_linva_hat = eqs_linva_hat.subs(Derivative(Omega1, Z).subs(t, tau), I*kz*Omega1_hat)
eqs_linva_hat = eqs_linva_hat.subs(Omega1.subs(t, tau), Omega1_hat)
eqs_linva_hat = eqs_linva_hat.subs(rhop[2, 0].subs(t, tau), rho31_hat)
eqs_linva_hat









    Out[113]:





$$\left[\begin{matrix}\frac{\partial}{\partial \tau} \hat{\rho}_{31}{\left (\tau,k_{z},v_{Z} \right )} & = & \left(- \frac{\gamma_{32}}{2} - \frac{i v_{Z}}{c} \varpi_{1} + \frac{i v_{Z}}{c} \varpi_{2} + \frac{\Omega_{2}{\left (\tau \right )}}{2 \Gamma_{21}} \overline{\Omega_{2}{\left (\tau \right )}}\right) \hat{\rho}_{31}{\left (\tau,k_{z},v_{Z} \right )} + \frac{\Omega_{2}{\left (\tau \right )}}{2 \Gamma_{21}} \hat{\Omega}_1{\left (\tau,k_{z} \right )}\\\frac{\partial}{\partial \tau} \hat{\Omega}_1{\left (\tau,k_{z} \right )} & = & - \frac{i k_{z}}{2} c \hat{\Omega}_1{\left (\tau,k_{z} \right )} + \frac{e^{2} n r_{+1;21}^{2} \varpi_{1} \hat{\Omega}_1{\left (\tau,k_{z} \right )}}{2 \Gamma_{21} \hbar \varepsilon_{0}} + \frac{e^{2} n r_{+1;21}^{2} \varpi_{1} \hat{\rho}_{31}{\left (\tau,k_{z},v_{Z} \right )}}{2 \Gamma_{21} \hbar \varepsilon_{0}} \overline{\Omega_{2}{\left (\tau \right )}}\end{matrix}\right]$$



In [114]:

    
M = Matrix([group(eqs_linva_hat[i, 2], [rho31_hat, Omega1_hat], return_coeffs=True) for i in range(2)])/I
M = M.expand()
M









    Out[114]:





$$\left[\begin{matrix}\frac{i \gamma_{32}}{2} - \frac{v_{Z} \varpi_{1}}{c} + \frac{v_{Z} \varpi_{2}}{c} - \frac{i}{2 \Gamma_{21}} \Omega_{2}{\left (\tau \right )} \overline{\Omega_{2}{\left (\tau \right )}} & - \frac{i}{2 \Gamma_{21}} \Omega_{2}{\left (\tau \right )}\\- \frac{i e^{2} n r_{+1;21}^{2} \varpi_{1}}{2 \Gamma_{21} \hbar \varepsilon_{0}} \overline{\Omega_{2}{\left (\tau \right )}} & - \frac{c k_{z}}{2} - \frac{i e^{2} n r_{+1;21}^{2} \varpi_{1}}{2 \Gamma_{21} \hbar \varepsilon_{0}}\end{matrix}\right]$$



In [115]:

    
S = Function("S")(tau, kz)
xx = M[1, 0]/Omega2_p.conjugate()*Gamma21/-I*2
ss_rescale = {Omega1_hat: S*sqrt(xx)}
ss_rescale_inv = {S: Omega1_hat/sqrt(xx)}
# xx, ss_rescale



In [116]:

    
eq1 = eqs_linva_hat[0, 0]-eqs_linva_hat[0, 2]
eq1 = eq1.subs(ss_rescale)
eq1 = eq1 - Derivative(rho31_hat, tau)
# eq1



In [117]:

    
eq2 = eqs_linva_hat[1, 0]-eqs_linva_hat[1, 2]
#eq2 = eq2.subs(ss_rescale).doit()
eq2 = (eq2.subs(ss_rescale).doit()/sqrt(xx)).expand()
eq2 = eq2 - Derivative(S, tau)
#eq2 = eq2.conjugate()
# eq2



In [118]:

    
kappa = symbols("kappa", positive=True)
Deltak = symbols(r"\Delta_k", real=True)
rr = -r_p[0][0, 1]
ss_kappa = {rr: kappa*sqrt(hbar*epsilon0*c/(n_p*omega_laser[0]*e**2))}
ss_kappa_inv = {kappa: rr/sqrt(hbar*epsilon0*c/(n_p*omega_laser[0]*e**2))}

ss_Dop = {vZ*omega_laser[0]/c-vZ*omega_laser[1]/c: vZ*Deltak}
ss_Dop_inv = {vZ*Deltak: vZ*omega_laser[0]/c-vZ*omega_laser[1]/c}

ss_kappa_inv, ss_Gamma21_inv, ss_rescale_inv, ss_Dop_inv









    Out[118]:





$$\left ( \left \{ \kappa : \frac{e \sqrt{n} r_{+1;21} \sqrt{\varpi_{1}}}{\sqrt{c} \sqrt{\hbar} \sqrt{\varepsilon_{0}}}\right \}, \quad \left \{ \Gamma_{21} : \gamma_{21} - 2 i \left(\delta_{1} - \frac{v_{Z} \varpi_{1}}{c}\right)\right \}, \quad \left \{ S{\left (\tau,k_{z} \right )} : \frac{\sqrt{\hbar} \sqrt{\varepsilon_{0}} \hat{\Omega}_1{\left (\tau,k_{z} \right )}}{e \sqrt{n} \sqrt{\varpi_{1}} \left|{r_{+1;21}}\right|}\right \}, \quad \left \{ \Delta_k v_{Z} : \frac{v_{Z} \varpi_{1}}{c} - \frac{v_{Z} \varpi_{2}}{c}\right \}\right )$$



In [119]:

    
Y = Matrix([rho31_hat, S])
# Y



In [120]:

    
eqss = [eq2, eq1]
M = Matrix([group(eqss[i], [S, rho31_hat], return_coeffs=True) for i in range(2)])/I
M = M.expand().subs(ss_kappa).subs(ss_Dop)
M









    Out[120]:





$$\left[\begin{matrix}\frac{c k_{z}}{2} + \frac{i c \kappa^{2}}{2 \Gamma_{21}} & \frac{i \sqrt{c} \kappa}{2 \Gamma_{21}} \overline{\Omega_{2}{\left (\tau \right )}}\\\frac{i \sqrt{c} \kappa}{2 \Gamma_{21}} \Omega_{2}{\left (\tau \right )} & \Delta_k v_{Z} - \frac{i \gamma_{32}}{2} + \frac{i}{2 \Gamma_{21}} \Omega_{2}{\left (\tau \right )} \overline{\Omega_{2}{\left (\tau \right )}}\end{matrix}\right]$$

Solution for square pulses.

Let's assume that the control field is constant for a duration $T$.



In [121]:

    
T, Omega2a = symbols("T, Omega_2a", positive=True)
phi = symbols("phi", real=True)
T, Omega2a, phi









    Out[121]:





$$\left ( T, \quad \Omega_{2a}, \quad \phi\right )$$



In [122]:

    
ss_square = {Omega2_p: Omega2a*exp(I*phi)}
M_square = M.subs(ss_square)
M_square









    Out[122]:





$$\left[\begin{matrix}\frac{c k_{z}}{2} + \frac{i c \kappa^{2}}{2 \Gamma_{21}} & \frac{i \Omega_{2a} \sqrt{c} \kappa}{2 \Gamma_{21}} e^{- i \phi}\\\frac{i \Omega_{2a} \sqrt{c} \kappa}{2 \Gamma_{21}} e^{i \phi} & \Delta_k v_{Z} - \frac{i \gamma_{32}}{2} + \frac{i \Omega_{2a}^{2}}{2 \Gamma_{21}}\end{matrix}\right]$$



In [123]:

    
p, q, r, s = symbols("p, q, r, s")
ss_M = {M_square[0, 0]: p, M_square[0, 1]: q, M_square[1, 0]: r, M_square[1, 1]: s}
ss_M_inv = {ss_M[i]: i for i in ss_M}
ss_M_inv









    Out[123]:





$$\left \{ p : \frac{c k_{z}}{2} + \frac{i c \kappa^{2}}{2 \Gamma_{21}}, \quad q : \frac{i \Omega_{2a} \sqrt{c} \kappa}{2 \Gamma_{21}} e^{- i \phi}, \quad r : \frac{i \Omega_{2a} \sqrt{c} \kappa}{2 \Gamma_{21}} e^{i \phi}, \quad s : \Delta_k v_{Z} - \frac{i \gamma_{32}}{2} + \frac{i \Omega_{2a}^{2}}{2 \Gamma_{21}}\right \}$$



In [124]:

    
print ss_M_inv









    



{s: \Delta_k*v_Z - I*gamma_32/2 + I*Omega_2a**2/(2*Gamma_21), p: c*k_z/2 + I*c*kappa**2/(2*Gamma_21), r: I*Omega_2a*sqrt(c)*kappa*exp(I*phi)/(2*Gamma_21), q: I*Omega_2a*sqrt(c)*kappa*exp(-I*phi)/(2*Gamma_21)}



In [125]:

    
M_simp = M_square.subs(ss_M)
M_simp









    Out[125]:





$$\left[\begin{matrix}p & q\\r & s\end{matrix}\right]$$

We have equations

$$ \dot{X} = i M X $$



In [126]:

    
sig = (s+p)/2
delt = (s-p)/2
alp = sig**2 - M_simp.det()
lamp, lamm = sig + sqrt(alp), sig - sqrt(alp)
lamp, lamm









    Out[126]:





$$\left ( \frac{p}{2} + \frac{s}{2} + \sqrt{- p s + q r + \left(\frac{p}{2} + \frac{s}{2}\right)^{2}}, \quad \frac{p}{2} + \frac{s}{2} - \sqrt{- p s + q r + \left(\frac{p}{2} + \frac{s}{2}\right)^{2}}\right )$$



In [127]:

    
xp, xm = q, q
yp, ym = delt+sqrt(alp), delt-sqrt(alp)



In [128]:

    
delt.subs(ss_M_inv)









    Out[128]:





$$\frac{\Delta_k v_{Z}}{2} - \frac{c k_{z}}{4} - \frac{i \gamma_{32}}{4} + \frac{i \Omega_{2a}^{2}}{4 \Gamma_{21}} - \frac{i c \kappa^{2}}{4 \Gamma_{21}}$$



In [129]:

    
M_square









    Out[129]:





$$\left[\begin{matrix}\frac{c k_{z}}{2} + \frac{i c \kappa^{2}}{2 \Gamma_{21}} & \frac{i \Omega_{2a} \sqrt{c} \kappa}{2 \Gamma_{21}} e^{- i \phi}\\\frac{i \Omega_{2a} \sqrt{c} \kappa}{2 \Gamma_{21}} e^{i \phi} & \Delta_k v_{Z} - \frac{i \gamma_{32}}{2} + \frac{i \Omega_{2a}^{2}}{2 \Gamma_{21}}\end{matrix}\right]$$



In [130]:

    
A = Matrix([[xp, xm], [yp, ym]])
Ainv = Matrix([[ym, -xm], [-yp, xp]])/(xp*ym-xm*yp)
D = Matrix([[lamp, 0], [0, lamm]])



In [131]:

    
A









    Out[131]:





$$\left[\begin{matrix}q & q\\- \frac{p}{2} + \frac{s}{2} + \sqrt{- p s + q r + \left(\frac{p}{2} + \frac{s}{2}\right)^{2}} & - \frac{p}{2} + \frac{s}{2} - \sqrt{- p s + q r + \left(\frac{p}{2} + \frac{s}{2}\right)^{2}}\end{matrix}\right]$$



In [132]:

    
simplify(A*D*Ainv)









    Out[132]:





$$\left[\begin{matrix}p & q\\r & s\end{matrix}\right]$$



In [133]:

    
x_p, x_m, y_p, y_m, lam_p, lam_m = symbols("x_+, x_-, y_+, y_-, \lambda_+, \lambda_-")
x_p, x_m, y_p, y_m, lam_p, lam_m









    Out[133]:





$$\left ( x_{+}, \quad x_{-}, \quad y_{+}, \quad y_{-}, \quad \lambda_+, \quad \lambda_-\right )$$



In [134]:

    
ss_xy = {xp: x_p, yp: y_p, ym: y_m, lamp: lam_p, lamm: lam_m}
ss_xy_inv = {x_p: xp, x_m: xm, y_p: yp, y_m: ym, lam_p: lamp, lam_m: lamm}
# ss_xy



In [135]:

    
A = Matrix([[x_p, x_m], [y_p, y_m]])
Ainv = Matrix([[y_m, -x_m], [-y_p, x_p]])/(x_p*y_m-x_m*y_p)
D = Matrix([[lam_p, 0], [0, lam_m]])
A, Ainv









    Out[135]:





$$\left ( \left[\begin{matrix}x_{+} & x_{-}\\y_{+} & y_{-}\end{matrix}\right], \quad \left[\begin{matrix}\frac{y_{-}}{x_{+} y_{-} - x_{-} y_{+}} & - \frac{x_{-}}{x_{+} y_{-} - x_{-} y_{+}}\\- \frac{y_{+}}{x_{+} y_{-} - x_{-} y_{+}} & \frac{x_{+}}{x_{+} y_{-} - x_{-} y_{+}}\end{matrix}\right]\right )$$



In [136]:

    
from sympy import eye
G = sum([A[:, i]*Ainv[i, :]*exp(I*T*D[i, i]) for i in range(2)], zeros(2, 2))
G = simplify(G)
G









    Out[136]:





$$\left[\begin{matrix}\frac{1}{x_{+} y_{-} - x_{-} y_{+}} \left(x_{+} y_{-} e^{i T \lambda_+} - x_{-} y_{+} e^{i T \lambda_-}\right) & \frac{x_{+} x_{-} \left(- e^{i T \lambda_+} + e^{i T \lambda_-}\right)}{x_{+} y_{-} - x_{-} y_{+}}\\\frac{y_{+} y_{-} \left(e^{i T \lambda_+} - e^{i T \lambda_-}\right)}{x_{+} y_{-} - x_{-} y_{+}} & \frac{1}{x_{+} y_{-} - x_{-} y_{+}} \left(x_{+} y_{-} e^{i T \lambda_-} - x_{-} y_{+} e^{i T \lambda_+}\right)\end{matrix}\right]$$



In [137]:

    
S_in = symbols("S_in")
rho31_in = symbols("rho_31in")
Yin = Y.subs({rho31_hat: 0, S: S_in})
Yint = Y.subs({rho31_hat: rho31_in, S: 0})
Yin, Yint









    Out[137]:





$$\left ( \left[\begin{matrix}0\\S_{in}\end{matrix}\right], \quad \left[\begin{matrix}\rho_{31in}\\0\end{matrix}\right]\right )$$



In [138]:

    
G*Yin









    Out[138]:





$$\left[\begin{matrix}\frac{S_{in} x_{+} x_{-}}{x_{+} y_{-} - x_{-} y_{+}} \left(- e^{i T \lambda_+} + e^{i T \lambda_-}\right)\\\frac{S_{in}}{x_{+} y_{-} - x_{-} y_{+}} \left(x_{+} y_{-} e^{i T \lambda_-} - x_{-} y_{+} e^{i T \lambda_+}\right)\end{matrix}\right]$$



In [139]:

    
G*Yint









    Out[139]:





$$\left[\begin{matrix}\frac{\rho_{31in}}{x_{+} y_{-} - x_{-} y_{+}} \left(x_{+} y_{-} e^{i T \lambda_+} - x_{-} y_{+} e^{i T \lambda_-}\right)\\\frac{\rho_{31in} y_{+} y_{-} \left(e^{i T \lambda_+} - e^{i T \lambda_-}\right)}{x_{+} y_{-} - x_{-} y_{+}}\end{matrix}\right]$$



In [140]:

    
Yint = G*Yin
Yint = Matrix([Yint[0, 0], 0])
Yout = G*Yint
Yout









    Out[140]:





$$\left[\begin{matrix}\frac{S_{in} x_{+} x_{-}}{\left(x_{+} y_{-} - x_{-} y_{+}\right)^{2}} \left(x_{+} y_{-} e^{i T \lambda_+} - x_{-} y_{+} e^{i T \lambda_-}\right) \left(- e^{i T \lambda_+} + e^{i T \lambda_-}\right)\\\frac{S_{in} x_{+} x_{-} y_{+} y_{-}}{\left(x_{+} y_{-} - x_{-} y_{+}\right)^{2}} \left(- e^{i T \lambda_+} + e^{i T \lambda_-}\right) \left(e^{i T \lambda_+} - e^{i T \lambda_-}\right)\end{matrix}\right]$$



In [141]:

    
Sout = Yout[1, 0]#.expand()
Sout









    Out[141]:





$$\frac{S_{in} x_{+} x_{-} y_{+} y_{-}}{\left(x_{+} y_{-} - x_{-} y_{+}\right)^{2}} \left(- e^{i T \lambda_+} + e^{i T \lambda_-}\right) \left(e^{i T \lambda_+} - e^{i T \lambda_-}\right)$$



In [142]:

    
sigma, delta, alpha = symbols("sigma, delta, alpha")
ss_amazing1 = {sigma: (s+p)/2, delta: (s-p)/2, alpha: sig**2 - M_simp.det()}
print ss_amazing1
ss_amazing1









    



{delta: -p/2 + s/2, sigma: p/2 + s/2, alpha: -p*s + q*r + (p/2 + s/2)**2}






    Out[142]:





$$\left \{ \alpha : - p s + q r + \left(\frac{p}{2} + \frac{s}{2}\right)^{2}, \quad \delta : - \frac{p}{2} + \frac{s}{2}, \quad \sigma : \frac{p}{2} + \frac{s}{2}\right \}$$



In [143]:

    
ss_amazing2 = {lam_p: sigma+sqrt(alpha), lam_m: sigma-sqrt(alpha)}
ss_amazing2









    Out[143]:





$$\left \{ \lambda_+ : \sqrt{\alpha} + \sigma, \quad \lambda_- : - \sqrt{\alpha} + \sigma\right \}$$



In [144]:

    
ss_amazing3 = {x_p: q, x_m: q}
ss_amazing3









    Out[144]:





$$\left \{ x_{+} : q, \quad x_{-} : q\right \}$$



In [145]:

    
ss_amazing4 = {y_p: delta+sqrt(alpha), y_m: delta-sqrt(alpha)}
ss_amazing4









    Out[145]:





$$\left \{ y_{+} : \sqrt{\alpha} + \delta, \quad y_{-} : - \sqrt{\alpha} + \delta\right \}$$



In [146]:

    
Sout_ = Sout.subs(ss_amazing2).factor()
Sout_ = Sout_.subs({exp(I*T*sqrt(alpha))-exp(-I*T*sqrt(alpha)): 2*I*sin(T*sqrt(alpha))})
Sout_ = Sout_.subs(ss_amazing3).factor()
Sout_ = Sout_.subs(ss_amazing4)
Sout_ = Sout_.subs({(delta-sqrt(alpha))*(delta+sqrt(alpha)): delta**2-alpha})
Sout_









    Out[146]:





$$\frac{S_{in}}{\alpha} \left(- \alpha + \delta^{2}\right) e^{2 i T \sigma} \sin^{2}{\left (T \sqrt{\alpha} \right )}$$



In [147]:

    
Fk = Sout_/S_in
Fk









    Out[147]:





$$\frac{1}{\alpha} \left(- \alpha + \delta^{2}\right) e^{2 i T \sigma} \sin^{2}{\left (T \sqrt{\alpha} \right )}$$

We check whether the phase in the exponential factor is real or imaginary.



In [148]:

    
sigma_phase = 2*I*T*ss_amazing1[sigma].subs(ss_M_inv).subs(ss_Gamma21_inv)
sigma_phase









    Out[148]:





$$2 i T \left(\frac{i \Omega_{2a}^{2}}{4 \gamma_{21} - 8 i \left(\delta_{1} - \frac{v_{Z} \varpi_{1}}{c}\right)} + \frac{\Delta_k v_{Z}}{2} + \frac{c k_{z}}{4} + \frac{i c \kappa^{2}}{4 \gamma_{21} - 8 i \left(\delta_{1} - \frac{v_{Z} \varpi_{1}}{c}\right)} - \frac{i \gamma_{32}}{4}\right)$$

Unsurprisingly, decoherence goes to zero if the the decay frequencies are zero.



In [149]:

    
re(sigma_phase).subs({gamma[1, 0]: 0, gamma[2, 1]: 0})









    Out[149]:





$$0$$



In [150]:

    
im(sigma_phase).subs({gamma[1, 0]: 0, gamma[2, 1]: 0})









    Out[150]:





$$2 T \left(- \frac{\Omega_{2a}^{2} \left(2 \delta_{1} - \frac{2 v_{Z}}{c} \varpi_{1}\right)}{4 \left(- 2 \delta_{1} + \frac{2 v_{Z}}{c} \varpi_{1}\right)^{2}} + \frac{\Delta_k v_{Z}}{2} + \frac{c k_{z}}{4} - \frac{c \kappa^{2} \left(2 \delta_{1} - \frac{2 v_{Z}}{c} \varpi_{1}\right)}{4 \left(- 2 \delta_{1} + \frac{2 v_{Z}}{c} \varpi_{1}\right)^{2}}\right)$$



In [168]:

    
F_guillaume = Fk.subs(ss_amazing1).subs(ss_M_inv)#.factor()
F_guillaume









    Out[168]:





$$\frac{e^{2 i T \left(\frac{\Delta_k v_{Z}}{2} + \frac{c k_{z}}{4} - \frac{i \gamma_{32}}{4} + \frac{i \Omega_{2a}^{2}}{4 \Gamma_{21}} + \frac{i c \kappa^{2}}{4 \Gamma_{21}}\right)} \sin^{2}{\left (T \sqrt{- \left(\frac{c k_{z}}{2} + \frac{i c \kappa^{2}}{2 \Gamma_{21}}\right) \left(\Delta_k v_{Z} - \frac{i \gamma_{32}}{2} + \frac{i \Omega_{2a}^{2}}{2 \Gamma_{21}}\right) + \left(\frac{\Delta_k v_{Z}}{2} + \frac{c k_{z}}{4} - \frac{i \gamma_{32}}{4} + \frac{i \Omega_{2a}^{2}}{4 \Gamma_{21}} + \frac{i c \kappa^{2}}{4 \Gamma_{21}}\right)^{2} - \frac{\Omega_{2a}^{2} c \kappa^{2}}{4 \Gamma_{21}^{2}}} \right )}}{- \left(\frac{c k_{z}}{2} + \frac{i c \kappa^{2}}{2 \Gamma_{21}}\right) \left(\Delta_k v_{Z} - \frac{i \gamma_{32}}{2} + \frac{i \Omega_{2a}^{2}}{2 \Gamma_{21}}\right) + \left(\frac{\Delta_k v_{Z}}{2} + \frac{c k_{z}}{4} - \frac{i \gamma_{32}}{4} + \frac{i \Omega_{2a}^{2}}{4 \Gamma_{21}} + \frac{i c \kappa^{2}}{4 \Gamma_{21}}\right)^{2} - \frac{\Omega_{2a}^{2} c \kappa^{2}}{4 \Gamma_{21}^{2}}} \left(\left(\frac{c k_{z}}{2} + \frac{i c \kappa^{2}}{2 \Gamma_{21}}\right) \left(\Delta_k v_{Z} - \frac{i \gamma_{32}}{2} + \frac{i \Omega_{2a}^{2}}{2 \Gamma_{21}}\right) + \left(\frac{\Delta_k v_{Z}}{2} - \frac{c k_{z}}{4} - \frac{i \gamma_{32}}{4} + \frac{i \Omega_{2a}^{2}}{4 \Gamma_{21}} - \frac{i c \kappa^{2}}{4 \Gamma_{21}}\right)^{2} - \left(\frac{\Delta_k v_{Z}}{2} + \frac{c k_{z}}{4} - \frac{i \gamma_{32}}{4} + \frac{i \Omega_{2a}^{2}}{4 \Gamma_{21}} + \frac{i c \kappa^{2}}{4 \Gamma_{21}}\right)^{2} + \frac{\Omega_{2a}^{2} c \kappa^{2}}{4 \Gamma_{21}^{2}}\right)$$



In [214]:

    
arg = F_guillaume.args[1].args[0].args[0]**2
arg = arg.factor().expand()
arg









    Out[214]:





$$\frac{\Delta_k^{2} v_{Z}^{2}}{4} T^{2} - \frac{\Delta_k k_{z}}{4} T^{2} c v_{Z} - \frac{i \Delta_k}{4} T^{2} \gamma_{32} v_{Z} + \frac{T^{2} k_{z}^{2}}{16} c^{2} + \frac{i \gamma_{32}}{8} T^{2} c k_{z} - \frac{T^{2} \gamma_{32}^{2}}{16} + \frac{i \Omega_{2a}^{2} \Delta_k v_{Z}}{4 \Gamma_{21}} T^{2} - \frac{i \Omega_{2a}^{2} c k_{z}}{8 \Gamma_{21}} T^{2} + \frac{\Omega_{2a}^{2} T^{2} \gamma_{32}}{8 \Gamma_{21}} - \frac{i \Delta_k c v_{Z}}{4 \Gamma_{21}} T^{2} \kappa^{2} + \frac{i T^{2} c^{2} k_{z}}{8 \Gamma_{21}} \kappa^{2} - \frac{T^{2} c \gamma_{32} \kappa^{2}}{8 \Gamma_{21}} - \frac{\Omega_{2a}^{4} T^{2}}{16 \Gamma_{21}^{2}} - \frac{\Omega_{2a}^{2} T^{2} c \kappa^{2}}{8 \Gamma_{21}^{2}} - \frac{T^{2} c^{2} \kappa^{4}}{16 \Gamma_{21}^{2}}$$



In [217]:

    
coefs = coeffs(arg, kz, 2)
#[coefs[i]*kz**i for i in reversed(range(3))]
coefs









    



T**2*\Delta_k**2*v_Z**2/4 - T**2*\Delta_k*c*k_z*v_Z/4 - I*T**2*\Delta_k*gamma_32*v_Z/4 + I*T**2*c*gamma_32*k_z/8 - T**2*gamma_32**2/16 - k_z*(-T**2*\Delta_k*c*v_Z/4 + I*T**2*c*gamma_32/8 - I*Omega_2a**2*T**2*c/(8*Gamma_21) + I*T**2*c**2*kappa**2/(8*Gamma_21)) + I*Omega_2a**2*T**2*\Delta_k*v_Z/(4*Gamma_21) - I*Omega_2a**2*T**2*c*k_z/(8*Gamma_21) + Omega_2a**2*T**2*gamma_32/(8*Gamma_21) - I*T**2*\Delta_k*c*kappa**2*v_Z/(4*Gamma_21) + I*T**2*c**2*k_z*kappa**2/(8*Gamma_21) - T**2*c*gamma_32*kappa**2/(8*Gamma_21) - Omega_2a**4*T**2/(16*Gamma_21**2) - Omega_2a**2*T**2*c*kappa**2/(8*Gamma_21**2) - T**2*c**2*kappa**4/(16*Gamma_21**2)






    Out[217]:





$$\left [ \frac{\Delta_k^{2} v_{Z}^{2}}{4} T^{2} - \frac{\Delta_k k_{z}}{4} T^{2} c v_{Z} - \frac{i \Delta_k}{4} T^{2} \gamma_{32} v_{Z} + \frac{i \gamma_{32}}{8} T^{2} c k_{z} - \frac{T^{2} \gamma_{32}^{2}}{16} - k_{z} \left(- \frac{\Delta_k v_{Z}}{4} T^{2} c + \frac{i \gamma_{32}}{8} T^{2} c - \frac{i \Omega_{2a}^{2} T^{2} c}{8 \Gamma_{21}} + \frac{i T^{2} c^{2} \kappa^{2}}{8 \Gamma_{21}}\right) + \frac{i \Omega_{2a}^{2} \Delta_k v_{Z}}{4 \Gamma_{21}} T^{2} - \frac{i \Omega_{2a}^{2} c k_{z}}{8 \Gamma_{21}} T^{2} + \frac{\Omega_{2a}^{2} T^{2} \gamma_{32}}{8 \Gamma_{21}} - \frac{i \Delta_k c v_{Z}}{4 \Gamma_{21}} T^{2} \kappa^{2} + \frac{i T^{2} c^{2} k_{z}}{8 \Gamma_{21}} \kappa^{2} - \frac{T^{2} c \gamma_{32} \kappa^{2}}{8 \Gamma_{21}} - \frac{\Omega_{2a}^{4} T^{2}}{16 \Gamma_{21}^{2}} - \frac{\Omega_{2a}^{2} T^{2} c \kappa^{2}}{8 \Gamma_{21}^{2}} - \frac{T^{2} c^{2} \kappa^{4}}{16 \Gamma_{21}^{2}}, \quad - \frac{\Delta_k v_{Z}}{4} T^{2} c + \frac{i \gamma_{32}}{8} T^{2} c - \frac{i \Omega_{2a}^{2} T^{2} c}{8 \Gamma_{21}} + \frac{i T^{2} c^{2} \kappa^{2}}{8 \Gamma_{21}}, \quad \frac{T^{2} c^{2}}{16}\right ]$$



In [216]:

    
from sympy import diff
def coeffs(expr, x, deg=1):
    coefs = []
    rem = expr
    for n in reversed(range(1, deg+1)):
        coef = diff(rem, x, n)/n
        term = coef*x**n
        rem = rem - term
        coefs += [coef]
    print rem
    coefs += [rem]
    return list(reversed(coefs))



In [ ]:



In [198]:

    
expr = 3*kz**2 + 2*kz+1
coeffs(expr, kz, 2)









    



3 3*k_z**2 2*k_z + 1
2 2*k_z 1
[3, 2, 1] 1



In [ ]:



In [173]:

    
aux = F_guillaume.args[0]*F_guillaume.args[1]*F_guillaume.args[2]
aux = aux.factor()
arg =









    Out[173]:





$$\left ( 4, \quad c, \quad \Omega_{2a}^{2}, \quad \kappa^{2}, \quad \frac{1}{4 \Gamma_{21}^{2} \Delta_k^{2} v_{Z}^{2} - 4 \Gamma_{21}^{2} \Delta_k c k_{z} v_{Z} - 4 i \Gamma_{21}^{2} \Delta_k \gamma_{32} v_{Z} + \Gamma_{21}^{2} c^{2} k_{z}^{2} + 2 i \Gamma_{21}^{2} c \gamma_{32} k_{z} - \Gamma_{21}^{2} \gamma_{32}^{2} + 4 i \Gamma_{21} \Omega_{2a}^{2} \Delta_k v_{Z} - 2 i \Gamma_{21} \Omega_{2a}^{2} c k_{z} + 2 \Gamma_{21} \Omega_{2a}^{2} \gamma_{32} - 4 i \Gamma_{21} \Delta_k c \kappa^{2} v_{Z} + 2 i \Gamma_{21} c^{2} k_{z} \kappa^{2} - 2 \Gamma_{21} c \gamma_{32} \kappa^{2} - \Omega_{2a}^{4} - 2 \Omega_{2a}^{2} c \kappa^{2} - c^{2} \kappa^{4}}, \quad \sin^{2}{\left (T \sqrt{\frac{\Delta_k^{2} v_{Z}^{2}}{4} - \frac{\Delta_k k_{z}}{4} c v_{Z} - \frac{i \Delta_k}{4} \gamma_{32} v_{Z} + \frac{c^{2} k_{z}^{2}}{16} + \frac{i \gamma_{32}}{8} c k_{z} - \frac{\gamma_{32}^{2}}{16} + \frac{i \Omega_{2a}^{2} \Delta_k v_{Z}}{4 \Gamma_{21}} - \frac{i \Omega_{2a}^{2} c k_{z}}{8 \Gamma_{21}} + \frac{\Omega_{2a}^{2} \gamma_{32}}{8 \Gamma_{21}} - \frac{i \Delta_k c v_{Z}}{4 \Gamma_{21}} \kappa^{2} + \frac{i c^{2} k_{z} \kappa^{2}}{8 \Gamma_{21}} - \frac{c \gamma_{32} \kappa^{2}}{8 \Gamma_{21}} - \frac{\Omega_{2a}^{4}}{16 \Gamma_{21}^{2}} - \frac{\Omega_{2a}^{2} c \kappa^{2}}{8 \Gamma_{21}^{2}} - \frac{c^{2} \kappa^{4}}{16 \Gamma_{21}^{2}}} \right )}\right )$$



In [162]:

    
ss_M_inv









    Out[162]:





$$\left \{ p : \frac{c k_{z}}{2} + \frac{i c \kappa^{2}}{2 \Gamma_{21}}, \quad q : \frac{i \Omega_{2a} \sqrt{c} \kappa}{2 \Gamma_{21}} e^{- i \phi}, \quad r : \frac{i \Omega_{2a} \sqrt{c} \kappa}{2 \Gamma_{21}} e^{i \phi}, \quad s : \Delta_k v_{Z} - \frac{i \gamma_{32}}{2} + \frac{i \Omega_{2a}^{2}}{2 \Gamma_{21}}\right \}$$



In [ ]:



In [ ]:



In [ ]:

Decoherence

Solving the second equation in the abscence of a control field we obtain a spin wave proportional to



In [151]:

    
Delta, sigma = symbols("Delta sigma", positive=True)
f = exp(I*Delta*t*vZ)
f









    Out[151]:





$$e^{i \Delta t v_{Z}}$$



In [152]:

    
g=1/sqrt(2*pi)/sigma*exp(-(vZ/sigma)**2/2)
g









    Out[152]:





$$\frac{\sqrt{2} e^{- \frac{v_{Z}^{2}}{2 \sigma^{2}}}}{2 \sqrt{\pi} \sigma}$$



In [153]:

    
Doppler_decoherence=integrate(g*f, (vZ, -oo, oo)).factor()
Doppler_decoherence









    Out[153]:





$$e^{- \frac{\Delta^{2} t^{2}}{2} \sigma^{2}}$$



In [154]:

    
decoherence = Doppler_decoherence*exp(-gamma[2, 1]*t/2)
decoherence









    Out[154]:





$$e^{- \frac{\gamma_{32} t}{2}} e^{- \frac{\Delta^{2} t^{2}}{2} \sigma^{2}}$$



In [155]:

    
solve(log(Doppler_decoherence**2)+1, t)









    Out[155]:





$$\left [ - \frac{1}{\Delta \sigma}, \quad \frac{1}{\Delta \sigma}\right ]$$



In [156]:

    
solve(log(decoherence**2)+1, t)









    Out[156]:





$$\left [ \frac{1}{2 \Delta^{2} \sigma^{2}} \left(- \gamma_{32} + \sqrt{4 \Delta^{2} \sigma^{2} + \gamma_{32}^{2}}\right), \quad - \frac{1}{2 \Delta^{2} \sigma^{2}} \left(\gamma_{32} + \sqrt{4 \Delta^{2} \sigma^{2} + \gamma_{32}^{2}}\right)\right ]$$



In [157]:

    
A, B, C = symbols("A, B, C")
omega87, omega97 = symbols("omega87, omega97", positive=True)



In [158]:

    
decoherence_hfs = decoherence*(A+B*exp(I*omega87*t)+C*exp(I*omega97*t))/(A+B+C)
decoherence_hfs









    Out[158]:





$$\frac{e^{- \frac{\gamma_{32} t}{2}} e^{- \frac{\Delta^{2} t^{2}}{2} \sigma^{2}}}{A + B + C} \left(A + B e^{i \omega_{87} t} + C e^{i \omega_{97} t}\right)$$



In [159]:

    
eta_hfs = decoherence_hfs*decoherence_hfs.conjugate()
eta_hfs









    Out[159]:





$$\frac{e^{- \gamma_{32} t} e^{- \Delta^{2} \sigma^{2} t^{2}}}{\left(A + B + C\right) \left(\overline{A} + \overline{B} + \overline{C}\right)} \left(A + B e^{i \omega_{87} t} + C e^{i \omega_{97} t}\right) \left(\overline{A} + e^{- i \omega_{97} t} \overline{C} + e^{- i \omega_{87} t} \overline{B}\right)$$



In [160]:

    
D = symbols("D", positive=True)



In [161]:

    
eq = eta_hfs
eq = eq.subs({(A+B+C)*(A+B+C).conjugate(): 1}) -exp(-1)*D
eq = eq.expand()
eq









    Out[161]:





$$A e^{- \gamma_{32} t} e^{- \Delta^{2} \sigma^{2} t^{2}} \overline{A} + A e^{- \gamma_{32} t} e^{- i \omega_{97} t} e^{- \Delta^{2} \sigma^{2} t^{2}} \overline{C} + A e^{- \gamma_{32} t} e^{- i \omega_{87} t} e^{- \Delta^{2} \sigma^{2} t^{2}} \overline{B} + B e^{- \gamma_{32} t} e^{i \omega_{87} t} e^{- \Delta^{2} \sigma^{2} t^{2}} \overline{A} + B e^{- \gamma_{32} t} e^{i \omega_{87} t} e^{- i \omega_{97} t} e^{- \Delta^{2} \sigma^{2} t^{2}} \overline{C} + B e^{- \gamma_{32} t} e^{- \Delta^{2} \sigma^{2} t^{2}} \overline{B} + C e^{- \gamma_{32} t} e^{i \omega_{97} t} e^{- \Delta^{2} \sigma^{2} t^{2}} \overline{A} + C e^{- \gamma_{32} t} e^{- \Delta^{2} \sigma^{2} t^{2}} \overline{C} + C e^{- \gamma_{32} t} e^{- i \omega_{87} t} e^{i \omega_{97} t} e^{- \Delta^{2} \sigma^{2} t^{2}} \overline{B} - \frac{D}{e}$$



In [ ]:



In [ ]:

[1] H.J. Metcalf and P. van der Straten. Laser Cooling and Trapping. Graduate Texts in Contempo- rary Physics. Springer New York, 2001.