text:
# Approximate Master Equations for Atom Optics

## I The Master Equation

Many investigations in the realm of quantum physics are based on a single two level atom coupled to an approximately resonant light field. One major aspect of this work is to investigate the quantum mechanical motion of such an atom. The work of Steck et al Steck  is a prime example. They investigate chaos assisted tunneling by studying the motion of cold cesium atoms in an amplitude modulated standing wave of light. Haycock et al Haycock  study cesium atoms in an effort to observe quantum coherent dynamics. The work of Hensinger et al investigates quantum chaos Hen01a  and quantum tunneling Hen01b . Proposals in QED to utilize interactions between atoms and their cavity for physical realizations of a quantum computer Pellizzari  and for communication of quantum states Cirac  also require good descriptions of quantum mechanical motion.

The full quantum system here consists of the light field and the atom. The evolution of such a system is governed by a Schrödinger equation. In most cases however, the evolution of the light field is not of interest and so the Schrödinger equation can be reduced to a master equation. The most general form of this master equation is given by

$$\dot{\rho }=\frac{i}{\mathrm{}}[H,\rho ]+\rho ,$$
(1)
where $``$ is a Lindbladian superoperator Lind .

The form we will use here, specific to a two level atom interacting with a light field, is

$$H=\frac{p_x^2}{2m}+\mathrm{}\delta \sigma ^{}\sigma +\frac{\mathrm{}}{2}\left[\mathrm{\Omega }(x,t)\sigma ^{}+\mathrm{\Omega }^{}(x,t)\sigma \right],$$
(2a)
$$\rho =\mathrm{\Gamma }\left(𝒥\left[\sigma \right]\rho 𝒜\left[\sigma \right]\rho \right).$$
(2b)
Here we are only interested in the motion of the atom in the one direction. Thus $`p_x`$ is the component of the atom’s momentum in the $`x`$-direction. The detuning of the light field is $`\delta `$. The atomic lowering operator is given by $`\sigma =|ge|`$ where $`|g`$ and $`|e`$ represent the ground and excited states of the atom respectively. $`\mathrm{\Omega }(x,t)`$ is the complex, and possibly time-dependent Rabi frequency operator for the light field, which, from here on, will be represented simply as $`\mathrm{\Omega }`$. $`k`$ and $`m`$ are the wavenumber of the incident photons and the mass of the atom respectively. Planck’s constant ($`\mathrm{}`$) will be set to one for the rest of the analysis. The superoperator $``$ describes random momentum kicks. It is defined for any arbitrary operator $`a`$ as

$$a=_1^1W(u)e^{ikux}ae^{ikux}𝑑u,$$
(3)
where $`W(u)`$ is the atomic dipole radiation distribution function produced by the electronic transition, reduced to one dimension Hen01a . For motion parallel to the direction of propagation of the laser light, it is given by

$$W(u)=\frac{3}{8}\left(1+u^2\right).$$
(4)
The superoperators $`𝒥`$ and $`𝒜`$ are defined for arbitrary operators $`c`$ and $`a`$ as

$$𝒥\left[c\right]a=cac^{},𝒜\left[c\right]a=\frac{1}{2}\{c^{}c,a\},$$
(5)
and define the general form of a Lindbladian Lind  superoperator.

The entire expression of Eq. (2b) describes the irreversible evolution of the system at rate $`\mathrm{\Gamma }`$. It is this part of the master equation that we wish to minimize by working in the regime $`\delta \mathrm{\Omega }_{\mathrm{max}}\mathrm{\Gamma }`$ (where $`\mathrm{\Omega }_{\mathrm{max}}`$ is the maximum modulus of the Rabi frequency operator). The problem with this regime is that when making $`\delta `$ very large, the full master equation still has to be solved using a timestep smaller than $`\delta ^1`$. This makes solving the full master equation numerically very difficult.

There are a number of ways to approximate the master equation, four of which are investigated here. These are

1. The standard adiabatic approximation;
2. A more sophisticated adiabatic approximation;
3. A secular approximation; and
4. A dressed-state approximation.

Approximation (1) is a standard approach used by many researchers in the field, both as a semi-classical treatment Parkins  and as a fully quantum approximation method Hen01a ; Graham . Unfortunately, this treatment is not valid in the regime of the work of Hen01a . These experiments were performed in the regime of $`\delta \mathrm{\Omega }_{\mathrm{max}}\mathrm{\Gamma }`$ but where $`\mathrm{\Gamma }`$ is the same order as $`\mathrm{\Omega }_{\mathrm{max}}^2/\delta `$. In their work, approximation (1) was used but on closer examination, this approximation was seen only to be valid in the regime $`\mathrm{\Gamma }\mathrm{\Omega }_{\mathrm{max}}^2/\delta `$. We will concentrate on the difficult $`(\mathrm{\Gamma }\mathrm{\Omega }_{\mathrm{max}}^2/\delta )`$ regime in this analysis. Approximation (2) is one way to correct the standard adiabatic approximation for this regime. Approximation (3) was proposed in Ref. dyrtingmilburn , and was used for the formation of quantum trajectory simulations in Ref.Hen01a . Approximation (4) is a fully quantum dressed-state treatment including the effect of spontaneous emission. Previously, only semi-classical treatments, both omitting Deutschmann , and including Dalibard ; Parkins , spontaneous emission have been done. We look at all four approximations, examining the validity and complexity, as well as the numerical accuracy (compared to the full simulation) and computational resource requirements, of each.

## II The standard adiabatic approximation

The method of adiabatic approximation applied to the full master equation serves to eliminate the internal state structure of the atom. One reason for wishing to have no internal states is so that the system has a clear classical analogue dyrtingmilburn . In performing this treatment, we also remove the Hamiltonian term of order $`\delta `$, removing the requirement that the master equation be solved on a timestep at least as small as $`\delta ^1`$. This adiabatic elimination technique is described in many different text books, see Meystre ; Barnett . The result we obtain was first derived by Graham, Schlautmann and Zoller Graham . To achieve this adiabatic elimination however, we follow a procedure similar to that of Hensinger et al Hen01a .

The density matrix can be written using the internal state basis as

$`\rho `$ $`=`$ $`\rho _{gg}|gg|+\rho _{ge}|ge|`$ (6)
$`\rho _{ge}^{}|eg|+\rho _{ee}|ee|,`$
where the $`\rho _{gg}`$ etc are still operators on the centre-of-mass Hilbert space L$`{}_{}{}^{(2)}()`$. From Eqs. (2), these obey

$`\dot{\rho }_{gg}`$ $`=`$ $`\mathrm{\Gamma }\rho _{ee}{\displaystyle \frac{i}{2}}\left(\mathrm{\Omega }^{}\rho _{eg}\rho _{ge}\mathrm{\Omega }\right)𝒦\rho _{gg},`$ (7a)
$`\dot{\rho }_{ge}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}\rho _{ge}{\displaystyle \frac{i}{2}}\left(\mathrm{\Omega }^{}\rho _{ee}\rho _{gg}\mathrm{\Omega }^{}\right)+i\delta \rho _{ge}`$ (7b)
$`𝒦\rho _{ge},`$
$`\dot{\rho }_{ee}`$ $`=`$ $`\mathrm{\Gamma }\rho _{ee}{\displaystyle \frac{i}{2}}\left(\mathrm{\Omega }\rho _{ge}\rho _{eg}\mathrm{\Omega }^{}\right)𝒦\rho _{ee}.`$ (7c)
Here the kinetic energy term is represented in the superoperator

$$𝒦\rho =i[\frac{p_x^2}{2m},\rho ].$$
(8)
If the standard adiabatic elimination procedure is valid, we can represent the system by the density matrix for the centre of mass (com) alone,

$$\rho _{\mathrm{com}}=\mathrm{Tr}_{\mathrm{int}}[\rho ]=\rho _{gg}+\rho _{ee},$$
(9)
where $`\mathrm{Tr}_{\mathrm{int}}`$ is the trace over the internal states of the atom. In this standard adiabatic approximation we further simplify this by noting that large detuning leads to very small excited state populations such that $`\rho _{gg}\rho _{ee}`$ and thus $`\rho _{\mathrm{com}}`$ is approximately $`\rho _{gg}`$. Thus denoting $`\rho _{\mathrm{com}}`$ simply as $`\rho `$, we can replace $`\rho _{gg}`$ in Eq. (7a) just with $`\rho `$.

Now we require expressions for $`\rho _{ge}`$ and $`\rho _{ee}`$ in terms of $`\rho `$. This is achieved by noting that from Eq. (7b), $`\rho _{ge}`$ comes to equilibrium on a timescale much shorter than $`\rho _{gg}`$, at a rate $`\mathrm{\Gamma }/2`$. Thus we set $`\dot{\rho }_{ge}=0`$. Also, if the kinetic energy term is much smaller than $`\delta `$ then it can be ignored. This will be the case if $`\delta p_x^2/m`$. This is typically true, and so we get

$$\rho _{ge}=i\frac{\rho _{gg}\mathrm{\Omega }^{}\mathrm{\Omega }^{}\rho _{ee}}{\mathrm{\Gamma }2i\delta }.$$
(10)
This can be substituted back into the equation for $`\rho _{ee}`$ to give

$`\dot{\rho }_{ee}`$ $`=`$ $`\mathrm{\Gamma }\left(\rho _{ee}+{\displaystyle \frac{\{\mathrm{\Omega }\mathrm{\Omega }^{},\rho _{ee}\}}{2\left(\mathrm{\Gamma }^2+4\delta ^2\right)}}+{\displaystyle \frac{i\delta [\mathrm{\Omega }\mathrm{\Omega }^{},\rho _{ee}]}{\mathrm{\Gamma }\left(\mathrm{\Gamma }^2+4\delta ^2\right)}}\right)`$ (11)
$`+{\displaystyle \frac{\mathrm{\Gamma }𝒥\left[\mathrm{\Omega }\right]\rho _{gg}}{\mathrm{\Gamma }^2+4\delta ^2}}𝒦\rho _{ee}.`$
From Eq. (11), we see that $`\rho _{ee}`$ equilibrates on a timescale much faster than $`\rho _{gg}`$ and so we also set $`\dot{\rho }_{ee}=0`$. This time, the kinetic energy term must be ignored compared with $`\mathrm{\Gamma }`$ rather than $`\delta `$. Allowing this approximation gives

$`\rho _{ee}+{\displaystyle \frac{\{\mathrm{\Omega }\mathrm{\Omega }^{},\rho _{ee}\}}{2\left(\mathrm{\Gamma }^2+4\delta ^2\right)}}+{\displaystyle \frac{i\delta [\mathrm{\Omega }\mathrm{\Omega }^{},\rho _{ee}]}{\mathrm{\Gamma }\left(\mathrm{\Gamma }^2+4\delta ^2\right)}}{\displaystyle \frac{\mathrm{\Omega }\rho _{gg}\mathrm{\Omega }^{}}{\mathrm{\Gamma }^2+4\delta ^2}}.`$ (12)
The first correction term on the left hand side scales as $`(\mathrm{\Omega }/\delta )^2`$ which in the regime chosen is negligible compared to the leading order term. The second correction term however scales as $`\left(\mathrm{\Omega }^2/\delta \right)/\mathrm{\Gamma }`$. Had we been working in the regime $`\mathrm{\Gamma }\mathrm{\Omega }^2/\delta `$, then this term could also be safely ignored compared to the leading term. This condition however, is not satisfied in the experiments of Hen01a  nor in our chosen regime, leaving this term the same order as the leading term. This term however was dropped in Ref. Hen01a  on the basis that in a more sophisticated approach (Sec. III), this term does not appear and the correction to the final master equation is small Hen01a . Knowing that this is an unjustified approximation, but in the interest of comparison to currently used techniques Hen01a , we will continue to follow this method as others have done. Thus, dropping the second correction term we are left with

$$\rho _{ee}\frac{\mathrm{\Omega }\rho _{gg}\mathrm{\Omega }^{}}{\mathrm{\Gamma }^2+4\delta ^2}.$$
(13)
Now substituting first Eq. (10) then Eq. (13) into Eq. (7a) and replacing $`\rho _{gg}`$ with $`\rho `$ gives the final adiabatically eliminated master equation:

$`\dot{\rho }`$ $`=`$ $`\mathrm{\Gamma }\left(𝒥\left[{\displaystyle \frac{\mathrm{\Omega }}{2\delta }}\right]\rho 𝒜\left[{\displaystyle \frac{\mathrm{\Omega }}{2\delta }}\right]\rho \right)i[{\displaystyle \frac{p_x^2}{2m}}{\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{}}{4\delta }},\rho ],`$
where here we have used the fact that $`\delta \mathrm{\Omega }\mathrm{\Gamma }`$ to eliminate some of the higher order terms. This master equation is of the Lindblad form.

## III A more sophisticated Adiabatic approximation

As noted, the standard adiabatic elimination method described in section II is not strictly valid in the regime of the experiments of Hensinger et al Hen01a , $`\mathrm{\Gamma }\mathrm{\Omega }_{\mathrm{max}}^2/\delta `$. There are a number of ways to try and develop a strictly valid version of the adiabatic approximation in this regime. One way would be to not drop any terms without justification and continue to plough through the mathematics. Another method which we believe to be neater and just as accurate is proposed here by using a slightly more sophisticated method similar to that in the appendix of howpolly .

The basis of the approach is to move into an interaction picture with respect to

$$H_0=\frac{\mathrm{\Omega }\mathrm{\Omega }^{}}{4\delta }.$$
(15)
This approach may seem counter-intuitive to most. Usually when moving to an interaction picture, it would be with respect to a $`H_0`$ that is already one of the terms in the Hamiltonian. In our case, if we investigate Eq. (LABEL:adiab1me), we are moving to an interaction picture with respect to the opposite of a term in the effective Hamiltonian and as such will actually be adding a term to the Hamiltonian. The reason we choose to do this is that the problem term we encountered in Sec. II was in the Hamiltonian for the excited state. The potential seen by the excited state of the atom is in fact inverted and so the chosen $`H_0`$ is designed to cancel the excited state potential.

After moving into this interaction picture, we then perform the adiabatic elimination process, and finally transform back into the Schrödinger picture. This method can give a different result because the approximations we make in the interaction picture may not have been valid in the Schrödinger picture.

With the unitary transformation operator $`U_0(t)=e^{iH_0t}`$, the interaction picture density matrix is

$$\stackrel{~}{\rho }=U_0^{}(t)\rho U_0(t).$$
(16)
Using this in Eqs. (2) gives an interaction picture master equation equation still of the Lindblad form

$$\dot{\stackrel{~}{\rho }}=i[\stackrel{~}{V},\stackrel{~}{\rho }]+\mathrm{\Gamma }\stackrel{~}{\rho },$$
(17)
but now with an extra Hamiltonian term such that $`\stackrel{~}{V}`$ is given by

$`\stackrel{~}{V}={\displaystyle \frac{\stackrel{~}{p}_x^2}{2m}}`$ $`+`$ $`\delta \sigma ^{}\sigma +{\displaystyle \frac{1}{2}}\left(\mathrm{\Omega }\sigma ^{}+\mathrm{\Omega }^{}\sigma \right)+{\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{}}{4\delta }}.`$ (18)
Where $`\stackrel{~}{p}_x`$ is just the component of the momentum in the $`x`$-direction, transformed into the interaction picture. The Lindbladian superoperator $``$ is unaffected by the interaction picture because Rabi frequency operator $`\mathrm{\Omega }`$ commutes with the position operator, $`x`$, as well as with the state operators, $`\sigma `$ and $`\sigma ^{}`$.

Following the same procedure as the standard adiabatic elimination process, the equations for the centre-of-mass operators can be extracted:

$`\dot{\stackrel{~}{\rho }}_{gg}`$ $`=`$ $`\mathrm{\Gamma }\stackrel{~}{\rho }_{ee}{\displaystyle \frac{i}{2}}\left(\mathrm{\Omega }^{}\stackrel{~}{\rho }_{eg}\stackrel{~}{\rho }_{ge}\mathrm{\Omega }\right){\displaystyle \frac{i}{4\delta }}[\mathrm{\Omega }\mathrm{\Omega }^{},\stackrel{~}{\rho }_{gg}]`$ (19a)
$`\stackrel{~}{𝒦}\stackrel{~}{\rho }_{gg},`$
$`\dot{\stackrel{~}{\rho }}_{ge}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}\stackrel{~}{\rho }_{ge}{\displaystyle \frac{i}{2}}\left(\mathrm{\Omega }^{}\stackrel{~}{\rho }_{ee}\stackrel{~}{\rho }_{gg}\mathrm{\Omega }^{}\right){\displaystyle \frac{i}{4\delta }}[\mathrm{\Omega }\mathrm{\Omega }^{},\stackrel{~}{\rho }_{ge}]`$ (19b)
$`+i\delta \stackrel{~}{\rho }_{ge}\stackrel{~}{𝒦}\stackrel{~}{\rho }_{ge},`$
$`\dot{\stackrel{~}{\rho }}_{ee}`$ $`=`$ $`\mathrm{\Gamma }\stackrel{~}{\rho }_{ee}{\displaystyle \frac{i}{2}}\left(\mathrm{\Omega }\stackrel{~}{\rho }_{ge}\stackrel{~}{\rho }_{eg}\mathrm{\Omega }^{}\right){\displaystyle \frac{i}{4\delta }}[\mathrm{\Omega }\mathrm{\Omega }^{},\stackrel{~}{\rho }_{ee}]`$ (19c)
$`\stackrel{~}{𝒦}\stackrel{~}{\rho }_{ee},`$
where obviously $`\stackrel{~}{𝒦}`$ uses the interaction picture momentum operator $`\stackrel{~}{p}_x`$.

In this more sophisticated approach, we still take the trace over the internal states of the atom, letting $`\stackrel{~}{\rho }_{\mathrm{com}}=\stackrel{~}{\rho }_{gg}+\stackrel{~}{\rho }_{ee}`$, but now we do not simplify this further and simply let $`\stackrel{~}{\rho }=\stackrel{~}{\rho }_{gg}+\stackrel{~}{\rho }_{ee}`$. This gives a master equation of the form

$`\dot{\stackrel{~}{\rho }}`$ $`=`$ $`\mathrm{\Gamma }\left(1\right)\stackrel{~}{\rho }_{ee}{\displaystyle \frac{i}{2}}[\mathrm{\Omega }^{},\stackrel{~}{\rho }_{eg}]{\displaystyle \frac{i}{2}}[\mathrm{\Omega },\stackrel{~}{\rho }_{ge}]`$ (20)
$`{\displaystyle \frac{i}{4\delta }}[\mathrm{\Omega }\mathrm{\Omega }^{},\stackrel{~}{\rho }]\stackrel{~}{𝒦}\stackrel{~}{\rho }.`$
Hence we again need to find expressions for $`\stackrel{~}{\rho }_{ee}`$ and $`\stackrel{~}{\rho }_{ge}`$ in terms of $`\stackrel{~}{\rho }`$. We can achieve this by noting that, as in the earlier treatment, $`\stackrel{~}{\rho }_{ee}`$ and $`\stackrel{~}{\rho }_{ge}`$ equilibrate on a timescale much faster than $`\stackrel{~}{\rho }_{gg}`$. By setting $`\dot{\stackrel{~}{\rho }}_{ge}=0`$, and again ignoring the kinetic energy term, we get

$`\stackrel{~}{\rho }_{ge}`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{\rho }\mathrm{\Omega }^{}}{2\delta }}{\displaystyle \frac{[\mathrm{\Omega }\mathrm{\Omega }^{},\stackrel{~}{\rho }\mathrm{\Omega }^{}]}{8\delta ^3}}+{\displaystyle \frac{\{\mathrm{\Omega }^{},\stackrel{~}{\rho }_{ee}\}}{2\delta }}{\displaystyle \frac{\mathrm{\Gamma }\stackrel{~}{\rho }\mathrm{\Omega }^{}}{4i\delta ^2}}`$ (21)
$`+{\displaystyle \frac{\mathrm{\Gamma }\{\mathrm{\Omega }^{},\stackrel{~}{\rho }_{ee}\}}{4i\delta ^2}}+{\displaystyle \frac{[\mathrm{\Omega }\mathrm{\Omega }^{},\{\mathrm{\Omega }^{},\stackrel{~}{\rho }_{ee}\}]}{8\delta ^3}}.`$
Setting $`\dot{\stackrel{~}{\rho }}_{ee}=0`$ and ignoring the kinetic energy term allows us to solve for $`\stackrel{~}{\rho }_{ee}`$. In the more sophisticated adiabatic treatment, instead of getting an equation of the form of Eq. (12), we get

$$\stackrel{~}{\rho }_{ee}=\frac{\mathrm{\Omega }\stackrel{~}{\rho }\mathrm{\Omega }^{}}{4\delta ^2}.$$
(22)
In Eq. (22), the superoperators $``$ and $``$ are defined

$``$ $``$ $`1+{\displaystyle \frac{i}{2\delta \mathrm{\Gamma }}}[\mathrm{\Omega }\mathrm{\Omega }^{},\mathrm{}],`$ (23)
$``$ $``$ $`+\text{ higher order terms}.`$ (24)
In the regime chosen ($`\mathrm{\Gamma }\mathrm{\Omega }^2/\delta `$), the second term of Eq. (23) is of order 1 and so cannot be ignored compared to the leading term. In Eq. (24), the higher order terms are of order much smaller than 1 and even smaller than $`\mathrm{\Omega }^4/\delta ^3`$, which is the smallest order kept by making Taylor approximations to get the expression for $`\stackrel{~}{\rho }_{ge}`$. Thus to simplify Eq. (22), we act on the left of both sides with the superoperator $`^1`$. Then to leading order we have

$`\stackrel{~}{\rho }_{ee}`$ $``$ $`{\displaystyle \frac{\mathrm{\Omega }\stackrel{~}{\rho }\mathrm{\Omega }^{}}{4\delta ^2}}.`$ (25)
As expected, the interaction picture chosen produced a term in $`\rho _{gg}`$ to counteract the unwanted term in $`\rho _{ee}`$. Notice here that we come to essentially the same result as in the standard adiabatic treatment in Eq. (13) but without making any unjustified approximations.

Now substituting Eq. (21) and Eq. (25) into Eq. (20), after simplification, leads to the final interaction picture master equation

$`\dot{\stackrel{~}{\rho }}`$ $`=`$ $`\mathrm{\Gamma }\left(𝒥\left[{\displaystyle \frac{\mathrm{\Omega }}{2\delta }}\right]\stackrel{~}{\rho }𝒜\left[{\displaystyle \frac{\mathrm{\Omega }}{2\delta }}\right]\stackrel{~}{\rho }\right)i[{\displaystyle \frac{\mathrm{\Omega }^2\mathrm{\Omega }_{}^{}{}_{}{}^{2}}{16\delta ^3}},\stackrel{~}{\rho }]`$ (26)
$`\stackrel{~}{𝒦}\stackrel{~}{\rho }.`$
Finally, all that remains is to transform the interaction picture master equation back to the Schrödinger picture by performing the opposite unitary transformation. This leaves the final master equation for the more sophisticated adiabatic elimination treatment as

$`\dot{\rho }`$ $`=`$ $`\mathrm{\Gamma }\left(𝒥\left[{\displaystyle \frac{\mathrm{\Omega }}{2\delta }}\right]\rho 𝒜\left[{\displaystyle \frac{\mathrm{\Omega }}{2\delta }}\right]\rho \right)`$ (27)
$`i[{\displaystyle \frac{p_x^2}{2m}}{\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{}}{4\delta }}+{\displaystyle \frac{\mathrm{\Omega }^2\mathrm{\Omega }_{}^{}{}_{}{}^{2}}{16\delta ^3}},\rho ].`$
Notice here that by explicitly accounting for the term in $`\mathrm{\Omega }^2/\mathrm{\Gamma }\delta `$, the master equation derived is the same as that of the standard adiabatic approach with an extra potential term of order $`\mathrm{\Omega }^4/\delta ^3`$. This term is very small and as such the statements made to justify the standard approach Hen01a  were correct in that the adjustment to the final master equation is small. The more sophisticated treatment however, although more algebraically intensive to produce the initial master equation, requires very little extra effort than the standard adiabatic treatment to simulate. More importantly though, the master equation derived in the more sophisticated adiabatic approach is valid in the regime $`\mathrm{\Omega }_{\mathrm{max}}^2/\delta \mathrm{\Gamma }`$.

## IV A Secular Approximation

A secular approximation to the full master equation is quite different from any adiabatic approximation. A secular approximation does not totally remove any dependence on the internal state of the atom. It only eliminates the coherences between the internal atomic states. The secular approximation was also used by Hensinger et al Hen01a , alongside the standard adiabatic approximation.

There are a number of ways to derive a secular approximation to the full master equation. One such method has been performed by Dyrting and Milburn dyrtingmilburn  but this method is quite complicated. A much simpler method which produces the same approximate master equation is based on the technique that is applied for the standard adiabatic approximation.

We start with the Eqs. (7) for $`\rho _{gg}`$ etc, and then adiabatically eliminate $`\rho _{eg}`$ as in Sec. II. However, instead of also trying to solve for $`\rho _{ee}`$, we just substitute the expression for $`\rho _{ge}`$, Eq. (10), into the equations for $`\rho _{gg}`$, Eq. (7a) and $`\rho _{ee}`$, Eq. (7c). This eliminates the coherences between the excited and ground states while still keeping much of the original master equation. One advantage of this is that it allows comparison to an analogous 2-state classical model dyrtingmilburn . More importantly to us, this approximation eliminates the evolution at rate $`\delta `$ as is necessary to simplify the simulations. We are left with the following equations for the ground and excited state density matrices

$`\dot{\rho }_{gg}`$ $`=`$ $`\mathrm{\Gamma }\rho _{ee}+\mathrm{\Gamma }\left(𝒥\left[{\displaystyle \frac{\mathrm{\Omega }^{}}{2\delta }}\right]\rho _{ee}𝒜\left[{\displaystyle \frac{\mathrm{\Omega }}{2\delta }}\right]\rho _{gg}\right)`$ (28)
$`+i[{\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{}}{4\delta }},\rho _{gg}]𝒦\rho _{gg},`$
$`\dot{\rho }_{ee}`$ $`=`$ $`\mathrm{\Gamma }\rho _{ee}+\mathrm{\Gamma }\left(𝒥\left[{\displaystyle \frac{\mathrm{\Omega }}{2\delta }}\right]\rho _{gg}𝒜\left[{\displaystyle \frac{\mathrm{\Omega }^{}}{2\delta }}\right]\rho _{ee}\right)`$ (29)
$`i[{\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{}}{4\delta }},\rho _{ee}]𝒦\rho _{ee}.`$
The final master equation is constructed by recombining the equations (28) and (29) to produce an equation which reproduces these equations for $`\dot{\rho }_{gg}`$ and $`\dot{\rho }_{ee}`$ while also only giving rapidly decaying terms for $`\dot{\rho }_{ge}`$ and $`\dot{\rho }_{eg}`$. This final master equation for the secular approximation is

$`\dot{\rho }`$ $`=`$ $`i[{\displaystyle \frac{p_x^2}{2m}}+{\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{}}{4\delta }}\sigma _z,\rho ]+\mathrm{\Gamma }\left(𝒥\left[\sigma \right]\rho 𝒜\left[\sigma \right]\rho \right)`$ (30)
$`+\mathrm{\Gamma }\left(𝒥\left[{\displaystyle \frac{\mathrm{\Omega }^{}\sigma }{2\delta }}\right]\rho 𝒜\left[{\displaystyle \frac{\mathrm{\Omega }^{}\sigma }{2\delta }}\right]\rho \right)`$
$`+\mathrm{\Gamma }\left(𝒥\left[{\displaystyle \frac{\mathrm{\Omega }\sigma ^{}}{2\delta }}\right]\rho 𝒜\left[{\displaystyle \frac{\mathrm{\Omega }\sigma ^{}}{2\delta }}\right]\rho \right),`$
where $`\sigma _z`$ is just the Pauli spin operator $`\sigma ^{}\sigma \sigma \sigma ^{}`$.

To compare this master equation with those we have already seen, the Hamiltonian terms derived here are the same as those derived in the standard adiabatic approximation, Eq. (LABEL:adiab1me). The spontaneous emission term exactly as in the full master equation, Eq. (2b), remains, while two extra jump terms involving a state change with no spontaneous emission have been derived.

There are no apparent problems in this derivation, or that in dyrtingmilburn . Nevertheless, as we will discuss in Sec VIII, Eq. (30) does not give accurate results in comparison with Eq. (1), in the regime $`\mathrm{\Omega }^2\mathrm{\Gamma }\delta `$.

## V A Dressed-State Approximation

A semi-classical dressed-state treatment of atomic motion was put forward by Dalibard and Cohen-Tannoudji Dalibard . The dressed-state approximation used here is a fully quantum version of that treatment.

The states we have been using for a basis so far, $`|g`$ and $`|e`$, are called bare states. We can also work with another basis of position-dependent states which we call dressed-states. These dressed-states are derived by considering the Hamiltonian of the full master equation Eq. (2a) and ignoring the kinetic energy component to get

$$V=\delta \sigma ^{}\sigma +\frac{1}{2}\left(\mathrm{\Omega }\sigma ^{}+\mathrm{\Omega }^{}\sigma \right).$$
(31)
The diagonalization of the Hamiltonian yields the form

$$V=E_1(x,t)a^{}a+E_2(x,t)aa^{},$$
(32)
where $`E_1`$ and $`E_2`$ are the eigenenergies of $`V`$. The corresponding eigenstates are the position dependent dressed-states we will be considering here, labeled $`|1`$ and $`|2`$. In terms of these dressed-states, $`a`$ is given by $`|21|`$. The energies and states are

$`E_1(x,t)`$ $`=`$ $`\frac{1}{2}\left(\delta +\mathrm{\Delta }\right),`$ (33)
$`E_2(x,t)`$ $`=`$ $`\frac{1}{2}\left(\delta \mathrm{\Delta }\right),`$ (34)
$`|1`$ $`=`$ $`\mathrm{sin}\theta |g+\mathrm{cos}\theta |e,`$ (35)
$`|2`$ $`=`$ $`\mathrm{cos}\theta |g\mathrm{sin}\theta |e,`$ (36)
where $`\theta `$ is defined by

$`\mathrm{sin}\theta `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{\sqrt{2}}}\left(\delta ^2+\delta \mathrm{\Delta }+\mathrm{\Omega }\mathrm{\Omega }^{}\right)^{\frac{1}{2}},`$ (37a)
$`\mathrm{cos}\theta `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }+\delta }{\sqrt{2}}}\left(\delta ^2+\delta \mathrm{\Delta }+\mathrm{\Omega }\mathrm{\Omega }^{}\right)^{\frac{1}{2}},`$ (37b)
and $`\mathrm{\Delta }=\left(\delta ^2+\mathrm{\Omega }\mathrm{\Omega }^{}\right)^{\frac{1}{2}}.`$

Now, we use this basis to get equations for $`\rho _{11}`$, $`\rho _{12}`$, $`\rho _{21}`$ and $`\rho _{22}`$. Without including the spontaneous emission or the kinetic energy, these equations are simply

$$\dot{\rho }_{jk}=i\left[E_j(x,t)\rho _{jk}\rho _{jk}E_k(x,t)\right].$$
(38)
We do however want to include the spontaneous emission in our treatment and so we must assess the effect of the raising and lowering operators for the bare states ($`\sigma =|ge|`$ and $`\sigma ^{}=|eg|`$) acting on the new basis states ($`|1`$ and $`|2`$). This is easily done.

So far there have been no approximations made. The essence of the dressed-state approximation is similar to the secular approximation in that we keep the internal dressed-states but allow the coherences to go to zero. In contrast to the secular approximation, we do not simply set $`\dot{\rho }_{12}=0`$. Instead, we notice from Eq. (38) that, if we ignore the operator nature of the eigenenergies, then the equations for $`\rho _{12}`$ would have a term of the form $`i\left(E_1E_2\right)\rho _{12}`$. To first order, $`E_1E_2`$ is approximately $`\delta `$. Thus $`\rho _{12}`$ will rotate very quickly such that only terms rotating at this very rapid pace will be able to contribute to its evolution. Thus only terms involving $`\rho _{12}`$ are kept in the equation for $`\rho _{12}`$

$`\dot{\rho }_{12}`$ $`=`$ $`\mathrm{\Gamma }\left(\mathrm{sin}\theta \mathrm{cos}\theta \rho _{12}\mathrm{cos}\theta \mathrm{sin}\theta \right)`$ (39)
$`{\displaystyle \frac{\mathrm{\Gamma }}{2}}\left(\mathrm{cos}^2\theta \rho _{12}+\rho _{12}\mathrm{sin}^2\theta \right)`$
$`i\left[E_1(x,t)\rho _{12}\rho _{12}E_2(x,t)\right].`$
Knowing that the last term in Eq. (39) serves to force $`\rho _{12}`$ to oscillate very rapidly and the other two terms force it to decay quickly, $`\rho _{12}`$ will average to zero. Thus we set $`\rho _{12}`$ and $`\rho _{21}`$ equal to zero in the population equations giving

$`\dot{\rho }_{11}`$ $`=`$ $`\mathrm{\Gamma }\left(\mathrm{sin}\theta \mathrm{cos}\theta \rho _{11}\mathrm{sin}\theta \mathrm{cos}\theta +\mathrm{sin}^2\theta \rho _{22}\mathrm{sin}^2\theta \right)`$ (40)
$`{\displaystyle \frac{\mathrm{\Gamma }}{2}}\{\mathrm{cos}^2\theta ,\rho _{11}\}i[E_1(x,t),\rho _{11}],`$
$`\dot{\rho }_{22}`$ $`=`$ $`\mathrm{\Gamma }\left(\mathrm{sin}\theta \mathrm{cos}\theta \rho _{22}\mathrm{sin}\theta \mathrm{cos}\theta +\mathrm{cos}^2\theta \rho _{11}\mathrm{cos}^2\theta \right)`$ (41)
$`{\displaystyle \frac{\mathrm{\Gamma }}{2}}\{\mathrm{sin}^2\theta ,\rho _{22}\}i[E_2(x,t),\rho _{22}].`$
If we could ignore the operator nature of the rates in Eqs. (39-41), we would obtain the same rates as given by Dalibard and Cohen-Tannoudji in Dalibard .

All that remains now is to design the master equation that forces $`\rho _{12}`$ to zero giving these population equations. The master equation which fulfills these requirements is

$`\dot{\rho }`$ $`=`$ $`𝒦\rho i[E_1(x,t)|11|+E_2(x,t)|22|,\rho ]`$ (42)
$`+\mathrm{\Gamma }𝒥\left[\mathrm{cos}\theta \mathrm{sin}\theta \left(a^{}aaa^{}\right)\right]\rho `$
$`\mathrm{\Gamma }𝒜\left[\mathrm{cos}\theta \mathrm{sin}\theta \left(a^{}aaa^{}\right)\right]\rho `$
$`+\mathrm{\Gamma }\left(𝒥\left[\mathrm{cos}^2\theta a\right]\rho 𝒜\left[\mathrm{cos}^2\theta a\right]\rho \right)`$
$`+\mathrm{\Gamma }\left(𝒥\left[\mathrm{sin}^2\theta a^{}\right]\rho 𝒜\left[\mathrm{sin}^2\theta a^{}\right]\rho \right),`$
where the kinetic energy term has been restored. This form of the master equation is still very complicated remembering the definitions of $`\mathrm{cos}\theta `$ and $`\mathrm{sin}\theta `$ from Eqs. (37). It is possible to simulate this master equation as it stands but the simulation would be very slow and inefficient. This master equation can, however, be approximated further with almost no loss in accuracy.

The definitions in Eqs. (37) can be approximated by remembering that we are working in the regime of $`\delta \mathrm{\Omega }`$. Thus, to leading order, $`\mathrm{cos}\theta `$ is simply 1, and $`\mathrm{sin}\theta `$ is $`\mathrm{\Omega }/2\delta `$. Also any terms of order $`\mathrm{sin}^4\theta `$ are extremely small and so are also ignored. If we propagate these approximations through our system, we find that the dressed-state $`|1`$ is very close to the bare excited state $`|e`$. Also, the dressed-state $`|2`$ is very near the bare ground state $`|g`$. Thus if we make the approximations

$$|1|e\text{ , }|2|g,$$
(43)
then we can similarly approximate

$$\sigma a,\sigma ^{}a^{}.$$
(44)
The last approximation is to expand the dressed-state eigenenergies in a Taylor series to second order giving

$`E_1(x,t)`$ $``$ $`\delta +{\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{}}{4\delta }}{\displaystyle \frac{\mathrm{\Omega }^2\mathrm{\Omega }_{}^{}{}_{}{}^{2}}{16\delta ^3}},`$ (45a)
$`E_2(x,t)`$ $``$ $`\left({\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{}}{4\delta }}{\displaystyle \frac{\mathrm{\Omega }^2\mathrm{\Omega }_{}^{}{}_{}{}^{2}}{16\delta ^3}}\right).`$ (45b)
This leaves us with the final master equation for the dressed-state approximation as

$`\dot{\rho }`$ $`=`$ $`i[{\displaystyle \frac{p_x^2}{2m}}+\left({\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{}}{4\delta }}{\displaystyle \frac{\mathrm{\Omega }^2\mathrm{\Omega }_{}^{}{}_{}{}^{2}}{16\delta ^3}}\right)\sigma _z,\rho ]`$ (46)
$`+\mathrm{\Gamma }\left(𝒥\left[\sigma \right]\rho 𝒜\left[\sigma \right]\rho \right)`$
$`+\mathrm{\Gamma }\left(𝒥\left[{\displaystyle \frac{\mathrm{\Omega }}{2\delta }}\sigma _z\right]\rho 𝒜\left[{\displaystyle \frac{\mathrm{\Omega }}{2\delta }}\sigma _z\right]\rho \right).`$
Here we have removed the $`H_0=\delta \sigma ^{}\sigma `$ Hamiltonian term by moving into an interaction picture.

Again comparing this master equation to those already seen, the Hamiltonian terms here are exactly analogous to those derived in the sophisticated adiabatic approximation Eq. (27). Again we have the same spontaneous emission term as in Eq. (2b), but this time, the higher order correction term includes a spontaneous emission without changing the internal state of the atom, the opposite of the case in the secular approximation.

This would be a useful approximate master equation, especially if the atom was initially in the excited state. In our case however, we can simplify this equation further by noting that the jump terms keep the atom in the ground state. Once the excited state populations are reduced to zero, this master equation reduces to exactly the same master equation as derived for the more sophisticated adiabatic approximation Eq. (27). The dressed state master equation results would thus lie exactly on top of those of the more sophisticated adiabatic approximation and as such are not included in the simulations.

## VI Simulation Of The Master Equations

Now that we have derived the forms of the master equations which we wish to compare, we need to set up a method of simulation for the different approximations and the full master equation. In this case, the numerical environment MATLAB turned out to be the most useful tool, combined with the Quantum Optics Toolbox produced by S.M. Tan Sze1 ; Sze2 . The simulation is designed by converting states and operators to vectors and matrices. Making this conversion requires a number of different adaptations of the theoretical master equations.

Firstly, we need to chose a form for the complex Rabi frequency operator, $`\mathrm{\Omega }`$. The form of the Rabi frequency operator we use here is $`\mathrm{\Omega }=\mathrm{\Omega }_{\mathrm{max}}\mathrm{sin}kx`$, such that there is no time dependence (standing wave) and the operator is Hermitian. Rewriting the sine function in terms of exponentials, $`e^{\pm ikx}`$, we can examine the action of these exponentials. This suggests using the momentum representation, because the exponentials, $`e^{\pm ikx}`$, simply represent single momentum kicks to the atom of $`1\mathrm{}k`$, or $`1`$ atomic unit of momentum.

Having chosen to use the momentum representation to evaluate our master equation solutions, we need to address concerns with the conversion process. Firstly, the momentum range must be truncated. The loss of probability due to this truncation should be kept below some small level, say $`10^5`$. Using this rule of thumb, it was found that a momentum Hilbert space ranging from $`25\mathrm{}k`$ up to $`+25\mathrm{}k`$ was sufficient for the parameters we chose (discussed later).

We still face a limitation problem in that the momentum Hilbert space is continuous. To simulate this on a computer, we need to discretize this space. We are able to discretize this space, again as a consequence of the choice of $`\mathrm{\Omega }`$. The exponential component nature of $`\mathrm{\Omega }`$ provides for momentum kicks of exactly one atomic unit of momentum in the direction of propagation of the laser light (the $`x`$-direction). Thus the only momentum kicks that are not in a single unit of atomic momentum in the $`x`$-direction are due to spontaneous emission of photons from the atom. This allows for a momentum kick of one atomic unit in any direction, which, when projected onto the $`x`$-direction, allows for a random kick of anywhere between $`1`$ and $`1`$ atomic unit. However the relative infrequency of the spontaneous emission allows us to approximate this by a kick of $`1`$, $`0`$, or $`1`$.

This approximation requires us to convert the integral over the atomic dipole distribution to a discrete sum. Following Ref. Hen01a , we let

$$_1^1W(u)𝑑u\underset{u=1,0,1}{}V(u),$$
(47)
where the discrete function $`V(u)`$ is obtained by making sure that it has the same normalization, mean momentum kick and mean squared momentum kick. This gives $`V(1)=\frac{1}{5}`$, $`V(0)=\frac{3}{5}`$ and $`V(1)=\frac{1}{5}`$.

We now need to define our matrix notation. The simulations here will set $`\mathrm{}`$, $`k`$ and $`m`$ to 1 so that the kinetic energy operator can be represented as

$$\frac{p_x^2}{2m}\frac{1}{2}\underset{n=25}{\overset{25}{}}n^2|nn|,$$
(48)
where $`|n`$ denotes a momentum eigenstate. The $`e^{ikx}`$ operator in a momentum Hilbert space just acts as a raising operator given by

$$R=\underset{n=25}{\overset{25}{}}|n+1n|.$$
(49)
The $`e^{ikx}`$ operator is just the lowering operator $`R^{}`$. For simplicity, we take the initial density matrix to be the zero momentum state $`|0`$, given by a $`51\times 51`$ matrix of zeros with a 1 in the very centre.

These are all the approximations required to allow the simulation of the adiabatic approximations. The full master equation and the secular and dressed-state approximations require the internal state information as well. This is achieved by constructing the tensor product of the $`51\times 51`$ momentum Hilbert with a $`2\times 2`$ internal state Hilbert space. This is easily done using the quantum optics toolbox.

The only problem remaining is to adopt a method of simulation. There are a number of different simulation methods available as part of the quantum optics toolbox as well as any number of methods available using regular ordinary differential equation techniques. The method we use here is a built-in function of the Quantum Optics Toolbox called odesolve. This allows options to be specified for use on both smooth and stiff ODE problems. The simulations presented here using odesolve were checked using first, a hybrid Euler and matrix exponential method, and second, a modified mid-point method combined with Richardson extrapolation from Press et al Press . The results all agreed well but the odesolve method was by far the fastest and easiest to use.

We wish to simulate the experimental regime of Hensinger et al, where $`\delta \mathrm{\Omega }_{\mathrm{max}}\mathrm{\Gamma }\mathrm{\Omega }_{\mathrm{max}}^2/\delta `$. However, the actual experimental parameters would be prohibitively time-consuming to simulate. This is both because of the separation in time scales between the fastest ($`\delta `$) and slowest dynamics ($`\mathrm{\Gamma }\mathrm{\Omega }_{\mathrm{max}}`$), and also because of the basis size required. The latter is determined by the fact that $`p_{\mathrm{max}}^2/2m`$ must be larger than the effective potential drop, $`\mathrm{\Omega }_{\mathrm{max}}^2/4\delta `$. If we were to use the parameters of the experiment of Hen01a  then we would need a basis size of more than $`\pm 70\mathrm{}k`$. However, we can scale the parameters down and still preserve the regime of the experiments.

As well as working in the regime $`\delta \mathrm{\Omega }_{\mathrm{max}}\mathrm{\Gamma }\mathrm{\Omega }_{\mathrm{max}}^2/\delta `$, we require for the validity of the adiabatic approximations that $`\mathrm{\Gamma }`$ be much larger than the oscillation frequency near the potential minimum, $`\omega _{\mathrm{osc}}`$. In our scaled units ($`\mathrm{}=k=m=1`$), the latter is of order $`\sqrt{\mathrm{\Omega }_{\mathrm{max}}^2/\delta }`$. On this basis, we have chosen paramaters of $`\delta =10^4`$, $`\mathrm{\Omega }_{\mathrm{max}}=2\times 10^3`$ and $`\mathrm{\Gamma }=200`$, leaving $`\mathrm{\Omega }_{\mathrm{max}}^2/\delta =400`$, giving $`\omega _{\mathrm{osc}}=20`$. For Rb as in Ref. Hen01b , we have $`k=2\pi /780`$nm and $`m=1.42\times 10^{25}`$kg, so in SI units, the frequency unit is $`\mathrm{}k^2/m=4.821\times 10^4`$ s<sup>-1</sup>. Note that the $`\mathrm{\Gamma }`$ we have chosen is not the true radiative decay rate for Rb.

The scaled time unit can be given meaning by examining the spontaneous emission rate. For each of the approximations, the fluorescence rate is (to leading order) $`\mathrm{\Gamma }\overline{\mathrm{\Omega }^2}/4\delta ^2`$. Here $`\overline{\mathrm{\Omega }^2}`$ is a time-averaged effective value, which would be somewhat less than $`\mathrm{\Omega }_{\mathrm{max}}^2`$. This means that after a time period of $`4\delta ^2/\mathrm{\Gamma }\mathrm{\Omega }_{\mathrm{max}}^2`$, we would expect there to have somewhat under one spontaneous emission. This time is 2 scaled time units. There is, however, a lot of evolution occurring in that time period. To compare the approximations, we look in detail at a period from 0 to 2 time units, and also look at the long time results at 8 time units.

## VII Results of the simulations

In comparing the results of the simulations, we first compare the accuracy of the approximations, and then the resources required to perform the calculations.

### VII.1 Accuracy

To compare the accuracy of the simulations, we look at the momentum distribution of the atom as it evolves through time. The interesting components of this evolution are the probability to have 0 momentum and the probability to have 1 atomic unit of momentum as the atom evolves through time. The other probabilities evolve similarly to one of these two.

Firstly, we will look at the probability to have zero momentum over a relatively short timescale. The approximations are so close to the full master equation that, at full size, they are almost impossible to distinguish from the full master equation. Fig. 1a shows an overall picture of how the probability to have zero momentum evolves through time.

Fig. 1b zooms in on a section of the full size figure to illustrate the differences between the approximations. As we can see, the secular and standard adiabatic approximation evolutions both significantly lead that of the full master equation in time. The dressed-state and sophisticated adiabatic approximations are very close to the full master equation solution even at this magnification. Fig 1c zooms in even closer to try to distinguish the sophisticated adiabatic approximation from the full master equation.

This trend continues as we investigate the probability to have 1 atomic unit of momentum in Fig. 2.

The inaccuracy of the secular and standard adiabatic approximations are again evident in Fig. 2b. The more sophisticated adiabatic approximation is very close to the full master equation evolution and is still indistinguishable at this magnification.

Fig. 2c provides a means of comparing the sophisticated adiabatic approximation to the full master equation evolution in detail. As we can see, the more sophisticated adiabatic approximation is again very close to the full master equation.

The last visual comparison to make is to see how the evolution described by the approximations matches that of the full master equation after a very long time period. The time period that has elapsed in Fig. 3 is 8 time units, after which, we would have expected there to be a number of spontaneous emissions. At this point in time, we compare the probability densities described by the approximations and the full master equation.

It is evident from Fig. 3 that the standard adiabatic and the secular approximations are quite poor methods for simulating the full master equation evolution over a long time period. The main reason for this is the leading behaviour such as we see in Fig. 1a. Even at long times, the probability to have zero momentum is still oscillating and the leading behaviour means that the standard adiabatic and secular approximations are not oscillating in phase with the full master equation. Thus, even though they follow roughly the same shape, they are not necessarily at the same point in the oscillation as the full master equation. What is even more striking though is that the secular approximation evolution seems to follow a slightly different trend for the probabilities to have larger momentum. The intriguing features of the secular approximation are discussed in Sec VIII.

To examine how close the sophisticated adiabatic and the dressed-state approximations are to the full master equation, Fig. 4 focuses on a smaller section to provide a comparison.

It is evident from Fig. 4 that the dressed-state and sophisticated adiabatic approximations are very close to the full master equation, even at this long time.

### VII.2 Resources

The main resource required to perform the simulations is time. Although each approximation has different memory and processing power requirements, these needs are reasonably accurately reflected in the time each simulation takes to run. The times quoted in table 1 are the times required to calculate the set of results from 0 to 8 time units and are quoted in seconds.

As we can see from the results in table 1, the adiabatic approximations are both clearly the fastest. The secular approximation is just over 4 times larger than the adiabatic approximations. This is not entirely unexpected. We would have expected at least a doubling in time by using the methods that included state information. The secular approximation is supposed to force the coherences to zero. Unfortunately, using our method to simulate this allows zero to be anywhere up to $`10^5`$. This, although small, still has to be processed explaining the four-fold increase in time. One reason it is just over 4 times the time required for the adiabatic approximations could be that after the evolution has been calculated, there is still a partial trace to be performed to obtain a solution of the same form as the adiabatic approximations. We could have limited this time by simulating two coupled equations instead of a full master equation and then would probably have only doubled the time taken.

## VIII Discussion

The standard and more sophisticated adiabatic approximations take similar times to simulate, and are much faster than the full master equation. The only difference between them is an extra potential term. It turns out, though, that this Hamiltonian term is quite important in accurately describing the motion of the atom. While the standard approach evolves too quickly and leads that of the full master equation evolution, the more sophisticated approach with the modified potential does not suffer this problem. The evolution described by this more sophisticated approach is very close to the full master equation, even at long times. Of course the dressed-state approximation offers the same accuracy as the sophisticated adiabatic approximation which may be useful if we wanted to simulate an initially excited atom.

These successes contrast the results from the secular approximation. As one can see from Figs. 1 and 2, the secular approximation not only leads the full master equation solution but it also predicts a lower probability to have either zero or one atomic unit of momentum. This result is surprising because the secular approximation master equation is quite similar to the others. Investigating this further, we find that the secular approximation master equation simulation shows that the probability to have 25 atomic units of momentum increases exponentially much faster than any of the other approximations. To analyse this in another manner, we find that for the secular approximation Tr$`[\rho ]`$ falls off from 1 exponentially much faster than any of the other approximations or the full master equation.

If we consider the dressed-state approximation master equation, Eq. (46), we notice that the Hamiltonian terms are, to leading order, the same as the final secular approximation master equation Eq. (30). The only advantage the dressed state approximation has in terms of it’s Hamiltonian is the presence of the higher order term. The leading order Lindbladian terms are also identical with only the higher order terms differing. The dressed-state master equation includes Lindbladian terms which give rise to a momentum kick to the atom (from the $``$ superoperator) without a change in the internal state of the atom. The higher order secular approximation master equation Lindbladian term provides the opposite. Here there is an internal state change without a momentum kick to the atom. Thus it is quite unexpected that the secular approximation solution should differ so greatly from the dressed-state approximation solution as they are identical in the leading order terms.

It is a fairly simple matter to investigate the effect of the differences between the two approximations by simulating only parts of the master equations. If we only included the Hamiltonian terms from each master equation, the difference would only be the higher order term in the dressed state master equation. This term has the same effect here as it does for the sophisticated adiabatic approximation, correcting the leading behaviour. It does not however account for the secular approximation predicting lower probabilities as it does in Figs. 1 and 2. Thus to investigate the difference between the higher order Lindbladian terms, we drop the $``$ superoperator from the dressed-state master equation and add it to the secular approximation master equation. This does not correct the problems with the secular approximation nor does it severely affect the dressed-state approximation. This only leaves the inherent difference that the dressed-state master equation provides a correction term without requiring a state change where the secular approximation forces a state change in it’s correction term. Thus we have to conclude that the best approximation to the full master equation involves a Lindbladian term that does not change the internal state of the atom. Lacking it, the secular approximation is the poorest.

One other point to note in this discussion is one of the conditions on which the adiabatic elimination techniques are based. This is that $`\mathrm{\Gamma }\rho _{ee}`$ be much larger than $`𝒦\rho _{ee}`$. This will be satisfied if

$$\mathrm{Tr}\left[\rho _{ee}\frac{p_x^2}{2m}\right]/\mathrm{\Gamma }\mathrm{Tr}\left[\rho _{ee}\right]1.$$
(50)
For the parameter regime of this investigation, it is not obvious that this holds. Fig 5 shows how the ratio in Eq. (50) evolves through time for the sophisticated adiabatic approximation.

Here we see that numerically, this ratio is around 0.1 which is of the same order as the other ratios (such as $`\mathrm{\Omega }_{\mathrm{max}}/\delta `$) which are required to be small for our approximations.

Finally, we discuss the possibility of simulating with the true experimental parameters. This is difficult because of the stiffness of the full master equation, and the basis size required for all methods. The latter problem can be avoided by using quantum trajectories MolCasDal93 ; DumZolRit92 ; Car93 . Hensinger et al Hen01a  actually use quantum trajectory simulations based on the secular approximation master equation. It is however, possible to convert any master equation of the Lindblad form to a quantum trajectory simulation Wiseman . All of the approximate master equations we have developed here have been written in the Lindblad form and as such all of these could be simulated using quantum trajectories.

## IX conclusion

There are a number of theoretical models for the motion of an atom as it interacts with a light field. This paper has investigated the possibility of using four different approximations as opposed to using the full master equation to simulate an experimental system. Two have been widely used in the past and two have not. We have given a detailed explanation of the mathematical principles to perform each of these approximations on a fairly general system. We have also compared them numerically to to the true dynamics from the full master equation.

In a regime of particular experimental interest, we have found that the most accurate results are obtained from two approaches that we have introduced here, a sophisticated adiabatic approach and a dressed-state approach. These give identical equations in the regime of interest, and in terms of resources, they are almost as fast to simulate as the standard adiabatic approximation. This has been most used in the past, but deviates significantly from the true dynamics for long times. The other approximation that has been used in the past, the secular approximation, is even poorer. On top of the failings of the standard adiabatic approximation, it takes longer to simulate and appears to produce anomalous momentum diffusion.

###### Acknowledgements.

We gratefully acknowledge discussions with W. Hensinger which prompted us to pursue the investigation. This work was supported by the Australian Research Council.