text:
# First Order Transition in the Ginzburg-Landau Model

\[

## Abstract

The $`d`$-dimensional complex Ginzburg-Landau (GL) model is solved according to a variational method by separating phase and amplitude. The GL transition becomes first order for high superfluid density because of effects of phase fluctuations. We discuss its origin with various arguments showing that, in particular for $`d=3`$, the validity of our approach lies precisely in the first order domain.
pacs: 05.70.Fh, 72.40.+k, 64.60.-i

\]

Ginzburg-Landau free energy functionals involving a $`n`$-component space and time dependent field have been widely used in order to describe different types of phase transitions in a semi-phenomenological way. The case $`n=2`$, corresponding to a complex field $`\psi `$, applies in particular to superconductivity, superfluidity, metal-insulator transitions or to magnetic systems with moments that are confined to a plane (XY moments). For such particular applications, the coefficients determining the functional can be derived from appropriate microscopic models. For superconductors, this was originally done by Gorkov  and has been refined since by numerous authors . In its simplest form, the GL functional involves a time independent field and thus describes (classical) thermodynamic and static phenomena of superconductors. For the description of dynamic phenomena, charging effects  or pairing fluctuations in strong coupling superconductors , one needs the generalisation to a time dependent field.

An important aspect of the case $`n=2`$ is the interplay between variations of amplitude $`|\psi |`$ and phase $`\varphi `$ of the corresponding complex field $`\psi =|\psi |e^{i\varphi }`$. In various approximate treatments, such as the mean field approximation or the Hartree decoupling of terms involving higher powers of $`\psi `$, the field is treated as a whole, without separating amplitude and phase. More accurate studies of the static GL problem, like the renormalisation group approach  focusing in particular on the region near the phase transition, show that the amplitude has no critical behaviour and is irrelevant at the transition. The phase transition scenario should then correspond to the one of the XY model with the same dimensionality. In the framework of the $`\epsilon `$-expansion , for $`d=3`$, the transition is second order and seems to have the same critical exponents as the XY model. On the other hand, Bormann and Beck  have shown that amplitude fluctuations, even though not being critical by themselves, might alter the cooperative phenomenon occurring with phases, in particular in dimension 2. Like the XY model, corresponding to a fixed value of the amplitude, the 2$`d`$GL model can be mapped onto a Coulomb gas describing vortex-antivortex pairs. As soon as one allows for amplitude variations, these topological excitations become energetically more favorable. Taking into account gaussian amplitude fluctuations, Bormann and Beck  have shown that the system may be driven into a regime where - according to Minnhagen’s phase diagram  \- a first order transition replaces the usual Kosterlitz-Thouless scenario.

As far as superconductors are concerned, BCS theory  predicts that the transition between the normal state and the superconducting state is a second order phase transition. However, it is well known that fluctuations can change the order of the transition. For example, fluctuations of the magnetic field change the GL-BCS transition to a first order transition for type I superconductors . The three state Potts model in two dimensions is an opposite example: mean field theory predicts a first order transition, whereas the actual transition is continuous. So the question of the order of the transition in the GL model, as well as the more detailed mutual influence between phase and amplitude, is still an open problem.

The aim of this letter is to show the reciprocal influence between phase and amplitude by separating self-consistently, from the outset, the GL functional into two parts: the amplitude part and the phase part.

According to Ginzburg-Landau theory, we define the effective hamiltonian functional

$$H[\psi ]=d^dr\left[at\left|\psi \right|^2+\frac{b}{2}\left|\psi \right|^4+\frac{\gamma }{2}\left|\psi \right|^2\right]$$
(1)
where $`a`$, $`b`$ and $`\gamma `$ are coefficients independent of the temperature derived from a microscopic model. $`t=T/T_01`$ is the reduced temperature and $`T_0`$ is the mean field critical temperature. We now introduce the amplitude $`|\psi |`$ and the phase $`\varphi `$ of the field $`\psi =|\psi |e^{i\varphi }`$. On the lattice, with lattice spacing $`\epsilon `$, we normalize the hamiltonian by setting $`R^2=|\psi |^2/(a/b)`$, $`\stackrel{}{u}=\stackrel{}{r}/\xi _0`$, where $`\xi _0^2=\gamma /a`$ is the mean field correlation length at zero temperature. The normalized hamiltonian is then:

$$H[R,\varphi ]=k_BV_0\left(H_R+\underset{i=1}{\overset{N}{}}R_i^2f_i\right)$$
(2)
where $`f_i:=_{j=1}^d\left[1\mathrm{cos}(\varphi _i\varphi _j)\right]`$ and

$$H_R:=\underset{i=1}{\overset{N}{}}\left[\sigma \left(tR_i^2+R_i^4/2\right)+\left(R_i\right)^2/2\right].$$

where $`R_i:=_j\stackrel{}{e}_j(R_iR_j)`$. $`j`$ points to the nearest neighbours of $`i`$ and $`\stackrel{}{e}_j`$ is a unit vector in direction $`j`$. We have set $`R_iR_j=R_i((R_jR_i)+R_i)R_i^2`$ in the XY part of the hamiltonian. Indeed, the term $`R_i(R_jR_i)`$ is not important for the present discussion because our approach will be reliable when amplitude fluctuations are small compared to phase fluctuations. Only two parameters, that are in competition, remain:

$$\sigma :=\epsilon ^2/\xi _0^2V_0:=\frac{1}{k_B}\frac{a}{b}\gamma \epsilon ^{d2}$$

$`\sigma `$ controls amplitude fluctuations. $`V_0`$ corresponds to the zero temperature phase stiffness. $`V_0`$ is proportional to the superfluid density $`a/b`$ and controls the general critical behaviour. When $`V_0/T_0`$ is large, the critical region is small and the material behaves according to mean field theory (as BCS superconductors). When $`V_0/T_0`$ is of the order of 1, phase fluctuations become very large and give an upper bound for the critical temperature . One has to bear in mind that $`V_0`$ is not a constant that only normalizes the temperature $`T`$ in the Boltzmann factor because of the occurrence of the temperature $`t=T/T_01`$ inside the hamiltonian $`H[\psi ]`$.
Let us write the partition function in polar coordinates:

$$Z=_0^{2\pi }D\varphi _0^{\mathrm{}}DRe^{\beta H_{eff}}$$

where the effective hamiltonian $`H_{eff}`$ is

$$H_{eff}=k_BV_0\left[H_R+\underset{i}{}\left(R_i^2f_i\frac{T}{V_0}\mathrm{log}R_i\right)\right]$$
(3)
We keep the factor $`R`$ of the Jacobian, $`RD\varphi DR`$, and take it in the exponential giving a contribution $`logR`$ to the potential. This incorporates the fact that small values of $`R`$ have a small statistical weight, due to the volume element in phase space, into the Boltzmann factor of the canonical ensemble. We compute now the partition function by integrating only the phase. For this purpose, we will drop the integration on $`R_i`$ and search for the minimum of the free energy with respect to $`R`$. The partition function becomes then:

$$Z_0^{2\pi }D\varphi e^{\beta H_{eff}}$$

The minimum of the free energy $`F=k_BT\mathrm{log}Z`$ is given through the equation $`\delta F/\delta R_i=0`$. Assuming that the gradient of the amplitude is zero, we have

$$\sigma \left(tR^2+R^4\right)\frac{T}{2V_0}+R^2f(K)=0$$
(4)
where we have multiplied the equation by $`R/2`$. $`f(K)=f_i`$ is the expectation value within the XY model of $`f_i`$ with a dimensionless coupling constant $`K=\frac{V_0}{T}R^2`$. Although $`K`$ has the same value at each lattice site and does not explicitely adapt to the vortex structure of the phase field, the average energy of the latter still determines the value of $`K`$ through the minimalisation of the free energy $`F`$. In this work, we take Monte-Carlo simulations to evaluate the function $`f(K)`$ which is just the energy of the XY model . The critical temperature $`T_c`$ is reached when the coupling $`K`$ equals the critical constant $`1/\omega `$:

$$1=\omega K_c=\omega \frac{V_0}{T_c}R^2(T_c)$$
(5)
where $`\omega =0,0.9,2.2`$ (for $`d=1,2,3`$ respectively) is the pure XY critical temperature on a square lattice with unit coupling constant $`K=1/T`$.

The solution of equation (4) is plotted in figure 1 for $`d=3`$. We use a normalization which is independent of the temperature $`t`$: with $`\stackrel{~}{R}^2=R^2/(t)`$, $`\stackrel{~}{T}=T/(tV_0)`$ and $`\stackrel{~}{\sigma }=t\sigma `$. The free energy $`F`$, calculated as: $`F=𝑑RF/R`$, has three branches. The crossing point of the two lower branches determines the location of the first order transition, marked by a vertical dashed line. The consequence is a first order transition for

σ~c=σtc<

 4.5,subscript~𝜎𝑐𝜎subscript𝑡𝑐<

4.5\tilde{\sigma}_{c}=-\sigma t_{c}{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 5.0pt\vbox{\hbox{$\sim$}}}}\ }4.5, (6)
i.e. the mean amplitude $`R^2`$ jumps at $`T_c`$. For $`d=2`$, we have also a first order transition for σ~c<

 1.25subscript~𝜎𝑐<

1.25\tilde{\sigma}_{c}{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 5.0pt\vbox{\hbox{$\sim$}}}}\ }1.25. For $`d=1`$, there is no transition.

A good approximation for the critical temperature due to phase fluctuations is given by (5) with $`R^2=t_c`$:

$$t_cT_c/(\omega V_0)$$
(7)
Combining this last equality and (6), the first order transition occurs then for σTc/V0<

C𝜎subscript𝑇𝑐subscript𝑉0<

𝐶\sigma{T_{c}/V_{0}}{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 5.0pt\vbox{\hbox{$\sim$}}}}\ }C where $`C=\omega 4.59.9`$ for $`d=3`$.

What are the influence of amplitude fluctuations around the saddle point that we found? Do they destroy the first order transition or not? To answer to this question, we include harmonic amplitude fluctuations in our computation. We separate phase and amplitude, introduce three variational parameters and derive two self consistent equations for these parameters.

The effective potential $`U_{eff}:=\sigma \left(tR_i^2+\frac{1}{2}R_i^4\right)\frac{T}{V_0}\mathrm{log}R_i`$ of the effective hamiltonian (3) has only one minimum for all temperatures. So one can expect that a gaussian approximation for the amplitude gives a good approximation for all temperatures. We want also to get an hamiltonian with no direct coupling between phase and amplitude, but with effective constants that keep the memory of their interaction. The idea is then to separate phase and amplitude as:

$$R_i^2f_iR_i^2f_i+R_i^2f_i$$
(8)
Therefore, we set the trial hamiltonian:

$$H_t[R,\varphi ]=\underset{i}{}[B(R_iR_0)^2+\frac{1}{2}\left(R_i\right)^2+Af_i]$$
(9)
Using the Bogoliubov inequality, we have $`FF_t+H_{eff}H_t_t=\stackrel{~}{F}`$, where $`F`$ is the free energy and $`\mathrm{}_t`$ is the canonical average with respect to $`H_t`$. The right hand side is thus to be minimized with respect to the constants $`A`$, $`R_0`$ and $`B`$ to give the best approximation of $`F`$. We introduce also the local amplitude fluctuation $`\eta _i=R_iR_0`$. The derivative of $`\stackrel{~}{F}`$ with respect to these parameters gives three equations:

$`A=R^2`$ (10)
$`\left[\sigma (t+3\eta ^2)+f\right]R_0^2+\sigma R_0^4{\displaystyle \frac{T}{V_0}}{\displaystyle \frac{R_0}{2}}{\displaystyle \frac{\mathrm{log}R}{R_0}}=0`$ (11)
$`\sigma (t+3\eta ^2)+fB+3\sigma R_0^2{\displaystyle \frac{T}{V_0}}{\displaystyle \frac{\mathrm{log}R}{B}}=0`$ (12)
where all indices $`t`$ and $`i`$ are dropped. $`\mathrm{log}R_i`$ is computed by a cumulant expansion. Equations (11) and (12) are the central result of this work. The mean square amplitude fluctuation is

$$\eta ^2=\frac{T}{V_0}1/V\underset{|\stackrel{}{k}|<\mathrm{\Lambda }}{}\frac{1}{2}\frac{1}{B+k^2/2}$$

where we extend the lower bound $`R_0`$ to $`\mathrm{}`$ assuming that amplitude fluctuations are small. $`\mathrm{\Lambda }`$ is the reduced cut-off parameter and is computed on the first Brillouin zone.
Equation (11) reduces to equation (4) when the additional amplitude fluctuations are zero: $`\eta ^2=0`$. The phase diagram is approximately the same as the one without amplitude integration, except that the critical temperature is smaller due to the additional harmonic amplitude fluctuations. The transition is first order for σ<

 9.9V0/Tc𝜎<

9.9subscript𝑉0subscript𝑇𝑐\sigma{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 5.0pt\vbox{\hbox{$\sim$}}}}\ }9.9\ V_{0}/T_{c} for $`d=3`$.

In figure 2, we compare the simulations of Nguyen and Sudbø  with our phase-amplitude separation method for a corresponding $`\sigma /(V_0/T_0)=18`$ (continuous regime). We see that the simulated $`R^2`$ shows an inflexion point at $`T_c`$ as the analytical $`R^2`$ does. This is perhaps the forerunner for a jump at a smaller $`\sigma `$. Both methods agree with the fact that there is no first order transition at large $`\sigma `$. To our knowledge, accurate simulations for smaller $`\sigma `$ are not avalable.

A remarkable thing is the fact that a negative $`\sigma `$ produces a first order transition ( p.340). Therefore, the critical line which is believed until now to separate the region $`\sigma <0`$ and $`\sigma >0`$ is pushed in the region with $`\sigma >0`$ by phase fluctuations. The limiting case $`\sigma <0`$ is then consistent with a first order transition for σ<

 9.9V0/Tc𝜎<

9.9subscript𝑉0subscript𝑇𝑐\sigma{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 5.0pt\vbox{\hbox{$\sim$}}}}\ }9.9\ V_{0}/T_{c}. We mention here that a variational approximation due to Halperin et al  for gauge field fluctuations yields an equation that is very similar to equation (4) except that our function $`f(K)=f_i`$ due to phase fluctuations is replaced by the expectation value $`\stackrel{}{A}^2`$. For $`d=3`$, magnetic fluctuations of the gauge field $`\stackrel{}{A}`$ produce a first oder transition and move the critical point into the positive $`\sigma `$ region. The domain of validity is also in the small $`\sigma `$ region as for our approach.

We now show that the origin of the first order transition of equations (11) and (12) is not due to the approximation itself and that the domain of validity of the approximation is precisely in the first order domain for $`d=3`$. We establish a quantitative criterion by comparing the size of the critical gaussian region according to the Ginzburg criterion with the critical temperature due to phase fluctuations. The Ginzburg criterion  measures the importance of correlated gaussian fluctuations that are mainly due to the amplitude, i.e. they are negligible if

$$|t|>>t_G:=\left(\frac{T_c}{2V_0}\right)^{\frac{2}{4d}}\sigma ^{\frac{d2}{4d}}$$

where $`t_G`$ defines the critical region. $`|t_c|`$ from equation (7 is a measure of the importance of phase fluctuations. Our approach is then valid when (amplitude) fluctuations are dominated by phase fluctuations, i.e. when $`|t_c|>>t_G`$. Setting $`\omega d1`$ which is almost the correct XY critical temperature, we get the criterion:

$$\frac{d1}{2}\left(\sigma \frac{T_c}{2V_0}\right)^{\frac{d2}{4d}}<<1$$
(13)
For $`d=1`$, the criterion is in agreement with our results (not shown here) which are almost identical to the exact solution . For $`d=2`$, the criterion is not clear. The approximation is reliable precisely in the first order domain for the 3$`d`$ case ($`\sigma T_c/V_0<<1`$), and everywhere for $`d=4`$. For $`d>4`$, the domain of validity is in the continuous domain. Therefore, we conclude that the first order transition is not a consequence of the approximation for $`d=3`$, whereas a doubt remains for $`d=2`$.

The first order transition is due to the fact that the energy of vortices can be lowered by a reduced value of the amplitude, which is energetically favorable when the potential energy of the amplitude is sufficiently soft. In this case, the usual transition scenario, given by an unbinding of topological excitations, seems to be replaced by a sudden proliferation of the latter. Such a process may be difficult to describe within the renormalisation group $`\epsilon `$-expansion where the interplay between amplitude and phase is not explicitely taken into account . A similar scenario has been found by Minnhagen et al.  for the 2$`d`$ Coulomb gas. They show that a large vortex fugacity $`y=\mathrm{exp}(\beta \gamma |\psi |^2\pi ^2/2)`$ can lead to a discontinuous transition produced by a proliferation of vortices, whereas small fugacity causes the usual Kosterlitz-Thouless transition. When $`\stackrel{~}{\sigma }`$ is small, our scenario corresponds to the case of large vortex fugacity. The mean amplitude can become much smaller than the one of the XY model (see figure 1 where $`\stackrel{~}{R}=1`$ corresponds to the XY model). For large $`\stackrel{~}{\sigma }`$ ($`\stackrel{~}{R}1`$) , the KT transition is recovered. Our result ($`\stackrel{~}{\sigma }_c1.25`$) is in good agreement with the one of Bormann and Beck  who found a first order transition for σ~<

 1~𝜎<

1\tilde{\sigma}{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 5.0pt\vbox{\hbox{$\sim$}}}}\ }1.

For real superconductors, $`\sigma `$ is smaller than 1 . Therefore, superconductors with low $`T_c/V_0`$ could have a first order transition. A candidate for a possible observation would be a 3$`d`$ superconductor with a critical region $`t_G`$ that is not too small but still in the first oder domain. In the BCS limit, the transition becomes very weakly first order such that there is no contradiction with the continuous behaviour of BCS superconductors. Underdoped cuprates are quasi 2 dimensional and have a low $`V_0/T_c`$. Their transition is then XY-like, whereas overdoped cuprates, that are almost 3$`d`$ and have a large $`V_0/T_c`$, could have an observable first order transition. It is, however, interesting to remark that measurements of the entropy change $`\mathrm{\Delta }S`$ at the vortex lattice melting transition of Bi2Sr2CaCu2O8  have shown a dramatic increase in $`\mathrm{\Delta }S`$ per vortex when the zero field transition is approached. This could be a hint that the superconducting transition remains first order even in zero field.

In this paper we have investigated the thermodynamic properties of the classical Ginzburg-Landau (GL) model. It is determined by two model parameters, $`\sigma `$ and $`V_0`$. $`\sigma `$ governs the strength of amplitude fluctuations and $`V_0`$ the overall strength of fluctuations of the complex GL field. We have treated the model by a variational approximation which takes into account the coupling between phase and amplitude through effective coupling constants. Minimizing the corresponding variational free energy leads to a set of self-consistent modified GL equations containing phase and amplitude fluctuations. The behaviour of the GL transition changes when the ratio $`\sigma /(V_0/T_c)`$ is varied: for σ<

CV0/Tc𝜎<

𝐶subscript𝑉0subscript𝑇𝑐\sigma{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 5.0pt\vbox{\hbox{$\sim$}}}}\ }C\ V_{0}/T_{c} the transition of first order with $`C=9.9`$ for $`d=3`$ and $`C=2.2`$ for $`d=2`$.

For $`d=3`$, we showed that phase fluctuations dominate the transition in the first order domain and amplitude fluctuations can be neglected. We thus conclude that a first order transition is indeed a valid scenario for the GL model, once the amplitude of the field can sufficiently adapt in order to lower the total energy of the system.

We thank D. Baeriswyl, X. Bagnoud, M. Dzierzawa, C. Verdon, H. Fort and C. Wirth for very useful discussions. This work has been supported by the Swiss National Science Foundation (project No. 2000-056803.99).


format: text

text:
# Resonance fluorescence in a band gap material: Direct numerical simulation of non-Markovian evolution

## I Introduction

Photonic band-gap (PBG) materials give rise to an altered density of states (DOS) for the electromagnetic field, which leads to a fundamentally different atom-radiation field interaction . When atoms emit light into these structures, they exhibit inhibition of spontaneous emission , localization of light, and the formation of an atom-photon bound state . These effects are due to entanglement between the atomic and photonic states, and as such represent a failure of the Born-Markov approximation, which allows the derivation of a master equation for the reduced density matrix of the atom and has been an extraordinarily useful tool in the field of theoretical quantum optics for treating spontaneous emission. The spontaneous emission of photons into a PBG has been described analytically, but describing such a system with the added nonlinearity of driving has remained a theoretical challenge . This paper outlines a numerical method that calculates the resonance fluorescence of a driven atom in a PBG material. It may be applicable to other weakly non-Markovian systems.

A photonic band gap is a range of frequencies in which no photon propagation can occur. An atom near one of these band edges will experience an interaction with an altered electromagnetic vacuum, and will therefore decay in a different manner. Trivially, there will be no emission of photons into the band gap, but the altered DOS near the band edge will have further effects on the behavior. The case of spontaneous emission was treated by John and Quang  who, due to the restricted state space in this case, were able to solve the equations of motion for the system analytically by Laplace transform methods. They found that the atom would not decay to the ground state due to the presence of a localized atom-photon bound state which was more strongly bound with stronger coupling. Near a band edge, the energy of the propagating photons varies quadratically with the wavevector in the same way as a massive particle propagating in free space. The model of this system is therefore virtually identical to that of a single mode atom laser .

The extension of the work of John and Quang to a driven atom has remained a theoretical challenge. When pumping was added to the atom laser model, it allowed access to a potentially infinite ladder of states within the lasing mode. Fortunately, the linearity of the problem allowed a closed set of equations for the relevant expectation values to be developed, and the steady state output spectrum of the system was calculated . In contrast, a driven two level atom in a PBG material has a much smaller system (two states), but the nonlinear behavior of the atom due to the driving, in conjuction with the output coupling to the field, does not allow a closed set of equations to be deduced. Quang and John attempted a treatment in which they assumed that no more than one photon was ever present in the output field , but in a subsequent comment this assumption was shown to be equivalent to the Born approximation, and the Monte Carlo wavefunction method employed to deal with this problem was shown to be invalid .

The memory of a system with a quadratic dispersion relation in one dimension is effectively infinite. This is evident in the undriven case where photons can be permanently localized near the atom, which can remain partially in the excited state forever. Driving the atom in the presence of this dressed atom-photon bound state will cause the bound state to become more highly occupied. A driven atom in a PBG material is analogous to the pumped atom laser model, which does not have a steady state without the inclusion of a loss mechanism for the bound state . In a real system the build-up of energy in the bound state would reach some limit as imperfections in both the PBG material and the quadratic approximation became significant. In other words, describing the resonance fluorescence of an atom in such a material requires inclusion of the losses in order to produce a physical result. This will result in a finite memory time for the evolution of the atomic system. In addition, the calculations so far have dealt with an isotropic model of a bandgap but more appropriate for artificially created bandgap materials is the anisotropic bandgap model  which leads to memory function that decays much more rapidly than that for an isotropic bandgap. In this paper we will present a method which is able to calculate the evolution of a driven atom near a PBG providing the memory time of the system is sufficiently short, and the coupling is sufficiently weak.

In the next section we will describe our model for the driven atom, and demonstrate how a formal solution for the reduced density matrix of the atom can be generated. In Sec. III we will describe a numerical method for producing results with this formal solution. Sec. IV shows how our method reproduces the exact result from a similar, soluble system, and then describes the spectrum of the resonance fluorescence near a band gap in a PBG material with varying levels of detuning of the atomic resonance from the bandgap edge.

## II Model

We consider a two level atom with atomic energy level separation $`\mathrm{}\omega _0`$. It is described by the Pauli spin matrices $`\sigma `$, $`\sigma ^{}`$ and $`\sigma _z`$, which satisfy the commutation relations $`[\sigma ^{},\sigma ]=\sigma _z`$, and $`\{\sigma ^{},\sigma \}=1`$. The atom is driven by a classical field of frequency $`\omega _L`$, and is modeled by the Hamiltonian

$$H_0=\frac{\mathrm{}\omega _0}{2}\sigma _z+\frac{\mathrm{}\mathrm{\Omega }}{2}\{\mathrm{exp}(i\omega _Lt)\sigma ^{}+\mathrm{exp}(i\omega _Lt)\sigma \},$$
(1)
where $`\mathrm{\Omega }`$ is the strength of the driving field experienced by the atom.

Under the rotating wave approximation  the Hamiltonian for the atom coupled to a large number of discrete modes of the electromagnetic field can be written as

$$H=H_0+\underset{k}{}\left[\mathrm{}\omega _kb_k^{}b_k+i\mathrm{}(g_k^{}b_k^{}\sigma g_kb_k\sigma ^{})\right]$$
(2)
where $`b_k`$ is the annihilation operator of the $`k`$-th mode, characterized by the frequency $`\omega _k`$, and is coupled to the atom with the coupling strength $`g_k`$. It is important to note that the coupling is linear in the field operators.

For the formal manipulations which we will perform below it is convenient to work in the interaction picture. In the interaction picture the Hamiltonian can be written as

$$H_I(t)=i\mathrm{}[\xi ^{}(t)\sigma (t)\xi (t)\sigma ^{}(t)],$$
(3)
where the driving field, $`\xi (t)`$, is defined by

$$\xi (t)=\underset{k}{}g_kb_ke^{i\omega _kt},$$
(4)
and the interaction picture atomic operator $`\sigma (t)`$ is given by $`\sigma (t)=U_0^{}(t)\sigma U_0(t)`$. The free evolution operator $`U_0`$ in this equation is defined by $`U_0(t)=e^{\frac{i}{\mathrm{}}H_0t}`$.

A quantity fundamental to any attempt to describe this interaction in terms of the atomic system alone is the memory function, defined as the commutator of the driving field with its complex conjugate;

$$f_\mathrm{m}(t,t^{})=[\xi (t),\xi ^{}(t^{})].$$
(5)
For $`t>t^{}`$, the memory function is the probability amplitude of absorbing a photon at time $`t`$ that was emitted previously at time $`t^{}`$. It is therefore the Feynman propagator of a virtual photon. Non-Markovian behavior arises when one cannot consider the lifetime of the virtual photons to be vanishingly small. The Born-Markov approximation can be made for free space resonance fluorescence, as (for sufficiently high frequencies) the linear dispersion relation ($`\omega _k=ck`$) produces an infinitely narrow memory function. Unfortunately, near a photonic band gap the dispersion relations can be quadratic, and this leads to a very long memory time.

The combination of non-Markovian evolution and the nonlinear behavior of a driven atom renders this problem intractable to either Born-Markov master equation methods applicable to free-space resonance fluorescence , or Laplace transform methods that have been applied to the case of spontaneous emission into a band gap material . Instead, we approach this problem by a direct numerical solution of the evolution of the reduced atomic system.

The first important step is the derivation of a formal equation for the reduced density matrix of the atomic system. The reduced density matrix is defined by tracing over the field degrees of freedom of the density matrix for the full interacting system. In the interaction picture it has the form

$$\rho (t)=\mathrm{Tr}_{\mathrm{field}}\left\{\stackrel{}{T}\left[e^{\frac{i}{\mathrm{}}_0^t𝑑sH_I(s)}\right]\rho _{\mathrm{tot}}\stackrel{}{T}\left[e^{\frac{i}{\mathrm{}}_0^t𝑑uH_I(u)}\right]\right\},$$
(6)
where $`(\stackrel{}{T})\stackrel{}{T}`$ denotes (anti-)timeordering. Assuming an initial factorized state of the atom and field, $`\rho _{\mathrm{tot}}=\rho \rho _{\mathrm{field}}`$, and that the field is initially in a generalized Gaussian state  (e.g. a thermal state), a formal expression in terms of only system variables can be determined for the right-hand side of Eq.(6). Here we will perform a derivation for the case of an initial vacuum field, $`\rho _{\mathrm{field}}=|\{0\}\{0\}|`$.

We use the identity

$$\stackrel{}{T}\left[e^{\frac{i}{\mathrm{}}_0^t𝑑sH_I(s)}\right]=\stackrel{}{T}_\sigma \left\{V_+V_0V_{}\right\}$$
(7)
where $`\stackrel{}{T}_\sigma `$ is a time ordering operator which denotes a time ordering on the atomic operators only, and where we have defined the operators

$`V_+`$ $`=`$ $`\mathrm{exp}\left({\displaystyle _0^t}𝑑s\xi ^{}(s)\sigma (s)\right)`$ (8)
$`V_0`$ $`=`$ $`\mathrm{exp}\left({\displaystyle _0^t}𝑑u{\displaystyle _0^u}𝑑vf_\mathrm{m}(u,v)\sigma ^{}(u)\sigma (v)\right)`$ (9)
$`V_{}`$ $`=`$ $`\mathrm{exp}\left({\displaystyle _0^t}𝑑s\xi (s)\sigma ^{}(s)\right)`$ (10)
where $`f_\mathrm{m}(u,v)`$ is the memory function, defined by Eq. (5). Similarly, we can write

$$\stackrel{}{T}\left[e^{\frac{i}{\mathrm{}}_0^t𝑑sH_I(s)}\right]=\stackrel{}{T}_\sigma \left\{V_{}^{}V_0^{}V_+^{}\right\}$$
(11)
for the anti-timeordered term.

In order to greatly simplify the manipulation of these time-ordered quantities we introduce a notational apparatus, first used by Feynman , whereby unprimed atomic operators are assumed to be time-ordered to the left of $`\rho `$ and primed operators anti-timeordered to the right. We can then write Eq. (6) as

$$\rho (t)=\mathrm{Tr}_{\mathrm{field}}\{V_+V_0V_{}(\rho |\{0\}\{0\}|)V_{}{}_{}{}^{}{}_{}{}^{}V_{0}^{}{}_{}{}^{}{}_{}{}^{}V_{+}^{}{}_{}{}^{}{}_{}{}^{}\},$$
(12)
where $`V_j{}_{}{}^{}=V_j(\xi ,\xi ^{},\sigma {}_{}{}^{},\sigma ^{}).`$ This notation allows us to manipulate the atomic operators as if they were $`c`$-numbers while keeping track of their correct ordering in the expression. We can use the cyclic properties of the trace over the field variables to write Eq.(12) in the form

$$\rho (t)=\{0\}|V_{}{}_{}{}^{}{}_{}{}^{}V_{0}^{}{}_{}{}^{}{}_{}{}^{}V_{+}^{}{}_{}{}^{}{}_{}{}^{}V_{+}^{}V_0V_{}|\{0\}\rho .$$
(13)
We now normally order the field operators in this expression by making use of the relations $`e^{A+B}=e^Ae^Be^{\frac{1}{2}[A,B]}=e^Be^Ae^{\frac{1}{2}[A,B]}`$ for two operators $`A`$ and $`B`$ satisfying $`[A,[A,B]]=[B,[A,B]]=0`$. The normally-ordered field operators are annihilated on the field vacuum and we are left with the result

$$\rho (t)=e^{(t)}\rho ,$$
(14)
where

$$(t)=_0^t𝑑u_0^u𝑑v\left\{[\sigma _{}^{}{}_{}{}^{}(u)\sigma ^{}(u)]f_\mathrm{m}(u,v)\sigma (v)+\left[\sigma (u)\sigma ^{}(u)\right]f_\mathrm{m}^{}(u,v)\sigma _{}^{}{}_{}{}^{}(v)\right\}.$$
(15)
This is a formal solution that contains no reference to the field operators. The only reference to the particular characteristics of the field is via the memory function. The $`\sigma ^{}(u)\sigma (v)`$ and $`\sigma ^{}(u)\sigma _{}^{}{}_{}{}^{}(v)`$ terms arise from the vacuum field to vacuum field evolution while the terms linking the two sides of the density matrix, $`\sigma ^{}(u)\sigma (v)`$ and $`\sigma (u)\sigma _{}^{}{}_{}{}^{}(v),`$ arise from tracing over all the photons that are irretrievably emitted into the field. Note that only under the Born-Markov approximation \[where it is possible to replace $`f_\mathrm{m}(t,s)`$ by $`\gamma \delta (ts)`$\] does the evolution matrix, $`e^{(t)},`$ factorize at each time, as in this case the double integral reduces to a single integral. It is then possible to write down a ordinary differential equation for $`\rho (t)`$. In general, however, this is not possible and it is necessary to use the formal solution Eq.(14).

The theoretical difficulty with using this formal solution is not the size of the system, which can be described by three real numbers, but with the implicit time ordering of the system operators, which in general leads to an infinite number of ordered terms. Physically, these higher order terms are the processes associated with increasing numbers of photons, and will become increasingly less significant in a system with a finite memory time. A driven system without a finite memory time may lead to unphysical behavior , so we are particularly interested in systems with a finite memory time. When this timescale is shorter than other dynamic timescales we are describing a Markovian system, and otherwise we are describing a non-Markovian system.

In this work we consider a numerical approach to solving the evolution matrix $`e^{(t)}`$. We discretize time and then make the approximation of a finite memory time $`T_\mathrm{m}`$. In other words, we assume that we can effectively put $`f_\mathrm{m}(t,s)=0`$ for all $`s<tT_\mathrm{m}`$. Our algorithm considers all possible processes that can occur with these restrictions. From Eq.(14) it is clear that it is not possible to determine $`\rho (t+\mathrm{\Delta }t)`$ from only $`\rho (t)`$ as at time $`t`$ there are photons in the field, represented by the presence of the memory function that will re-interact with the atom at time $`t+\mathrm{\Delta }t`$. These photons have to be taken into account in any algorithm. It turns out that under the finite memory time approximation we can follow the evolution of these photons by introducing a finite number of additional, “virtual”, density matrices, $`\rho ^k,`$ $`k0`$. The evolution equation for the reduced density matrix together with the virtual density matrices has the form

$$\rho ^k(t+\mathrm{\Delta }t)=\underset{j}{}D_{kj}(t)\rho ^j(t)$$
(16)
where $`D_{kj}(t)=D_{kj}[\sigma (t),\sigma ^{}(t),\sigma ^{}(t),\sigma _{}^{}{}_{}{}^{}(t)]`$ is an evolution superoperator and $`\rho ^0(t+\mathrm{\Delta }t)\rho (t+\mathrm{\Delta }t),`$ approximates the reduced density matrix.

The method is based on transforming the numerical algorithm for propagating a non-Markovian differential equation into a Markovian form. To motivate this idea, consider the following elementary example of an integro-differential equation,

$`{\displaystyle \frac{dc(t)}{dt}}=i\omega _0c(t){\displaystyle _{tT_\mathrm{m}}^t}f_\mathrm{m}(ts)c(s).`$

with a finite lower-bound on the integral. This integro-differential equation can be propagated numerically via the algorithm

$$c_{n+1}=(1i\omega _0\mathrm{\Delta }t)c_n\frac{\mathrm{\Delta }t^2}{2}\left[f_0c_n+2\underset{j=nM+2}{\overset{n}{}}f_{nj+1}c_j+f_Mc_{nM+1}\right]+O(\mathrm{\Delta }t^2)$$
(17)
where $`c_n=c(n\mathrm{\Delta }t),f_n=f_\mathrm{m}(n\mathrm{\Delta }t)`$ and $`M+1=T_\mathrm{m}/\mathrm{\Delta }t.`$ We have used the trapezoidal rule for the integral and Euler’s method to propagate the differential equation forward. Due to the finite memory time, this algorithm can be written in an explicitly Markovian form by introducing $`M1`$ extra variables $`b^k`$. The equivalent Markovian algorithm is

$`c_{n+1}`$ $`=`$ $`(1i\omega \mathrm{\Delta }t)c_n{\displaystyle \frac{\mathrm{\Delta }t^2}{2}}f_0c_n+b_n^1+O(\mathrm{\Delta }t^2)`$
$`b_{n+1}^1`$ $`=`$ $`b_n^2\mathrm{\Delta }t^2f_2c_n`$
$`b_{n+1}^2`$ $`=`$ $`b_n^3\mathrm{\Delta }t^2f_3c_n`$
$`\mathrm{}`$
$`b_{n+1}^{M1}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }t^2}{2}}f_Mc_n`$
with the initial condition $`b_0^k=0.`$ In this particular case, this transformation does not offer any advantage, however, this example is analogous to the algorithm for propagating the reduced density matrix that we will present in the next section. In that case, due to the presence of time-ordering, reformulation as a Markovian algorithm becomes invaluable.

## III Numerical Method

As a first step in determining an algorithm let us evaluate the integrals in the exponential in Eq.(14) on a grid,

$$(t)=\underset{j=0}{\overset{j=N}{}}\mathrm{\Delta }tL_j+O(\mathrm{\Delta }t^2)$$
(18)
where

$$L_j=\underset{k=\mathrm{max}(jM+1,0)}{\overset{j}{}}\mathrm{\Delta }t\left\{(\sigma _{}^{}{}_{j}{}^{}\sigma _j^{})f_{jk}\sigma _k+(\sigma _j\sigma _j^{})f_{jk}^{}\sigma _{}^{}{}_{k}{}^{}\right\}.$$
(19)
For convenience of notation we have written $`t_j=j\mathrm{\Delta }t`$ where $`t_N=t`$, $`a_j=a(t_j)`$ and $`f_n=W_nf(n\mathrm{\Delta }t).`$ The weights are given by

$$W_n=\{\begin{array}{cc}\frac{1}{2},& n=0\\ \frac{1}{2},& n=M1\\ 1,& \mathrm{otherwise}\end{array},$$
(20)
as we are using the trapezoidal rule (as the time ordering produces an integrand that is effectively discontinuous at each time step there is good reason to believe that a higher order rule will not improve accuracy). Note that we have also explicitly truncated the sums after $`MT_\mathrm{m}/\mathrm{\Delta }t`$ as the memory function is assumed to be zero for all additional terms. We can now write

$`e^{(t)}`$ $`=`$ $`{\displaystyle \underset{j}{}}e^{\mathrm{\Delta }tL_j}+O(\mathrm{\Delta }t^2),`$ (21)
$`=`$ $`{\displaystyle \underset{j}{}}(1+\mathrm{\Delta }tL_j)+O(\mathrm{\Delta }t^2),`$ (22)
where in the second line we have expanded out the exponentials at each time, $`t_j`$, to first order in $`\mathrm{\Delta }t`$. Note that it is still necessary to timeorder the operators in this expression. Let us attempt to write this expression in terms of the above mentioned density matrices $`\rho ^k`$ at one particular time $`n1`$, where we choose $`NM>n>M`$ as this is the most common case. First split the evolution matrix at $`n+M`$ and $`n`$,

$$e^{(t)}=\underset{k=n+M}{\overset{N}{}}(\mathrm{})\underset{j=n}{\overset{n+M1}{}}(1+\mathrm{\Delta }tL_j)\overline{\rho }+O(\mathrm{\Delta }t^2),$$
(23)
where it is convenient to separate out the factor $`\overline{\rho }=_{j=0}^{n1}(1+\mathrm{\Delta }tL_j)\rho (0)`$ as it does not take part in our subsequent manipulations. Then split the sums over $`k`$ in $`L_j`$ for $`j<n+M`$ into parts containing interaction picture operators before and after $`n`$,

$$(1+\mathrm{\Delta }tL_j)=x_{j,n}+\mathrm{\Delta }t(\sigma _{}^{}{}_{j}{}^{}\sigma _j^{})y_{j,n}+\mathrm{\Delta }t(\sigma _j\sigma _j^{})z_{j,n},$$
(24)
where

$`x_{j,n}`$ $`=`$ $`1+\mathrm{\Delta }t{\displaystyle \underset{k=n}{\overset{j}{}}}\mathrm{\Delta }t\left[(\sigma _{}^{}{}_{j}{}^{}\sigma _j^{})f_{jk}\sigma _k+(\sigma _j\sigma _j^{})f_{jk}^{}\sigma _{}^{}{}_{k}{}^{}\right]`$ (25)
$`y_{j,n}`$ $`=`$ $`{\displaystyle \underset{k=jM+1}{\overset{n1}{}}}\mathrm{\Delta }tf_{jk}\sigma _k`$ (26)
$`z_{j,n}`$ $`=`$ $`{\displaystyle \underset{k=jM+1}{\overset{n1}{}}}\mathrm{\Delta }tf_{jk}^{}\sigma _{}^{}{}_{k}{}^{}.`$ (27)
Consider the $`M`$ factors from $`j=n`$ to $`j=n+M1`$. Due to the truncation of the sums at $`M`$ terms, it is only these factors that contain interaction picture operators at time $`n`$. We can expand this product of $`M`$ factors, each of which has three terms, as a sum over $`3^M`$ terms

$$\underset{j=n}{\overset{n+M1}{}}(1+\mathrm{\Delta }tL_j)=\underset{m=0}{\overset{3^{M1}1}{}}g_{m,n}h_{m,n},$$
(28)
where $`g_{m,n}`$ contains only operators after $`n1`$ and $`h_{m,n}`$ contains only operators before $`n`$. Note that since $`y_{n+M1,n}=z_{n+M1,n}=0`$ we only need to consider a sum over $`3^{M1}`$ terms. To write down the explicit form for $`g_{m,n}`$ and $`h_{m,n}`$ it is convenient to introduce a new notation to label the terms. We write the index $`m`$ in trinary (i.e. base three), where the $`k^{\mathrm{th}}`$ trit (bit in base three) tells us which of the three factors to choose from the $`(1+\mathrm{\Delta }tL_k)`$ term (shown in Eq. (24). Our new index then lists the trits which correspond to the first term, $`\{A_1,\mathrm{},A_p\}`$; the second term, $`\{B_1,\mathrm{},B_q\}`$; or the third term, $`\{C_1,\mathrm{},C_{Mqp}\}`$. These lists replace the $`m`$ index. In this notation we can write

$`g_n^{\{A_1,\mathrm{},A_p;B_1,\mathrm{},B_q;C_1,\mathrm{},C_{Mqp}\}}`$ $`=`$ $`x_{n+M1,n}{\displaystyle \underset{j=1}{\overset{p}{}}}x_{n+A_j1,n}{\displaystyle \underset{k=1}{\overset{q}{}}}\mathrm{\Delta }t(\sigma _{}^{}{}_{n+B_k1}{}^{}\sigma _{n+B_k1}^{})`$ (30)
$`\times {\displaystyle \underset{l=1}{\overset{Mqp}{}}}\mathrm{\Delta }t(\sigma _{n+C_l1}\sigma _{n+C_l1}^{})`$
and

$$h_n^{\{B_1,\mathrm{},B_q;C_1,\mathrm{},C_{Mqp}\}}=\underset{k=1}{\overset{q}{}}y_{n+B_k1,n}\underset{l=1}{\overset{Mqp}{}}z_{n+C_l1,n}.$$
(31)
where, if the lists $`\{B_1,\mathrm{},B_q;C_1,\mathrm{},C_{Mqp}\}`$ are the empty set, $`\mathrm{}`$, by construction we have $`h_n^{\mathrm{}}=1`$. This allows us to define a set of $`3^{M1}`$ density matrices at time $`n1`$ of the form

$$\rho _{n1}^{\{B_1,\mathrm{},B_q;C_1,\mathrm{},C_{Mqp}\}}=h_n^{\{B_1,\mathrm{},B_q;C_1,\mathrm{},C_{Mqp}\}}\overline{\rho },$$
(32)
Note that $`\rho _{n1}^{\mathrm{}}`$ is defined such that it corresponds to $`\rho (t_{n1})`$ up to order $`\mathrm{\Delta }t^2`$.

If we can determine the transformation from $`\rho _{n1}^k`$ to $`\rho _n^j`$ we will have achieved our aim of writing the non-Markovian evolution in terms of a Markovian evolution in an extended space. This is done by writing $`\rho _n^j`$ in terms of $`\rho _{n1}^k`$. Noting that $`y_{j,n+1}=y_{j,n}+\mathrm{\Delta }tf_{jn}\sigma _n`$ and $`z_{j,n+1}=z_{j,n}+\mathrm{\Delta }tf_{jn}^{}\sigma _{}^{}{}_{n}{}^{}`$ we can write

$`\rho _n^{\{B_1,\mathrm{}B_q;C_1,\mathrm{}C_{Mqp}\}}`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{q}{}}}[y_{n+B_k,n}+\mathrm{\Delta }tf_{B_k}\sigma _n]{\displaystyle \underset{j=1}{\overset{Mqp}{}}}[z_{n+C_j,n}+\mathrm{\Delta }tf_{C_j}^{}\sigma _{}^{}{}_{n}{}^{}]`$ (34)
$`\times [x_{n,n}+\mathrm{\Delta }t(\sigma _{}^{}{}_{n}{}^{}\sigma _n^{})y_{n,n}+\mathrm{\Delta }t(\sigma _n\sigma _n^{})z_{n,n}]\overline{\rho }.`$
Expanding out the first two sets of factors and keeping only first-order terms the right hand side can now be written in terms of the density matrices at $`n1`$,

$`\rho _n^{\{B_1,\mathrm{},B_q;C_1,\mathrm{},C_{Mqp}\}}`$ $`=`$ $`\left[1+\mathrm{\Delta }t^2(\sigma _{n}^{}{}_{}{}^{}\sigma _n^{})f_0\sigma _n+\mathrm{\Delta }t^2(\sigma _n\sigma _n^{})f_0^{}\sigma _{n}^{}{}_{}{}^{}\right]\rho _{n1}^{\{B_1+1,\mathrm{},B_q+1;C_1+1,\mathrm{},C_{Mqp}+1\}}`$ (39)
$`+\mathrm{\Delta }t\sigma _n{\displaystyle \underset{k=1}{\overset{q}{}}}f_{B_k}\rho _{n1}^{\{B_1+1,\mathrm{},B_{k1}+1,B_{k+1}+1,\mathrm{},B_q+1;C_1+1,\mathrm{},C_{Mqp}+1\}}`$
$`+\mathrm{\Delta }t\sigma _{}^{}{}_{n}{}^{}{\displaystyle \underset{k=1}{\overset{Mqp}{}}}f_{C_k}^{}\rho _{n1}^{\{B_1+1,\mathrm{},B_q+1;C_1+1,\mathrm{},C_{k1}+1,C_{k+1}+1,\mathrm{},C_{Mqp}+1\}}`$
$`+\mathrm{\Delta }t(\sigma _{}^{}{}_{n}{}^{}\sigma _n^{})\rho _{n1}^{\{1,B_1+1,\mathrm{},B_q+1;C_1+1,\mathrm{},C_{Mqp}+1\}}`$
$`+\mathrm{\Delta }t(\sigma _n\sigma _n^{})\rho _{n1}^{\{B_1+1,\mathrm{},B_q+1;1,C_1+1,\mathrm{},C_{Mqp}+1\}}+O(\mathrm{\Delta }t^2).`$
Physically, during a short enough interval only one photon can be emitted or absorbed at a time. The higher order terms in $`\mathrm{\Delta }t`$ correspond to two or more particle emissions or absorptions during this interval and hence we have dropped them. That is, except for the terms $`\mathrm{\Delta }t^2(\sigma _{n}^{}{}_{}{}^{}\sigma _n^{})f_0\sigma _n+\mathrm{\Delta }t^2(\sigma _n\sigma _n^{})f_0^{}\sigma _{n}^{}{}_{}{}^{},`$ which we include because often the memory function diverges at the origin and $`f_01/\mathrm{\Delta }t`$, so this term is effectively of order $`\mathrm{\Delta }t.`$ The additional density matrices are virtual in the sense that they represent states that are never observed in the system behavior and also in the related sense that they follow the paths of virtual photons. This equation has the form of Eq.(16) and therefore the right hand side can be put in the form of a superoperator acting on a density matrix. In this way it is possible to propagate the non-Markovian evolution forward in time.

Any expectation values are determined from $`\rho _n^{\mathrm{}}`$ alone;

$$a_n=\mathrm{Tr}\{a_n\rho _n^{\mathrm{}}\}.$$
(40)
For this reason it is convenient to introduce the projection operator $`P`$ defined by $`\rho _n^{\mathrm{}}=P\underset{¯}{\rho }_n,`$ where $`\rho _n^{\mathrm{}}`$ is projected out of the “vector” of matrices $`\underset{¯}{\rho }_n`$ which includes the reduced density matrix and all of the virtual density matrices. In this notation we can write the expectation value Eq.(40) as

$$a_n=\mathrm{Tr}\{a_nP\underset{¯}{\rho }_n\}$$
(41)
and the evolution equation (39) as

$$\underset{¯}{\rho }_n=D_n\underset{¯}{\rho }_{n1},$$
(42)
where $`D_n=D(\sigma _n,\sigma _n^{},\sigma _n^{},\sigma _{n}^{}{}_{}{}^{}).`$ Often one is not only interested in expectation values at a single time but also in temporal correlations. For example, the correlation function (with $`N_2>N_1`$)

$$\sigma _{N_2}^{}\sigma _{N_1}=\mathrm{Tr}\left\{\sigma _{N_2}^{}\varrho _{N_2,N_1}\right\}$$
(43)
where, assuming that the initial state of the field is the vacuum state, we have written

$$\varrho _{N_2,N_1}=\sigma _{N_1}e^{(N_2\mathrm{\Delta }t)}\rho .$$
(44)
In general, techniques for calculating the correlation function such as the quantum regression theorem (see, for example, ) will not work in the non-Markovian case as there is no factorization of the evolution matrix $`e^{(N_2\mathrm{\Delta }t)}`$ at time $`N_1\mathrm{\Delta }t`$. However, it is not difficult to see that this type of correlation function is easily calculated numerically within the above formalism. One simply propagates the set of density matrices $`\underset{¯}{\rho }_n`$ forward from $`n=0`$ as usual; then at $`n=N_1`$ multiply by $`\sigma _{N_1}`$; taking this as the new initial state one can propagate up to time $`N_2`$ giving $`\varrho _{N_2,N_1}`$ to the accuracy of the algorithm. In symbolic notation this procedure is written as

$$\varrho _{N_2,N_1}P\stackrel{}{}_{k=N_1+1}^{N_2}D_k\sigma _{N_1}\stackrel{}{}_{j=1}^{N_1}D_j\underset{¯}{\rho }_0,$$
(45)
where $`\stackrel{}{}`$ denotes a time-ordered product.

For numerical simulations it is often more useful to work in the Schrödinger picture. Equation (39) can easily be put in the Schrödinger picture because the operators are all acting at the same time. We simply multiply by $`U_0(n\mathrm{\Delta }t)U_{0}^{}{}_{}{}^{}(n\mathrm{\Delta }t)`$ which transforms the density matrices at time $`n`$ to the Schrödinger picture and then write the interaction picture operators on the right hand side in terms of the Schrödinger picture operators. The resulting equation has the same form as Eq.(39) (with the density matrices in the Schrödinger picture) except the right hand side is multiplied by $`U_0(\mathrm{\Delta }t)U_{0}^{}{}_{}{}^{}(\mathrm{\Delta }t)`$;

$$\underset{¯}{\rho }_n^S=D(\sigma ,\sigma ^{},\sigma ^{},\sigma _{}^{}{}_{}{}^{})U_0(\mathrm{\Delta }t)U_{0}^{}{}_{}{}^{}(\mathrm{\Delta }t)\underset{¯}{\rho }_{n1}^S.$$
(46)
The advantage of the Schrödinger picture is that the superoperator $`D`$ does not change with time and so only needs to be constructed once at the start of the simulation.

It is worth noting that the above superoperator representation does not need to be abandoned in numerical simulations. We can, for example, flatten the density matrix into a column vector so that both the primed and unprimed operators (which consequently become $`4\times 4`$ matrices rather than $`2\times 2`$ matrices) act only on the left-hand side of this flattened density matrix . The right-hand side of Eq.(46) then represents the $`4(3^M)\times 4(3^M)`$ sparse matrix $`D(\sigma ,\sigma ^{},\sigma ^{},\sigma _{}^{}{}_{}{}^{})`$ $`U_0(\mathrm{\Delta }t)U_{0}^{}{}_{}{}^{}(\mathrm{\Delta }t)`$ multiplying the $`(3^M)4\times 1`$ column vector $`\underset{¯}{\rho }_n^S`$. Note that the initial vector $`\underset{¯}{\rho }_0`$ has all elements zero except for $`\rho _0^{\mathrm{}}=\rho `$. Also note that we have written $`M`$ not $`M1`$, as some virtual density matrices are initialized at each time step.

## IV Results

Before we present the results of numerical simulations using the above algorithm, it is important to make some general comments about the limitations of the algorithm. As mentioned above, to propagate the non-Markovian evolution requires a vector of size $`(3^M)4\times 1`$ to be updated at each time step, thus computer memory considerations put a severe limitation on the number of time steps per memory time, $`M`$, that the algorithm can deal with. In our case, we created the time independent evolution matrix at the beginning and were in fact limited by the number of non-zero entries of this matrix that could be stored in computer memory. On a computer with 128MB of RAM we were limited to a maximum of $`M=11`$. This obviously poses quite a problem when one wishes to consider quite non-Markovian situations but also require high accuracy. However, for test cases where the system was not strongly non-Markovian, the results of our simulations show good agreement with calculations using standard methods of analysis.

In the rest of this article we consider the case of on-resonance driving ($`\omega _0=\omega _L`$). In a frame rotating at the atomic frequency the Hamiltonian for the atom reduces to

$`H_0={\displaystyle \frac{\mathrm{}\mathrm{\Omega }}{2}}\left(\sigma ^{}+\sigma \right).`$

This Hamiltonian gives rise to the following behavior

$`\sigma (t)`$ $`=`$ $`e^{\frac{i}{\mathrm{}}H_0t}\sigma e^{\frac{i}{\mathrm{}}H_0t}`$
$`=`$ $`{\displaystyle \frac{1}{2}}(\sigma _x+|+|e^{i\mathrm{\Omega }t}|+|e^{i\mathrm{\Omega }t})`$
where $`\sigma _x=\sigma ^{}+\sigma `$ and the states $`|\pm `$ are defined by $`\sigma _x|\pm =\pm |\pm .`$ When the atom is coupled to a propagating field these three terms give rise to three distinct peaks in the spectrum at $`\omega \omega _0=0,`$ $`\pm \mathrm{\Omega }`$. This is referred to as the resonance fluorescence triplet .

With a structured electromagnetic field, these emission peaks may be altered. In the following sections we use our algorithm to consider the case when the atom is in a cavity tuned to one of the side peaks and the case when one of the side-peaks is at a frequency where no photons propagate in a bandgap material.

### A Atom in a Cavity

In order to demonstrate the validity of our algorithm we consider a simple system that can also be described by a Born-Markov master equation (and so can be solved by traditional methods) and compare the results of the two methods. The system that we consider is an atom coupled to a single mode of a cavity which is in turn coupled to a continuum of propagating modes of the electromagnetic field via a partially transmissive mirror. The standard method of finding the evolution of this system is to treat the atom and the cavity as a system and the external field as a reservoir which can be traced over. This produces a Born-Markov master equation for the atom and cavity modes. This can then be solved by standard methods . One can also make a different choice of system and reservoir; consider the evolution of the atom alone while treating both the cavity and the external field as the reservoir; giving rise to a non-Markovian problem.

The usual description of this system is in terms of the Hamiltonian

$$H_1=H_0+\stackrel{~}{H}+i\mathrm{}\sqrt{\frac{\gamma }{2\pi }}\left[a^{}\sigma \sigma ^{}a\right],$$
(47)
where

$$\stackrel{~}{H}=\mathrm{}\nu a^{}a+\mathrm{}_{\mathrm{}}^{\mathrm{}}𝑑\omega \left\{\omega b^{}(\omega )b(\omega )+i\frac{\kappa }{\sqrt{2\pi }}\left[b^{}(\omega )ab(\omega )a^{}\right]\right\},$$
(48)
and where $`a`$ is the annihilation operator of the cavity mode, $`b(\omega )`$ is the annihilation operator for the external propagating electromagnetic field mode of energy $`\mathrm{}\omega `$, $`\sigma `$ is the lowering operator for the atom, and $`H_0`$ acts only on the atomic operators.

Choosing the atom as the system and eliminating the cavity mode it is possible to describe this system by the alternative Hamiltonian (see Appendix. A )

$$H_2=H_0+\mathrm{}_{\mathrm{}}^{\mathrm{}}𝑑\omega \left\{\omega c^{}(\omega )c(\omega )+i[\lambda ^{}(\omega )c^{}(\omega )\sigma \lambda (\omega )c(\omega )\sigma ^{}]\right\},$$
(49)
where $`c(\omega )`$ is a new bose annihilation operator and the, now frequency dependent, coupling constant is

$$\lambda (\omega )=\sqrt{\frac{\gamma }{2\pi }}\frac{\kappa }{i(\omega \nu )\frac{\kappa ^2}{2}}.$$
(50)
The memory function in this case is therefore given by

$`f_\mathrm{m}(ts)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega ^{}\lambda ^{}(\omega ^{})\lambda (\omega )[c(\omega ),c^{}(\omega ^{})]e^{i\omega t+\omega ^{}s}`$ (51)
$`=`$ $`{\displaystyle \frac{\gamma }{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega {\displaystyle \frac{\kappa ^2}{(\omega \nu )^2+\left(\frac{\kappa ^2}{2}\right)^2}}e^{i\omega (ts)}`$ (52)
$`=`$ $`\gamma \mathrm{exp}\left\{i\nu (ts){\displaystyle \frac{\kappa ^2}{2}}|ts|\right\}.`$ (53)
Transforming to a frame rotating at the atomic frequency appends a factor $`e^{i\omega _0(ts)}`$ to this memory function. This model differs from that of the bandgap only by the functional form of the memory function.

We compare our results for the system described by $`H_2`$ using the non-Markovian algorithm introduced in Sec. III with the solution of the Born-Markov master equation for the atom plus cavity system as described by $`H_1`$. In Fig.1 we plot the time evolution of the probability for the atom to be in the excited state. The internal spectrum, $`S(\omega )`$, defined as the Fourier transform of the steady-state correlation function $`\sigma ^{}(\tau )\sigma _{\mathrm{ss}}\sigma ^{}_{\mathrm{ss}}\sigma _{\mathrm{ss}}`$ , is plotted in Fig.2 showing a significant asymmetry in the height of the side peaks due to the tuning of the cavity to only one side band of the resonance fluorescence triplet. This is a completely non-Markovian effect as it arrises directly from the structured nature of the field, i.e., the fact that some frequencies of the field are easier to emit into than others. Both the probability and the spectrum show good agreement with the results of the traditional method (we used the numerical routines in the quantum optics toolbox ).

### B A Bandgap Material

In this section we apply our algorithm to resonance fluorescence in a bandgap material; a problem that cannot be solved by any traditional methods. Following the discussion of Vats and John  we consider an anisotropic bandgap model under the effective mass approximation. This approximation is appropriate to fabricated band gap materials. In this model the band edge is associated with a specific point in $`k`$-space, $`𝐤=𝐤_0`$, and the dispersion relation close to the upper edge of the gap has the form

$$\omega _𝐤=\omega _g+\frac{\mathrm{}\left(𝐤𝐤_0\right)^2}{2\mu }.$$
(54)
where $`\mu `$ is called the effective mass in analogy with the dispersion relations of massive particles.

This dispersion relation gives rise to a memory function with the asymptotic behavior $`f_\mathrm{m}(\tau )\tau ^{3/2}`$ so that, in contrast to an isotropic bandgap (where $`f_\mathrm{m}(\tau )\tau ^{1/2})`$ we are able to define a finite memory time.

The memory function given by $`f_\mathrm{m}(\tau )\tau ^{3/2}`$ suffers from a singularity at $`\tau =0.`$ This is due to the absence of a high-frequency cutoff, which must be included to model a physical system. We include this high frequency cutoff by considering the slightly modified memory function,

$$f_\mathrm{m}(\tau )=\beta \lambda ^{3/2}\frac{e^{i\left[\delta \tau +\pi /4\right]}}{\left(1+\lambda \tau \right)^{3/2}},$$
(55)
where (shifting to a frame rotating at the atomic frequency) $`\delta =\omega _0\omega _g`$ is the detuning between the band edge and the atomic resonance frequency and $`\beta `$ is a measure of the coupling strength between the atom and the electromagnetic field. The parameter $`\lambda `$ is infinite when there is no upper frequency limit, so for large $`\lambda `$ this memory function is effectively the same as that derived by Vats and John. It is useful to introduce the quantity

$`\gamma `$ $`=`$ $`2\mathrm{}{\displaystyle _0^{\mathrm{}}}𝑑\tau f_\mathrm{m}(\tau )|_{\delta =0}`$
$`=`$ $`2^{3/2}\beta \sqrt{\lambda }`$
which corresponds to the Born-Markov damping rate for this memory function. We use this quantity as a scale against which we can measure the various rates.

For large $`\lambda ,`$ the factor $`\left(1+\lambda \tau \right)^{3/2}`$ has a sharply rising behavior as $`\tau 0`$ on a time scale much faster than the system variables. This is a common feature of physical memory functions and requires special treatment in numerical calculations. We generalize the trapezoidal rule used for the atom in a cavity case by calculating weights that depend on the first two moments of $`\left(1+\lambda \tau \right)^{3/2}`$ and derive an extended two-point integration scheme on a uniform mesh as described in Ref.. Explicitly, we used the weights

$`W_n=\{\begin{array}{cc}w_0^0w_1^0,& n=0\\ (M2)w_0^{M2}+w_1^{M2},& n=M1\\ (n+1)w_0^nw_1^n(n1)w_0^{n1}+w_1^{n1},& \mathrm{otherwise}\end{array}`$

where

$`w_j^n={\displaystyle \frac{1}{\mathrm{\Delta }t^j}}{\displaystyle _{n\mathrm{\Delta }t}^{(n+1)\mathrm{\Delta }t}}𝑑\tau \tau ^j\left(1+\lambda \tau \right)^{3/2}.`$

Since

$`{\displaystyle \underset{n=0}{\overset{n=M1}{}}}W_n={\displaystyle _0^{(M1)\mathrm{\Delta }t}}𝑑\tau \left(1+\lambda \tau \right)^{3/2},`$

this integration scheme treats the semi-singular part of the memory function exactly and allows us to treat a range of non-Markovian situations, including the Markovian limit, within the same algorithm. As mentioned previously we do not consider schemes of higher order than the two-point scheme as the numerical difficulty of evolving the system a single timestep increases dramatically due to the combinatorics introduced by the time-ordering.

Let us first consider the undriven case where an initially excited atom simply undergoes decay inside a bandgap material. Since this case can be treated by an alternative method it serves as another test of the algorithm. In the case of simple decay there is at most only one photon in the electomagnetic field so we can write the state of the atom and field as

$`|\mathrm{\Psi }(t)=a(t)|1|\{0\}+{\displaystyle \underset{k}{}}c_k(t)b_k^{}|0|\{0\},`$

where $`|1`$ and $`|0`$ are the excited and ground states of the atom and $`|\{0\}`$ is the field vacuum. It is not difficult to establish that

$`{\displaystyle \frac{da(t)}{dt}}={\displaystyle _0^t}𝑑sf_\mathrm{m}(ts)a(s),`$

with $`f_\mathrm{m}(ts)`$ given by Eq.(55). This integro-differential equation can be numerically simulated by, for example, the routine (17). In actual fact we used the weighting scheme described above to approximate the integral. The probability of the atom being in the excited state is given by $`P(t)=|a(t)|^2.`$ In Fig. 3 we have compared the results of our algorithm to a direct numerical integration of the above integro-differential equation; the two methods show good agreement. In contrast to the isotropic band gap case, previously treated analytically in Ref. , the decay is always exponential. This can be seen more clearly in the corresponding plot of the logarithm of $`P(t)`$, Fig. 4. Note that in this plot the discrepancy between the two methods is apparent. Our algorithm gets more accurate as $`\lambda `$ is increased as this allows a smaller stepsize to be taken. These plots were made with the parameter value $`\lambda =3\times 10^2\gamma `$ and the variation of the damping rate with detuning from the bandgap is evident. However, when this parameter is increased the situation becomes more and more Markovian, despite the long time tail of the memory function. For $`\lambda =10^5\gamma `$ there is no noticeable difference from the Markovian case; the decay rate is exactly that predicted by the Born-Markov result, and the detuning becomes irrelevant. This is shown in Fig. 5. Note that the $`\delta <0`$ case requires the driving laser to be propagating in a direction that is perpendicular to the direction of $`𝐤_0`$ for the anisotropic bandgap.

For non-zero driving there are no traditional methods available and the results can only be determined by a quantum non-Markovian treatment. In Fig. 6 we have plotted the probability amplitude for the atom to be in its excited state. The driving induces oscillations in this probability and the atom tends to a steady state close to equal probability of being in the excited state and in the ground state due to the reasonably large driving. This behavior is qualitatively equivalent to the Markovian case. The steady-state correlation function, $`C(\tau )=\sigma ^{}(\tau )\sigma _{\mathrm{ss}}\sigma ^{}_{\mathrm{ss}}\sigma _{\mathrm{ss}},`$ is plotted in Fig. 7. Taking the Fourier transform of the steady-state correlation function, $`C(\tau ),`$ gives the so called internal spectrum of the atom, $`S(\omega ).`$ As in the atom in a cavity case this quality allows us to clearly see the effects of the non-trivial structure in the field on the behavior of the atom. In Fig. 8 we have plotted the internal spectrum with $`\lambda =3\times 10^2\gamma `$ for various detunings from the band edge. A definite asymmetry arises in the magnitude of the side peaks directly due to the fact that emission at frequencies inside the gap is prohibited. Note that as $`\lambda `$ is increased this effect becomes less and less until the spectrum becomes equivalent to the Markovian case. In Fig.9 we have plotted the internal spectrum with the same detunings but with $`\lambda =10^5\gamma ;`$ there is no evidence of asymmetry and the values of detuning become irrelevant just as in the case of the undriven atom.

## V Conclusions

A number of different methods of dealing with non-Markovian systems have appeared in the literature recently. These methods can be broadly separated into two main categories: methods that involve extending the system size by the addition of real or imaginary modes  such that the new extended system evolves via a Born-Markov master equation, and methods that attempt to evolve the non-Markovian evolution directly . Note that the approach in Ref.  takes advantage of the special case when the free evolution of the system has an analytical solution to derive a Markov equation for the system. Conceptually, the present work falls into the last category, however, in order to propagate the non-Markovian evolution we have effectively extended the system size by introducing virtual density matrices that follow the paths of the photons in the field before they are either irretrievably emitted or are reabsorbed by the atom. The two approaches have very different origins and the way in which the extended system is constructed is fundamentally different. The present method seems to be slightly more general as no restrictions are made on the form of the memory function, however, time discretization plays a much more fundamental role in our algorithm so it would be very interesting to see on what level these two approaches are related.

The strength of the present algorithm lies in its ability to accurately model weakly non-Markovian systems. An appropriate application would be to determining the validity of the Born-Markov approximation in borderline cases. The biggest constraint on the algorithm is the power law increase in the number of virtual density matrices that need to be stored to propagate the non-Markovian evolution. This is a fundamental constraint and arises from the inherent difficulty in numerically simulating quantum mechanical systems in general. Note that the algorithm presented here was derived without using the algebraic relations of the atomic operators. It is therefore applicable to other system-bath couplings that are linear in the bath operators.

With our algorithm we have determined the evolution of a driven atom radiating into a bandgap material with an anisotropic bandgap; a situation that cannot be treated by any traditional methods. The asymmetry in the internal spectrum of the atom clearly demonstrates the effect of the non-trivial structure of the field on the behavior of the atom. This behavior cannot occur under the Born-Markov approximation. In addition, we have demonstrated the transition to the Born-Markov limit, which occurs despite the presence of a long time tail on the memory function.

In this article we have presented a density matrix version of the algorithm. In this case the system size was so small that a density matrix treatment was appropriate, however, a wave function version can also be derived in a similar fashion. The wavefunction version is applicable to simulating non-Markovian quantum trajectories that also solve for the reduced density matrix . In the non-Markovian case increasing system size is obviously an even more difficult constraint to numerical commutation than in the Markovian case so wave function methods will be important.

## ACKNOWLEDGMENTS

The authors would like to thank M. Collett for useful discussions and M.J. would like to acknowledge the support of the Japan Society for the Promotion of Science. J.H. was supported by the University of Auckland Research Committee and the Marsden Fund of the Royal Society of New Zealand.

## A Non-Markovian model for an atom in a cavity

This appendix shows how the system in section IV A can be written in terms of a non-Markovian model.

Equation (48) is valid in the high-frequency limit : it assumes the rotating wave approximation and that the fundamental frequency $`\nu `$ is much greater than the spread of frequencies about the fundamental frequency such that we can replace the lower bound of the integral over the bath modes by $`\mathrm{}`$. Note that as we are considering a one dimensional field in the continuum limit it is convenient to work in terms of the frequencies rather than the mode labels.

The proof that the two Hamiltonians (47) and (49) are equivalent consists in demonstrating that they are related by a unitary transformation. Let the initial time of interest be $`t=0`$. For the moment let us restrict our consideration to the cavity mode-field system described by the Hamiltonian $`\stackrel{~}{H}`$ \[see Eq.(48)\]. Consider the following scenario: let the coupling between the cavity and the field be time dependent and let it increase slowly such that at a long time in the past $`(t<0)`$ the coupling is turned off and at the initial time $`(t=0)`$ it is equal to $`\kappa `$,

$$\kappa (t)=\{\begin{array}{cc}\kappa e^{ϵt},& t<0\\ \kappa ,& t0\end{array},$$
(A1)
where $`ϵ`$ is a small positive number. The adiabatic limit of an infinitely slow switching on of the coupling is given by the limit where $`e^{ϵt}0`$ as $`t\mathrm{}.`$ and $`ϵ0`$ The Heisenberg equations of motion for the cavity and field modes for $`t<0`$ are given by

$`{\displaystyle \frac{d\stackrel{~}{a}(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\stackrel{~}{a}(t),\stackrel{~}{H}(t)]=(i\nu ϵ)\stackrel{~}{a}(t){\displaystyle \frac{\kappa }{\sqrt{2\pi }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega b(\omega ,t)`$ (A2)
$`{\displaystyle \frac{db(\omega ,t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[b(\omega ,t),\stackrel{~}{H}(t)]=i\omega b(\omega )+{\displaystyle \frac{\kappa }{\sqrt{2\pi }}}\stackrel{~}{a}(t),`$ (A3)
where we have simplified the equations by making the replacement $`\stackrel{~}{a}(t)=a(t)e^{ϵt}.`$ Note that if $`lim_{t_0\mathrm{}}`$ $`a(t_0)`$ is finite this implies that

$$\underset{\mathrm{t}_0\mathrm{}}{lim}\stackrel{~}{a}(t_0)=0.$$
(A4)
Solving Eq.(A3) formally for $`b(\omega ,t)`$ in terms of $`b(\omega ,t_0),`$ where $`t_0<t,`$ and substituting this solution into Eq. (A2) we derive

$$\frac{d\stackrel{~}{a}(t)}{dt}=\left(i\nu +ϵ+\frac{\kappa ^2}{2}\right)\stackrel{~}{a}(t)\frac{\kappa }{\sqrt{2\pi }}_{\mathrm{}}^{\mathrm{}}𝑑\omega b(\omega ,t_0)e^{i\omega (tt_0)},$$
(A5)
where we have used the relation

$$\delta (ts)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑\omega e^{i\omega (ts)}.$$
(A6)
Solving this formally for $`\stackrel{~}{a}(0)`$ in terms of $`\stackrel{~}{a}(t_0)`$ we find

$$\stackrel{~}{a}(0)=\stackrel{~}{a}(t_0)e^{(i\nu +ϵ+\frac{\kappa ^2}{2})t_0}\frac{\kappa }{\sqrt{2\pi }}_{\mathrm{}}^{\mathrm{}}𝑑\omega \stackrel{~}{b}(\omega ,t_0)\frac{1e^{[i(\nu \omega )+ϵ+\frac{\kappa ^2}{2}]t_0}}{\left[i(\nu \omega )+ϵ+\frac{\kappa ^2}{2}\right]}.$$
(A7)
where $`\stackrel{~}{b}(\omega ,t_0)=`$ $`b(\omega ,t_0)e^{i\omega t_0}.`$ Letting $`t_0\mathrm{}`$ and $`ϵ0`$ (and assuming that $`\stackrel{~}{b}(\omega ,t_0)`$ exists in this limit) this solution becomes

$$a(0)=\stackrel{~}{a}(0)=\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^{\mathrm{}}𝑑\omega \stackrel{~}{b}(\omega ,\mathrm{})\frac{\kappa }{\left[i(\omega \nu )\frac{\kappa ^2}{2}\right]}.$$
(A8)
As the system contains no bound states the Hamiltonian in the same limit is given by

$$\stackrel{~}{H}(\mathrm{})=\mathrm{}_{\mathrm{}}^{\mathrm{}}𝑑\omega \omega \stackrel{~}{b}^{}(\omega ,\mathrm{})\stackrel{~}{b}(\omega ,\mathrm{}).$$
(A9)
Therefore, the unitary transformation given by

$$\stackrel{~}{U}=\underset{t_0\mathrm{},ϵ0}{lim}\left[\stackrel{}{T}\mathrm{exp}\left\{\frac{i}{\mathrm{}}_{t_0}^0𝑑s\stackrel{~}{H}(s)\right\}\right].$$
(A10)
diagonalizes the Hamiltonian $`\stackrel{~}{H}(0)`$ and gives $`a(0)`$ in terms of the diagonalizing modes. From Eq.(A4), (A8) and (A9) we see that this is just the transform required to write the Hamiltonian Eq.(47) as Eq.(49) with $`c(\omega )=\stackrel{~}{b}(\omega ,\mathrm{}),`$ thus demonstrating their physical equivalence.