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
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)
Out[5]:
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)
Out[6]:
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)
Out[7]:
In [8]:
fprint(omega_laser,print_ascii=print_ascii)
Out[8]:
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)
Out[9]:
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)
Out[10]:
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)
Out[11]:
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)
Out[13]:
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)
Out[14]:
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)
Out[15]:
In [16]:
fprint(omega, print_ascii=print_ascii)
Out[16]:
In [17]:
fprint(gamma, print_ascii=print_ascii)
Out[17]:
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)
Out[18]:
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)
Out[19]:
and the complete hamiltonian is
In [20]:
H=H0+H1
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)
Out[21]:
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)
Out[22]:
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)
Out[23]:
In [24]:
fprint(r_I[1], print_ascii=print_ascii)
Out[24]:
In [25]:
fprint(r_I[2], print_ascii=print_ascii)
Out[25]:
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)
Out[26]:
In [27]:
fprint(r_I_p[1], print_ascii=print_ascii)
Out[27]:
In [28]:
fprint(r_I_p[2], print_ascii=print_ascii)
Out[28]:
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)
Out[29]:
In [30]:
fprint(r_I_m[1], print_ascii=print_ascii)
Out[30]:
In [31]:
fprint(r_I_m[2], print_ascii=print_ascii)
Out[31]:
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)
Out[32]:
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)
Out[33]:
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)
Out[34]:
The matrix element $(\hat{H}_{1I,RWA})_{21}$ element is
In [35]:
fprint(H1IRWA[1,0].expand(),print_ascii=print_ascii)
Out[35]:
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)
Out[36]:
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)
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)
Out[38]:
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)
Out[39]:
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)
Out[40]:
In [41]:
fprint( [r[p]-(r_p[p]+r_m[p]) for p in range(3)] , print_ascii=print_ascii)
Out[41]:
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)
Out[42]:
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)
Out[43]:
In [44]:
cc,cctilde,phase=define_psi_coefficients(Ne)
fprint([cc,cctilde,phase], print_ascii=print_ascii)
Out[44]:
In [45]:
phase=Matrix([ Function("theta_"+str(i+1),real=True)(t,Z) for i in range(Ne)])
phase
Out[45]:
In [46]:
psi=Matrix([ exp(I*phase[i])*cctilde[i] for i in range(Ne)])
fprint(psi, print_ascii=print_ascii)
Out[46]:
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)
Out[47]:
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)
Out[49]:
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)
Out[50]:
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
Out[51]:
Thus the equations become
In [52]:
fprint(lhs_new, print_ascii=print_ascii)
Out[52]:
In [53]:
rhs_new=simplify(rhs_new.subs(pt)).expand().simplify()
fprint(rhs_new, print_ascii=print_ascii)
Out[53]:
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)
Out[54]:
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)
Out[55]:
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]:
In [57]:
rho=define_density_matrix(Ne, variables=[t, Z])
fprint( rho , print_ascii=print_ascii)
Out[57]:
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)
Out[58]:
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)
Out[59]:
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]:
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]:
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]:
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]:
In [67]:
lhs=lhs.subs(svea_subs)
lhs.simplify()
lhs
Out[67]:
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
Out[68]:
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
Out[69]:
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
Out[71]:
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
Out[74]:
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
Out[75]:
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]:
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]:
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]:
In [82]:
vh=Matrix(symbols("r_-,r0,r_+",real=True))
vh,helicity_to_cartesian(vh)
Out[82]:
In [83]:
rpl1_cartesian=[(rpl1[0]-rpl1[2])/sqrt(2),(rpl1[0]+rpl1[2])*I/sqrt(2),rpl1[1]]
rpl1_cartesian
Out[83]:
In [84]:
rpl2_cartesian=[(rpl2[0]-rpl2[2])/sqrt(2),(rpl2[0]+rpl2[2])*I/sqrt(2),rpl2[1]]
rpl2_cartesian
Out[84]:
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
Out[85]:
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
Out[86]:
In [87]:
rhs1=rhs1.subs({P0[0]:Ppl1})
rhs2=rhs2.subs({P0[1]:Ppl2})
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)
eqs_list=eqs_list.subs({mu0:1/epsilon0/c**2})
fprint(eqs_list, print_ascii=print_ascii)
Out[88]:
In [89]:
vX, vY, vZ=symbols("v_X, v_Y, v_Z", real=True)
v = Matrix([vX, vY, vZ])
v
Out[89]:
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 [90]:
rhop=define_density_matrix(Ne, variables=[t,Z,vZ])
g=Function("g")(vZ)
g
Out[90]:
In [91]:
delta1p=delta1-omega_laser[0]*vZ/c
delta2p=delta2+omega_laser[1]*vZ/c
delta1p,delta2p
Out[91]:
So our final equations are:
In [92]:
rho_subs={rho[i,j]:rhop[i,j] for i in range(Ne) for j in range(Ne)}
rho_subs
eqs_list = eqs_list.subs({delta1: delta1p, delta2: delta2p})
eqs_list = eqs_list.subs(rho_subs)
eqs_list[6,:] = eqs_list[6,:].subs({rhop[1,0]: Integral(g*rhop[1,0],(vZ,-oo,oo))})
eqs_list[7,:] = eqs_list[7,:].subs({rhop[2,1]: Integral(g*rhop[2,1],(vZ,-oo,oo))})
eqs_list
Out[92]:
In [93]:
m, kB, T = symbols("m, k_B, T", positive=True)
vZ = symbols("v_Z", real=True)
m, kB, T
Out[93]:
In [94]:
sigma_vZ=sqrt(kB*T/m)
g=sqrt(1/(2*pi))/sigma_vZ*exp(-(vZ/sigma_vZ)**2/2)
g
Out[94]:
In [95]:
from sympy import integrate, oo
integrate(g,(vZ,-oo,oo))
Out[95]:
In [ ]:
[1] H.J. Metcalf and P. van der Straten. Laser Cooling and Trapping. Graduate Texts in Contempo- rary Physics. Springer New York, 2001.
[2] Daniel Adam Steck. Quantum and Atom Optics. Oregon Center for Optics and Department of Physics, University of Oregon Copyright © 200
In [ ]: