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.
Let us check how many hyperfine states each atom has:
In [3]:
atoms = [["Rb", 85, 5], ["Rb", 87, 5], ["Cs", 133, 6]]
for ii in atoms:
atom = Atom(ii[0], ii[1])
e1 = State(ii[0], ii[1], ii[2], 0, 1/Integer(2))
e2 = State(ii[0], ii[1], ii[2], 1, 3/Integer(2))
e3 = State(ii[0], ii[1], ii[2], 2, 5/Integer(2))
fine_states = [e1, e2, e3]
Nhyp = len(make_list_of_states(fine_states, "hyperfine", verbose=0))
Nmag = len(make_list_of_states(fine_states, "magnetic", verbose=0))
print atom, "has", Nhyp, "hyperfine states, and",Nmag, "magnetic states."
ii = atoms[2]
e1 = State(ii[0], ii[1], ii[2], 0, 1/Integer(2))
e2 = State(ii[0], ii[1], ii[2], 1, 3/Integer(2))
e3 = State(ii[0], ii[1], ii[2], 2, 5/Integer(2))
fine_states = [e1, e2, e3]
hyperfine_states = make_list_of_states(fine_states, "hyperfine", verbose=0)
magnetic_states = make_list_of_states(fine_states, "magnetic", verbose=0)
Nhyp = len(hyperfine_states)
Nmag = len(magnetic_states)
Ne = Nhyp
In [4]:
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 [5]:
Nl=2
We define a few important symbols.
In [6]:
c,hbar,e,mu0, epsilon0=symbols("c hbar e mu0 varepsilon0",positive=True)
# fprint([c,hbar,e,mu0, epsilon0],print_ascii=print_ascii)
In [7]:
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 [8]:
E0,omega_laser=define_laser_variables(Nl, variables=[t,Z])
# fprint(E0,print_ascii=print_ascii)
In [9]:
# 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 [10]:
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 [11]:
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 [12]:
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 [13]:
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 [14]:
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 [15]:
r=define_r_components(Ne,helicity=True,explicitly_hermitian=True)
for i, ei in enumerate(hyperfine_states):
for j, ej in enumerate(hyperfine_states):
if not Transition(ei, ej).allowed:
r=[ri.subs({r[0][i,j]:0,r[1][i,j]:0,r[2][i,j]:0}) for ri in r]
The frequencies of the energy levels, the resonant frequencies, and the decay frequencies.
In [16]:
omega_level,omega,gamma=define_frequencies(Ne,explicitly_antisymmetric=True)
for i, ei in enumerate(hyperfine_states):
for j, ej in enumerate(hyperfine_states):
if not Transition(ei, ej).allowed:
gamma=gamma.subs({gamma[i,j]: 0})
# fprint(omega_level,print_ascii=print_ascii)
The atomic hamiltonian is
In [17]:
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 [18]:
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 [19]:
H=H0+H1
In [20]:
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 [21]:
# 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 [22]:
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 [23]:
# fprint(r_I[1], print_ascii=print_ascii)
In [24]:
# fprint(r_I[2], print_ascii=print_ascii)
Which can be decomposed in positive and negative frequencies as
In [25]:
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 [26]:
# fprint(r_I_p[1], print_ascii=print_ascii)
In [27]:
# fprint(r_I_p[2], print_ascii=print_ascii)
In [28]:
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 [29]:
# fprint(r_I_m[1], print_ascii=print_ascii)
In [30]:
# fprint(r_I_m[2], print_ascii=print_ascii)
that summed equal $\vec{\hat{r}}_I$
In [31]:
# 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 [32]:
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 [33]:
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 [34]:
# 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 [35]:
# 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 [36]:
Lij=[[1,2,[1]],[2,3,[2]]]
Lij1 = [[[i+1, j+1, [1]] for j in range(i)
if Transition(hyperfine_states[i], hyperfine_states[j]).allowed] for i in range(6)]
Lij2 = [[[i+1, j+1, [2]] for j in range(i)
if Transition(hyperfine_states[i], hyperfine_states[j]).allowed] for i in range(6, Ne)]
Lij1 = sum(Lij1, [])
Lij2 = sum(Lij2, [])
Lij = Lij1 + Lij2
# print Lij
Lij=formatLij(Lij,Ne)
print array(Lij)
Thus the interacion hamiltonian in the interaction picture can be approximated as
In [37]:
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 [38]:
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 [39]:
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 [40]:
# 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 [41]:
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 [42]:
HRWA=H0+H1RWA
# fprint(HRWA, print_ascii=print_ascii)
In [43]:
cc,cctilde,phase=define_psi_coefficients(Ne)
# fprint([cc,cctilde,phase], print_ascii=print_ascii)
In [44]:
phase=Matrix([ Function("theta_"+str(i+1),real=True)(t,Z) for i in range(Ne)])
# phase
In [45]:
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 [46]:
lhs=Matrix([(I*hbar*Derivative(psi[i],t).doit()).expand() for i in range(Ne)])
# fprint(lhs, print_ascii=print_ascii)
In [47]:
rhs=HRWA*psi
We multiply each of these equations by $e^{-i \theta_i}$ and substracting $i \theta_i \tilde{c}_i$
In [48]:
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 [49]:
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 [50]:
from fast import bloch
import numpy as np
In [51]:
xi = np.zeros((Nl, Ne, Ne))
for l in range(Nl):
for i in range(Ne):
for j in range(Ne):
if l+1 in Lij[i][j]:
xi[l, i, j] = 1
In [52]:
pt = bloch.phase_transformation(Ne, Nl, r_m, xi)
#pt = {phase[i]: pt[i] for i in range(Ne)}
ss ={t*omega_laser[0]: (t-Z/c)*omega_laser[0],
t*omega_laser[1]: (t+Z/c)*omega_laser[1]}
pt = {phase[i]: (pt[i]*t).expand().subs(ss).expand() for i in range(Ne)}
Thus the equations become
In [53]:
rhs_new2=simplify(rhs_new.subs(pt)).doit().expand().simplify()
# fprint(rhs_new2, 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_new2[i],cctilde[j]).doit().simplify() for j in range(Ne)] for i in range(Ne)])
# fprint(Htilde, print_ascii=print_ascii)
In [55]:
Om131, Om141, Om142, Om151, Om152, Om162 = symbols("Omega_131, Omega_141, Omega_142, Omega_151, Omega_152, Omega_162")
Om276, Om286, Om296, Om285, Om295, Om2105 = symbols("Omega_276, Omega_286, Omega_296, Omega_285, Omega_295, Omega_2105")
Om294, Om2104, Om2114, Om2103, Om2113, Om2123 = symbols("Omega_294, Omega_2104, Omega_2114, Omega_2103, Omega_2113, Omega_2123")
delta1, delta2 = symbols("delta1, delta2", real=True)
ssOm ={r[2][2,0]*E0[0]: Om131/e*hbar,
r[2][3,0]*E0[0]: Om141/e*hbar,
r[2][3,1]*E0[0]: Om142/e*hbar,
r[2][4,0]*E0[0]: Om141/e*hbar,
r[2][4,1]*E0[0]: Om142/e*hbar,
r[2][5,1]*E0[0]: Om152/e*hbar,
r[0][6, 5]*E0[1]: Om276 /e*hbar,
r[0][7, 5]*E0[1]: Om286 /e*hbar,
r[0][8, 5]*E0[1]: Om296 /e*hbar,
r[0][7, 4]*E0[1]: Om285 /e*hbar,
r[0][8, 4]*E0[1]: Om295 /e*hbar,
r[0][9, 4]*E0[1]: Om2105/e*hbar,
r[0][8, 3]*E0[1]: Om294 /e*hbar,
r[0][9, 3]*E0[1]: Om2104/e*hbar,
r[0][10,3]*E0[1]: Om2114/e*hbar,
r[0][9, 2]*E0[1]: Om2103/e*hbar,
r[0][10,2]*E0[1]: Om2113/e*hbar,
r[0][11,2]*E0[1]: Om2123/e*hbar}
ss_delta = {omega_laser[0]: delta1 +(omega_level[2]-omega_level[1]),
omega_laser[1]: delta2 +(omega_level[6]-omega_level[5])}
ssOm2 = {i.conjugate(): ssOm[i].conjugate() for i in ssOm}
ssOm.update(ssOm2)
In [56]:
from sympy import eye
Htilde = Htilde - eye(Ne)*hbar*(omega_level[1]-omega_level[0])
Htilde = Htilde.expand()
# Htilde = Htilde.subs(ssOm).expand()
Htilde = Htilde.subs(ss_delta).expand()
Htilde
Out[56]:
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 [57]:
vX, vY, vZ=symbols("v_X, v_Y, v_Z", real=True)
rho=define_density_matrix(Ne, variables=[t, Z, vZ])
# fprint( rho , print_ascii=print_ascii)
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]:
lindblads = lindblad_terms(gamma, rho, Ne)
The Optical Bloch equations are thus.
In [60]:
eqs=hamiltonian_terms + lindblads
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)
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
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
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
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,Z), Function("P_0^2")(t,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)
# rhs
In [73]:
eqs_wave=lhs-rhs
# 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
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]:
Om1 = Function("Omega1")(t, Z)
Om2 = Function("Omega2")(t, Z)
{E0[0]:Om131*hbar/r[2][2,0]/e}
fact1=e*r[2][2,0]/hbar
# fact1, lhs1, (fact1*lhs1.subs({E0[0]:Om1*hbar/r[2][2,0]/e}).doit()).expand()
In [78]:
fact1=e*r[2][2,0]/hbar
# lhs1=(fact1*lhs1.subs({E0[0]:Om1*hbar/r[2][2,0]/e}).doit()).expand()
# rhs1=(fact1*rhs1.subs({E0[0]:Om1*hbar/r[2][2,0]/e}).doit()).expand()
# Matrix([lhs1,eqsign,rhs1]).transpose()
In [79]:
fact2=e*r[0][6,5]/hbar
# lhs2=(fact2*lhs2.subs({E0[1]:Om2*hbar/r[0][6,5]/e}).doit()).expand()
# rhs2=(fact2*rhs2.subs({E0[1]:Om2*hbar/r[0][6,5]/e}).doit()).expand()
# Matrix([lhs2,eqsign,rhs2]).transpose()
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 [80]:
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 [81]:
rpl1=r_p_component[0]
rpl2=r_p_component[1]
In [82]:
n=Function("n")(R,Z)
n
Out[82]:
In [83]:
vh=Matrix(symbols("r_-,r0,r_+",real=True))
# vh,helicity_to_cartesian(vh)
In [84]:
rpl1_cartesian=[(rpl1[0]-rpl1[2])/sqrt(2),(rpl1[0]+rpl1[2])*I/sqrt(2),rpl1[1]]
# rpl1_cartesian
In [85]:
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 [86]:
Ppl1_v=-2*n*e*Matrix([ (rpl1_cartesian[i]*rho).trace() for i in range(3)])
Ppl2_v=-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 [87]:
Ppl1=ep[0].conjugate().dot(Ppl1_v).expand()
Ppl2=ep[1].conjugate().dot(Ppl2_v).expand()
Ppl1.factor()
Out[87]:
In [ ]:
In [88]:
rhs1=rhs1.subs({P0[0]:Ppl1})
rhs2=rhs2.subs({P0[1]:Ppl2})
In [89]:
wave_eq1 = rhs1
rhs1
Out[89]:
In [90]:
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)
Question: how to use calculations from detunings as a replacement for calculations for velocity classes.
In [91]:
epsilon=symbols("varepsilon", real=True)
lin_subs1={rho[i, i]: 0 for i in range(Ne) if hyperfine_states[i].l != 0}
lin_subs1.update({rho[0, 0]: 0, rho[1, 1]: 1})
In [92]:
lin_subs2 = {}
for i in range(Ne):
for j in range(i):
lin_subs2.update({rho[i, j]: epsilon*rho[i, j]})
lin_subs2.update({rho[j, i]: epsilon*rho[j, i]})
lin_subs2.update({E0[0]: epsilon*E0[0]})
# lin_subs2
In [93]:
lin_subs = {}
lin_subs.update(lin_subs1)
lin_subs.update(lin_subs2)
In [94]:
hyp_bound, mag_bound = calculate_boundaries(fine_states, magnetic_states)
In [95]:
hyp_bound
Out[95]:
In [96]:
rho21 = rho[2: 6, 0:2]
rho31 = rho[6:Ne , 0: 2]
rho32 = rho[6:Ne ,2: 6]
# rho21
In [97]:
lims21 = [[2, 6], [0,2]]
lims31 = [[6, Ne], [0,2]]
lims32 = [[6, Ne], [2,6]]
rho21 = rho[lims21[0][0]: lims21[0][1], lims21[1][0]: lims21[1][1]]
rho31 = rho[lims31[0][0]: lims31[0][1], lims31[1][0]: lims31[1][1]]
rho32 = rho[lims32[0][0]: lims32[0][1], lims32[1][0]: lims32[1][1]]
# rho32
In [98]:
# We make a vector of components of the density matrix rho21, rho31, and E1
rhov = sum([[rho[j, i] for j in range(i+1, Ne) if ((rho[j, i] in rho21)or(rho[j, i] in rho31))] for i in range(Ne)], [])
xx = rhov + [E0[0]]
Nxx = len(xx)
In [99]:
# We make a similar vector with the non linear right hand sides.
eqs_mem = sum([[eqs[j, i] for j in range(i+1, Ne) if ((rho[j, i] in rho21)or(rho[j, i] in rho31))] for i in range(Ne)], [])
eqs_mem += [wave_eq1.subs({E0[0]:0}).doit()]
In [100]:
# We get the linear right hand sides.
A = []
eqs_aprox = []
for mu in range(Nxx):
eq_mu = eqs_mem[mu]
eq_approx = eq_mu.subs(lin_subs).subs({epsilon**2: 0}).subs({epsilon: 1})
eqs_aprox += [eq_approx]
eq_lin_mu = [Derivative(eq_approx, xx[nu]).doit() for nu in range(Nxx)]
A+=[eq_lin_mu]
eqs_aprox = Matrix(eqs_aprox)
A = Matrix(A)
A
Out[100]:
In [101]:
# We prove that the linear aproximation produced linear equations.
zero = eqs_aprox - A*Matrix(xx)
zero = zero.expand()
pprint(zero.T)
We rewrite the equations in a factorized form
In [102]:
from sympy import diff, I
def semi_factor(expr0, f1, f2):
expr = expr0
aux = symbols("aux")
expr = expr.subs({f1: aux})
d1 = diff(expr, aux).factor()
r1 = (expr- d1*aux).expand()
expr = expr.subs({aux: f1})
expr = expr.subs({f2: aux})
d2 = diff(expr, aux)
r2 = (expr-d2*aux).expand()
return (d1*f1+d2*f2)
In [103]:
eqs_aprox2 =[semi_factor(eqs_aprox[i], I*e/hbar/2, xx[i]) for i in range(Nxx)]
eqs_aprox2[-1] = eqs_aprox[-1].factor()
zero = [(eqs_aprox[i]-eqs_aprox2[i]).expand() for i in range(Nxx)]
zero
Out[103]:
We get code to calculate the right hand sides:
In [104]:
fact = symbols("fact")
def format_expr(expr):
s = str(expr)
s = s.replace("\\", "")
s = s.replace("{", "")
s = s.replace("}", "")
s = s.replace(",", "_")
s = s.replace("_+1;", "_1_")
s = s.replace("_-1;", "_2_")
s = s.replace("E_0^2", "E02")
s = s.replace("E_0^1", "E01")
s = s.replace("I", "1j")
s = s.replace("/2", "*0.5")
s = s.replace("n(R_ Z)", "n_atomic")
s = s.replace("(t_ Z_ v_Z)", "")
s = s.replace("(t_ Z)", "")
s = s.replace("conjugate(E02)", "E02c")
return s
for i, eq in enumerate(eqs_aprox2[:-1]):
eqi = eq.subs({e: 2*hbar*fact/I})
codei = "eqs_rho["+str(i)+", jj] = "
codei += format_expr(eqi)
print codei
print
print "eqs_rho[20, jj] =", format_expr(eqs_aprox[-1].factor())
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 [105]:
rhop=define_density_matrix(Ne, variables=[t,Z,vZ])
g=Function("g")(vZ)
g
Out[105]:
In [106]:
k1, k2 = symbols("k1 k2", positive=True)
delta1_Doppler=delta1 + k1*vZ
delta2_Doppler=delta2 - k2*vZ
delta1_Doppler,delta2_Doppler
ss_Doppler = {delta1: delta1_Doppler, delta2: delta2_Doppler}
ss_Doppler
Out[106]:
The equations involving states $|1\rangle$ and $|2\rangle$ are decoupled. If we choose the higher hyperfine state to be our ORCA ground state, we finally get equations:
In [107]:
ss_higher = {xx[0]: 0, xx[1]: 0, xx[2]: 0}
ss_higher
Out[107]:
In [108]:
ss_wave = {xx[11]: Integral(g*xx[11], vZ), xx[12]: Integral(g*xx[12], vZ), xx[13]: Integral(g*xx[13], vZ)}
ss_wave
Out[108]:
In [109]:
eq_sign = symbols("=")
eqs_aprox3 = [[diff(xx[i], t), eq_sign, eqs_aprox2[i].subs(ss_higher).subs(ss_Doppler)] for i in range(Nxx)]
eqs_aprox3 = Matrix(eqs_aprox3)[10:, :]
eqs_aprox3[-1, 2] = -c*diff(xx[-1], Z)+eqs_aprox3[-1, 2].subs(ss_wave)
eqs_aprox3
Out[109]:
[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