- Email: [email protected]
A Maxwell–Schrödinger solver for quantum optical few-level systems
Computer Physics Communications 182 (2011) 739–747
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
A Maxwell–Schrödinger solver for quantum optical few-level systems ✩ Robert Fleischhaker, Jörg Evers ∗ Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, D-69117 Heidelberg, Germany
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 14 April 2010 Received in revised form 18 August 2010 Accepted 16 October 2010 Available online 28 October 2010 Keywords: Maxwell–Schrödinger equation Quantum optical few-level system Cross-coupling Magnetic field component Second order Lax–Wendroff algorithm Adams predictor
The msprop program presented in this work is capable of solving the Maxwell–Schrödinger equations for one or several laser fields propagating through a medium of quantum optical few-level systems in one spatial dimension and in time. In particular, it allows to numerically treat systems in which a laser field interacts with the medium with both its electric and magnetic component at the same time. The internal dynamics of the few-level system is modeled by a quantum optical master equation which includes coherent processes due to optical transitions driven by the laser fields as well as incoherent processes due to decay and dephasing. The propagation dynamics of the laser fields is treated in slowly varying envelope approximation resulting in a first order wave equation for each laser field envelope function. The program employs an Adams predictor formula second order in time to integrate the quantum optical master equation and a Lax–Wendroff scheme second order in space and time to evolve the wave equations for the fields. The source function in the Lax–Wendroff scheme is specifically adapted to allow taking into account the simultaneous coupling of a laser field to the polarization and the magnetization of the medium. To reduce execution time, a customized data structure is implemented and explained. In three examples the features of the program are demonstrated and the treatment of a system with a phase-dependent cross coupling of the electric and magnetic field component of a laser field is shown. Program summary Program title: msprop Catalogue identifier: AEHR_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEHR_v1_0.html Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 507 625 No. of bytes in distributed program, including test data, etc.: 10 698 552 Distribution format: tar.gz Programming language: C (C99 standard), Mathematica, bash script, gnuplot script Computer: Tested on x86 architecture Operating system: Unix/Linux environment RAM: Less than 30 MB Classification: 2.5 External routines: Standard C math library, accompanying bash script uses gnuplot, bc (basic calculator), and convert (ImageMagick) Nature of problem: We consider a system of quantum optical few-level atoms exposed to several near-resonant continuous-wave or pulsed laser fields. The complexity of the problem arises from the combination of the coherent and incoherent time evolution of the atoms and its dependence on the spatially varying fields. In systems with a coupling to the electric and magnetic field component the simultaneous treatment of both field components poses an additional challenge. Studying the system dynamics requires solving the quantum optical master equation coupled to the wave equations governing the spatio-temporal dynamics of the fields [1,2]. Solution method: We numerically integrate the equations of motion using a second order Adams predictor method for the time evolution of the atomic density matrix and a second order Lax–Wendroff scheme for iterating the fields in space [3]. For the Lax–Wendroff scheme, the source function is adapted such
✩ This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/ science/journal/00104655). Corresponding author. E-mail addresses: robert.fl[email protected] (R. Fleischhaker), [email protected] (J. Evers).
*
0010-4655/$ – see front matter doi:10.1016/j.cpc.2010.10.018
©
2010 Elsevier B.V. All rights reserved.
740
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
that a simultaneous coupling to the polarization and the magnetization of the medium can be taken into account. Restrictions: The evolution of the fields is treated in slowly varying envelope approximation [2] such that variations of the fields in space and time must be on a scale larger than the wavelength and the optical cycle. Propagation is restricted to the forward direction and to one dimension. Concerning the description of the atomic system, only a finite number of basis states can be treated and the laser-driven transitions have to be near-resonant such that the rotating-wave approximation can be applied [2]. Unusual features: The program allows the dipole interaction of both the electric and the magnetic component of a laser field to be taken into account at the same time. Thus, a system with a phasedependent cross coupling of electric and magnetic field component can be treated (see Section 4.2 and [4]). Concerning the implementation of the data structure, it has been optimized for faster memory access. Compared to using standard memory allocation methods, shorter run times are achieved (see Section 3.2). Additional comments: Three examples are given. They each include a readme file, a Mathematica notebook to generate the C-code form of the quantum optical master equation, a parameter file, a bash script which runs the program and converts the numerical data into a movie, two gnuplot scripts, and all files that are produced by running the bash script. Running time: For the first two examples the running time is less than a minute, the third example takes about 12 minutes. On a Pentium 4 (3 GHz) system, a rough estimate can be made with a value of 1 second per million grid points and per field variable. References: [1] Z. Ficek, S. Swain, Quantum Interference and Coherence: Theory and Experiments, Springer, Berlin, 2005. [2] M.O. Scully, M.S. Zubairy, Quantum Optics, Cambridge University Press, Cambridge, 1997. [3] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes, Cambridge University Press, Cambridge, 1992. [4] R. Fleischhaker, J. Evers, Phys. Rev. A 80 (2009) 063816.
© 2010 Elsevier B.V. All rights reserved.
1. Introduction In quantum optics one usually deals with a system of fewlevel atoms exposed to one or several driving laser fields. In a semi-classical description the dynamics of the atomic density matrix is then given by a quantum optical master equation [1]. It combines the coherent evolution due to the laser fields with the incoherent evolution due to spontaneous emission processes. The interaction with the laser fields is formulated in rotating-wave approximation and depends on the classical field amplitudes which are represented by the corresponding Rabi frequencies in the master equation [2]. In cases where the laser fields have an approximately constant intensity over the spatial region of the atomic sample, the dynamics of the system can be described by the master equation alone. But if spatial variations of the laser fields play a role, such as for pulsed laser fields, strong absorption or gain, in addition the wave equation governing the spatio-temporal dynamics of the laser fields has to be taken into account. For the description of typical quantum optical experiments, it is usually sufficient to treat the propagation in one dimension and to use the slowly varying envelope approximation [2]. This approximation reduces the general wave equation which is of second order in time and space to a set of first order differential equations, one for each laser field envelope. The master equation for the atomic density matrix together with the first order wave equations form a closed set of equations, the Maxwell–Schrödinger equations (MSE). When the initial conditions for the atomic density matrix and the boundary conditions for the laser fields are specified, the MSE define a unique solution for the spatio-temporal evolution of the system. In some simple cases such a solution can be expressed analytically. For example, in a medium of two-level atoms exists a group of soliton solutions that rely on self-induced transparency [3]. Similarly, in a three-level system under conditions of electromagnetically induced transparency (EIT) [4], the propagation of a weak, adiabatic probe pulse can be described analytically [5]. A more general class of so-
lutions was found also in the adiabatic regime of an EIT system, which holds for a strong probe field as well [6,7]. However, for systems where additional atomic states have to be considered and/or the driving fields exhibit a more complex behavior, an analytical treatment can easily be rendered impossible. This holds especially true for a new class of systems which received considerable interest in recent years [8–12]. These systems are characterized by a simultaneous interaction with the electric and magnetic component of the same laser field. Here, the additional coupling to the magnetic component adds a new degree of freedom to the propagation dynamics. For example, we have studied a system in which a combination of a closed interaction loop [13–16] and a simultaneous coupling of electric and magnetic field component leads to a dependence on the relative phase of the control fields. Changing this phase during propagation time gives rise to a pronounced spatio-temporal dependence of the fields [12]. To study the dynamics of such systems, one is dependent on a numerical solution of the MSE. Thus, we have implemented a numerical algorithm that integrates the MSE in one space dimension and in time. We use a Lax–Wendroff type [17] method for the space integration in combination with an Adams predictor scheme for the time evolution. We have specifically adapted the Lax–Wendroff scheme, such that a simultaneous coupling to the polarization and to the magnetization of the medium can be taken into account. The communication is organized as follows. In the next section (Section 2) we present the MSE and bring them into a dimensionless form suitable for the numerical treatment (Sections 2.1 and 2.2). We address the basic integration steps and give the explicit formulas to evolve the equations in space and time (Section 2.3). In Section 3 we present our program msprop. First, we discuss the structure of the program (Section 3.1) and the way a time efficient memory usage is realized (Section 3.2). Then, we shortly document the structure of the input and output file (Sections 3.3 and 3.4). In Section 4 we give three examples for a system in which the msprop program is used to calculate a numerical solution. The first example propagates a slow light pulse in an EIT medium. For
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
741
In addition to the coherent laser driving, we also consider spontaneous emission, which in Born–Markov approximation leads to the master equation [1]
∂t =
1 ih¯
[ H , ] + L ,
(3a)
L = L31 + L32 + L34 + L21 + L42
+ L54 + L51 , Fig. 1. (Color online.) Examples for quantum optical few-level systems discussed in this communication. The atomic eigenstates are represented by the black horizontal bars which are labeled with |1–|5. The arrows represent laser fields driving optical transitions between the eigenstates. Red solid arrows denote the probe fields, whereas blue dotted arrows indicate coupling fields. Our analysis focuses on setups in which both the electric (Rabi frequency Ω E ) and the magnetic (Ω B ) component of a single laser field couple to the atoms, as in example (b). Ω1 , Ω2 and ΩC are Rabi frequencies of auxiliary laser fields coupling only with their respective electric components. The discussed program can easily be adopted for other level schemes.
testing purpose, the numerically calculated data is compared to an analytical solution. The second example deals with an EIT system in which a switching of one laser field is used to store and retrieve a second pulsed laser field (Section 4.1). In the third example a system with a simultaneous coupling to the electric and magnetic component of a laser field is treated. Here, the relative phase of the additional laser fields is used to switch the refractive index during propagation time (Section 4.2). Finally, Section 5 summarizes our results.
The numerical method introduced in the following applies to general atomic few-level schemes and laser field couplings. Two examples are shown in Fig. 1. To illustrate the notation and the abstract set of equations solved in the numerical part, we start by considering the specific example in Fig. 1(b) in Section 2.1, and then proceed with an abstract description of the general case in Section 2.2. 2.1. Motivation of the problem The free atomic Hamiltonian for the system in Fig. 1(b) consisting of n = 5 states with energy h¯ ωl (l ∈ {1, . . . , n}) is given by [1, 12]
Ha =
5
h¯ ωl S ll ,
with S lm = |lm|. The atom–field interaction Hamiltonian in rotating-wave approximation can be written as
h¯
Ω B S 21 e −i ν B t + Ω E S 34 e −i ν E t + ΩC S 32 e −i νC t + Ω1 S 51 e −i ν1 t + Ω2 S 54 e −i ν2 t + H.c. . 2
(2)
Here, νl and Ωl denote the frequency and the Rabi frequency of field l ∈ {1, 2, B , E , C }. Note that in this example, all Rabi frequencies except for Ω B are electric field components coupling to the atoms in electric dipole approximation, and Ω E is the coupling of the electric component of the probe field. In addition, Ω B describes the coupling of the magnetic component of the probe field to a magnetic dipole transition of the atom. This generalization of the electric dipole approximation is motivated by the assumption that the atoms have suitable near-degenerate electric and magnetic dipole transitions such that both components of the probe field can couple [12]. In contrast, only the electric dipole coupling is considered for the other fields, for which there are typically only near-resonant electric dipole transitions.
(3c)
with total Hamiltonian H = H a + H I . Here, Llm is the standard Lindblad form for spontaneous emission from state |l to |m, and these equations can easily be generalized to also include dephasing of incoherent pump processes. The propagation of the electric (E) and magnetic (B) field component of a single laser field is described by wave equations following from Maxwell’s equation,
2 ∂ E = μ0 ∂t2 P + μ0 ∂t ∇ × M, c2 t 1 2 − 2 ∂t B = μ0 M − μ0 ∂t ∇ × P, −
1
(4a) (4b)
c
with polarization P and magnetization M acting as source terms. Specializing to a one-dimensional propagation in positive zdirection, we can write X ∈ E, B, P, M} as,
1 2
X 0 ( z, t )e X e i (k0 z−ω0 t ) + c.c.,
(5)
with polarization e X , and wave number of the carrier wave in vacuum k0 with ω0 = ck0 as the corresponding frequency. Assuming that variations in space and time of the envelope function are on a much larger scale than the wavelength λ0 = 2π /k0 and the oscillation period T 0 = 2π /ω0 of the carrier wave, the slowly varying envelope approximation can be applied. In this approximation derivatives of the envelope function in space and time are neglected compared to derivatives of the carrier wave. As a consequence, SVEA reduces the wave equations to first order equations for the envelope functions [2],
1 ik0 k0 ∂z + ∂t E 0 ( z, t ) = P 0 ( z, t ) ∓ M 0 ( z, t ), c 2ε0 2ε0 c 1 ik0 μ0 k0 ∂z + ∂t B 0 ( z, t ) = M 0 ( z, t ) ± P 0 ( z, t ). c
(1)
l =1
HI = −
Llm = γlm (2S ml S lm − S lm S ml − S lm S ml ),
X=
2. Theoretical background
(3b)
2
2ε0 c
(6a) (6b)
Throughout the light propagation, the atoms described by Eqs. (3) act as source terms for the wave propagation equations (6). Using the relations P 0 = Nd E ρ E and M 0 = Nd B ρ B , where N is the atom density, d E [d B ] is the dipole moment of the transition coupling to the electric [magnetic] probe field component, and ρ E [ρ B ] the density matrix element describing the coherence on the electric [magnetic] probe field transition, one can rewrite Eqs. (6) as
1 ∂z + ∂t Ω E = ξ(i ρ E ∓ αρ B ), c 1 ∂z + ∂t Ω B = ξ i α 2 ρ B ± αρ E ,
c
(7a) (7b)
with a constant ξ to be specified later. Here, we have assumed that the coupling strength of the magnetic probe field transition is suppressed by a factor of the fine structure constant α compared to the electric probe field transition. The sign of the magnetization contribution to the electric component and the polarization contribution to the magnetic component each depends on the helicity of the field such that the two different signs apply for leftor right-circular polarization [12].
742
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
It can be seen that the source terms on the right-hand side of Eqs. (6) are a linear combination of probe field coherences. The coupling of light and matter arises since these coherences are determined by Eqs. (3), which in turn depends on the Rabi frequencies determined by Eqs. (7). Together, Eqs. (3) and (7) constitute a closed set of equations that combine the time evolution of the atomic degrees of freedom with a one-dimensional propagation of the electromagnetic fields. As a result, we deal with a problem with two coordinates, the time t and the position z. For the numerical treatment, we define dimensionless time t and position z by
t = tγ ,
(8a)
z = zγ / v ,
(8b)
where 1/γ is a characteristic time scale for the atomic evolution and v is the expected propagation velocity. More specifically, for γ we used the fastest decay rate appearing in the respective level scheme and v can usually be estimated from an analytical formula for the group velocity (see examples in Section 4). On the basis of Eqs. (8) all other physical quantities are scaled to numerical units and Eqs. (3) and (7) become dimensionless equations. Since in what follows, we only deal with dimensionless quantities, we drop the primed notation for simplicity. Then, Eqs. (7) can be written as [12]
[∂z + u ∂t ]Ω E = g E (),
(9a)
[∂z + u ∂t ]Ω B = g B (),
(9b)
and
g B () = η i α 2 B ± α E .
(10)
Here, η = 3π L v /λγ is the dimensionless coupling constant, L = N λ3 /4π 2 , N is the density, and λ is the transition wave length. Note that in the system in Fig. 1(a) with only the electric component of the probe field coupling to the atoms, the corresponding propagation equation is simply
[∂z + u ∂t ]Ω E = i η E .
The numerical integration of the EOM [Eqs. (12) and (13)] is performed on a two-dimensional grid with M grid points in space direction and N grid points in time direction. The value of the density matrix and the fields at a grid point is denoted by
kj := ( j z, kt ) and Ω kj := Ω( j z, kt ),
(14)
where j and k are the indices for space and time direction, z and t are the respective stepsizes, and for notational clarity, we exemplary describe the descritization scheme for only one Rabi frequency Ω . In the same way as for Ω , we denote the space and time dependence of the RHS of the master equation by f kj and of the source function in the wave equation by g kj . The master equation (12) is evolved in time by a standard Adams predictor formula of second order accuracy. This linear multistep method has the advantage that it refers to known values at previous steps, rather than calculating intermediate steps such as in Runge–Kutta methods. This is of particular importance here, since the calculation of intermediate steps requires the values of the light fields at intermediate steps, which are priori are unknown. The Adams predictor formula can be written as,
kj +1 = kj +
t 2
3 f kj − f kj −1 + O t 3 .
(15)
For the first time step (k = 1), the values of kj −1 and Ω kj −1 are outside the numerical grid boundaries and we use a first order formula instead,
2j = 1j + t f j1 + O t 2 .
where the source functions take the form
g E () = η(i E ∓ α B )
2.3. Discretization scheme
(11)
(16)
To evolve the wave equation, we use a second order Lax– Wendroff scheme [17], modified to also include a source term [18]. First, Eq. (13) is discretized with a leap frog method centered around the grid point ( j + 1/2, k),
Ω kj+1 − Ω kj z
k+1/2
= −u
k−1/2
Ω j +1/2 − Ω j +1/2 t
k−1/2
(17)
k+1/2
Then, the unknown grid points Ω j +1/2 and Ω j +1/2 are calculated from a Lax scheme centered at ( j , k − 1/2) and ( j , k + 1/2), k−1/2
Ω j +1/2 − 12 (Ω kj −1 + Ω kj )
2.2. General equations of motion
+ g kj+1/2 .
z /2
= −u
k+1/2
Ω j +1/2 − 12 (Ω kj + Ω kj +1 )
Ω kj − Ω kj −1 t Ω kj +1 − Ω kj
k−1/2
+ gj
,
(18a)
.
(18b)
We now leave to specific examples of Fig. 1 and proceed with an abstract formulation of the modeled system to allow for an easy adaptation to other level schemes. As in Eqs. (3), in the general case, the time evolution of the atomic degrees of freedom is governed by a master equation, which in a compact dimensionless form can be written as,
The unknown grid points at half integer steps of the source function g in Eqs. (17) and (18) are replaced by the mean value of the nearest neighbors,
∂t ( z, t ) = f , {Ω} ,
g kj+1/2 =
(12)
where {Ω} denotes the set of dimensionless Rabi frequencies of the different applied fields. We use the function f to abbreviate the right-hand side (RHS) of the master equation. Its exact form depends on the respective level scheme and is not relevant for the algorithm. The wave equations for the respective dimensionless Ω become [see Eqs. (9)]
[∂z + u ∂t ]Ω = g (),
(13)
where u = v /c is the expected velocity divided by the velocity of light c.
z /2
k−1/2
gj
k+1/2
gj
= =
1 2 1 2 1 2
= −u
t
k+1/2
+ gj
g kj + g kj+1 ,
g kj −1 + g kj ,
g kj + g kj +1 .
(19a) (19b) (19c)
Note that this averaging introduces numerical dissipation. For example, the average in Eq. (19a) can be rewritten as,
1
1 k g kj + g kj+1 = g kj+1/2 + g j +1 − 2g kj+1/2 + g kj 2 2 k 1 k = g j +1/2 + z2 ∂z2 g j +1/2 + O z4 . 2
(20)
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
Here, (∂z2 g )kj+1/2 denotes the second derivative of g with respect to z. Thus, the replacement of g at intermediate values by averages in Eqs. (19) effectively adds a numerical dissipation contribution with diffusion constant z2 /2. The dissipation is restricted to small wavelengths, and can be neglected for the solution of the differential equation on scales much larger than z. On the other hand, the averaging has the advantage that it can lead to better numerical stability, see Chapter 19.1 in [17]. Using Eqs. (17), (18), and (19), we finally arrive at the formula for a space step,
Ω kj+1 = Ω kj + + +
uz
2t
Ω kj −1 − Ω kj +1
uz
Ω kj +1 − 2Ω kj + Ω kj −1
t
u z2
4t
g kj −1 − g kj +1 +
z 2
g kj + g kj+1 .
(21)
+1 At the final time k = N, Ω N and Nj +1 are not known. Here, j Eq. (21) cannot be applied and we use a nearly identical scheme centered around ( j + 1/2, N − 1/2) to derive
Ω Nj+1 =
1 1+
+
u z 2t
z 2
Ω Nj +
uz 2t
g Nj − g Nj+1
.
Ω Nj +−11 + Ω Nj −1 − Ω Nj
(22)
At the initial time k = 1, Ω kj −1 is outside the numerical grid boundaries. To still use Eq. (21) we replace it by the known vacuum value which, at the initial time, is usually zero. We conclude the discussion of the numerical implementation by noting that our numerical implementation is of second order precision in space and time [17]. 3. msprop program 3.1. Program structure The msprop program uses the formulas we have derived in the last section to iteratively calculate the numerical solution for the density matrix and the fields. It is written in C code and uses the C99 standard to include variables of complex type. In Fig. 2 we give an overview of the program structure. First, the program initializes all data structures. It reads the necessary parameters from the input file, allocates memory, opens an output file for each storage time, and sets the boundary conditions for the fields. Then, the density matrix kj and the fields Ω kj are calculated line-byline. The algorithm evolves along the time direction, taking time step k = 1 . . . N for a fixed space index j. This is what we refer to as a time line. After one time line is finished, the algorithm proceeds to the next space index. In the first time line, the values of 1k for k = 2 . . . N can be calculated from the initial condition 11
using only Eq. (15), since Ω1k for k = 1 . . . N is given by the boundary conditions of the fields. After the first time line is calculated, it is written to the output files. All following time lines are calculated in a loop over the space direction (see dark colored frame in Fig. 2). At time index k = 1, the density matrix 1j is given by the initial condition. Then, Ω 1j can be calculated from the last time line in a space step using Eq. (21). All following values of Ω kj are calculated in a second, nested loop (see light colored frame in Fig. 2). Since kj is necessary for the source function, it has to be calculated from the just gained value of Ω kj −1 and kj −1 using Eq. (15). In this manner, Eqs. (21) and (15) are used in alternation until the
743
last space step at k = N is performed by Eq. (22). After the calculation of a time line is finished, it is written to the output files if a storage interval is reached. Once the loop over the space direction is finished, the final time line is also written to an additional output file. 3.2. Memory structure The data structures for the density matrix and the fields each consist of a three-dimensional array. The first index represents the time step, the second index distinguishes between the last two space steps, and the third index runs over the elements of the density matrix or the field envelopes. Since the number of time steps depends on the values given in the input file, the actual size of the data structure is not known at compilation time such that it has to be allocated dynamically. Rather than allocating memory for each time step and both stored space steps separately, this is done in one block. This is advantageous because then, the data for one calculation step is close together in the physical memory which decreases memory access time. It requires, however, to set up the corresponding pointer structure by hand. The scheme is shown in Fig. 3. Each of the three layers of the data structure is allocated in one block and afterwards the connecting pointers are set. In the algorithm, further speed optimization is achieved by interchanging pointers rather than copying data once the new values for the elements of the density matrix and the field envelopes have been calculated. The new values are stored into the data structure addressed by the second index equal to 1 which overwrites the one before last values in space direction. Then, the pointers with second index 0 and 1 are interchanged, updating the data structure to the new values in the correct order for the next space step. 3.3. Input file The input file is a plain text file, used to set the physical and numerical parameters necessary for the numerical integration. Its name can be freely chosen and must be provided as the first command line argument when the msprop program is executed. The input file starts with an initial line describing the purpose of the file. Then, for each block of physical or numerical parameters a single comment line is followed by the actual data. First, the different decay rates, detunings, and field strength values are given. In the case of the system with a cross coupling of the electric and magnetic field component (see the third example, Section 4.3), there are additional values for the multi-photon phase, the helicity of the probe field, and the parameters of the dark state. Then, the peak and switch times, the respective width of pulsed fields, and the type of pulse envelope function are given. These parameters fix the initial conditions of the fields at the medium entry. It follows the value for L, necessary for the dimensionless coupling constant in the source function (see Section 2.2). Next, the parameters determining the numerical grid are specified, namely the maximum propagation time and distance as well as the numerical time and space step. The number of time and space steps are calculated in the program. But to reduce the amount of data, additional values for the number of time and space steps written to the output files are set. Finally, a frames-per-second value for the conversion to a movie file by the accompanying bash script is given. 3.4. Output files The output files store the calculated numerical data. They are plain text files where the name consists of a first part provided by the second command line argument and a second part of the form “_t = time.dat”. The program writes one output file for each
744
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
Fig. 2. (Color online.) Overview of the msprop program structure. The program flow is represented by the arrows and the boxes combine related processes in groups. The two larger frames represent the loop for the space direction nested with the loop for the time direction. Conditional statements are represented by diamond-shaped boxes.
Fig. 3. (Color online.) Memory structure for the density matrix and the fields. The structure consists of three arrays (light, medium and dark colored frame) which are connected by pointers (arrows). Data cells with one index contain a pointer to an array of pointers (R [i ]), with two indices a pointer to an array (R [i ][ j ]), and with three indices the actual data (R [i ][ j ][k]).
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
745
Fig. 4. (Color online.) First example, time line at the medium exit. The numerical probe pulse (solid red line) is compared to the analytic solution (dashed green line) for constant control field (dotted blue line). All fields are normalized to their initial value and the dimensionless time has been used.
storage time containing the numerical data for all storage points in space direction. These output files are written into a “data” subfolder. Another single output file with the name “timeline.dat” is written into the main folder. It contains the numerical data at the medium exit (last grid point in space direction) for all storage times. The structure of the output files is column based. The output files in the “data” subfolder contain the numerical value of the space grid point in the first column. Similarly the “timeline.dat” file contains the value of the time grid points in the first column. The other columns contain the values of the fields and the elements of the density matrix in the order of the corresponding data structures. Since this order depends on the name conventions of the respective quantum optical level scheme, a short specification of this data structures is given in the beginning of each source code file. 4. Examples With this manuscript, we provide three examples to illustrate the capabilities of the msprop program. The first example deals with a slow light pulse in an EIT system with constant control field. It can be solved analytically to a very good approximation and serves as a test case for the correct functioning of the program. The second example takes the same system to a different parameter range. Here, the numerical solution displays the full spatio-temporal dynamics of the light storage mechanism which in the general parameter regime is not accessible to analytical solutions. Finally, the third example deals with a system in which both the electric and magnetic components of the same probe field couple to the medium. It demonstrates how both components can be taken into account simultaneously in the numerical treatment. The numerical solution shows the specific phase-dependent dynamics which result from a cross-coupling of the electric and magnetic component. For all examples we show a plot of the time line at the medium exit, the numerical data from the last grid point in space direction and for all stored grid points in time direction. These plots can be generated directly from the numerical data with the gnuplot scrips “gp_timeline.sc”.
4.1. Slow light pulse The dynamics of a slow light pulse propagating in a three-level EIT medium can be solved analytically in the limit of a weak, adiabatic probe pulse [4,5]. We use the analytic formula for the group velocity to fix the numerically expected propagation velocity. This leads to the following expressions for the dimensionless propagation velocity u and the dimensionless coupling η [see Eqs. (13) and (11)]
u= 1+
3L γ p ω
Ωc2
−1 (23)
,
3
η = L ωu ,
(24)
2
where L = N λ3 /4π 2 is defined by the density N and the transition wavelength λ, and γ p , ω , Ωc are the probe field decay rate, probe field frequency, and the control field Rabi frequency. Here, and in the following, we express all numerical quantities in units of γ , the overall decay rate of the excited state. For a Gaussian probe pulse the dimensionless solution for the envelope function reads
Ω p ( z, t ) = Ω p
(t − z)2 σ , exp − σ˜ 2σ˜ 2
(25)
where Ω p and σ are the initial probe field amplitude and temporal width and σ˜ 2 = σ 2 + 4z/Ωc2 is the propagation dependent width of the pulse. To compare this solution to numerical data generated by the msprop program, we use a probe field amplitude of Ω p = 0.01, a control field amplitude of Ωc = 2, and a pulse length of σ = 50 to propagate the pulse over a length of z = 400. With these parameters the regime of a weak adiabatic probe field is well fulfilled. In Fig. 4 we show the data from the time line at the medium exit and compare it to the analytic solution. The exit time of the pulse is given by adding the time when the pulse enters the medium (t = 200) and the propagation time through the medium. Due to the scaling, the numerical velocity z/t = 1 and the propagation time through the medium (t = 400) is numerically
746
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
Fig. 5. (Color online.) Second example, time line at the medium exit. The probe field (solid red line) and the control field (dashed blue line) are shown. Both fields are scaled to their initial value and the dimensionless time has been used.
equivalent to the medium length (z = 400). Thus, the exit time of the pulse is given by t = 200 + 400 = 600 which is confirmed by the data shown in Fig. 4. The complete propagation dynamics can be studied with the “movie.avi” file which is generated by the accompanying bash script “run.bash”. 4.2. Storage of light The second example features the storage of a light pulse by switching the control field during the propagation time of a slow light pulse in an EIT medium. Since the EOM are the same as in the first example, the same expressions for u and η apply. To demonstrate the evolution of the spin coherence, we use a probe field strength of Ω p = 1.5 which is comparable to the control field strength of Ωc = 2 and a small ground state decoherence of γdec = 0.0025. The medium is now z = 150 long and the pulse with a width of σ = 20 enters the medium at t = 75. The control field is switched off at t = 150 and switched back on at t = 300 such that we have a storage time of t = 150. Since the exit time of the pulse is given by the sum of the entering time (t = 75), the propagation time (t = 150), and the storage time (t = 150), the retrieved pulse exits the medium at t = 375. In Fig. 5 we show the time line at the medium exit for the probe and control field. One can clearly see the signature of the energy transfer from the probe field to the control field at the moment where the pulse enters the medium (peak on top of control field at t = 75), the switch off time of the control field, and the adiabaton dip in the control field accompanying the pulse. The complete dynamics of the storage process can be studied with the movie file generated by the bash script. Here, in addition to the probe and control field, the build up and decay of the spin coherence are shown. 4.3. Phase controlled propagation in a system with a coupling of the electric and magnetic field component In the third example we deal with a system that exhibits a simultaneous coupling to the electric and magnetic component of a laser field. It serves to demonstrate the capabilities of the msprop
program to take into account both components and their coupling to the polarization and magnetization of the medium. In particular, it shows how to implement the modified Lax–Wendroff source function [Eq. (10)] into the numerical algorithm. A detailed discussion of the example system and its theoretical modeling are given in [12]. There, it was shown that the group velocity for a pulsed probe field is in principle given by an EIT-type expression and we find for the numerically expected propagation velocity,
u= 1+
3L γ34 ω44
Ωc2
−1 ,
(26)
where γ34 is the decay rate on the electric probe field transition and 44 is the population fraction in the lower state of the electric probe field transition. The scaling of all time units is done with the overall decay rate out of the upper state of the electric probe field transition. The physical parameters are closed to the example which is discussed in [12]. A Gaussian probe pulse of temporal width σ = 50 enters a z = 1200 long medium at t = 400. We assumed a leftcircular polarization of the probe field (helicity value +1) and an initial multi-photon phase of π /2. During the propagation time at t = 800 the multi-photon phase is switched to zero and at t = 1200 back to π /2. This changes the effective refractive index of the probe field. While in the beginning and at the end the probe field accumulates phase, it experiences gain during the intermediate time. In Fig. 6 we show the time line at the medium exit for the probe field and the multi-photon phase. One can clearly see the switching of the multi-photon phase as well as the change in phase and increase in amplitude of the probe pulse. The exit time (t = 1600) is given by the sum of the entering time and the propagation time (t = 1200). For the complete dynamics, we refer to the movie file provided with the example. 5. Summary We have implemented a numerical integration of the Maxwell– Schrödinger equations capable of solving the propagation dynamics of several laser fields through a medium of quantum optical fewlevel atoms in time and in one space dimension. In particular, our
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
747
Fig. 6. (Color online.) Third example, time line at medium exit. The absolute value (thick solid black line), the real part (dashed blue line), and the imaginary part (solid red line) are shown. Their value is scaled to the initial value of the probe field on the left vertical axis and the dimensionless time has been used. The phase of the probe pulse (straight dashed green line) and the multi-photon phase (dash-dotted cyan line) are also shown and their value can be read off from the right vertical axis. For a clear view the probe field phase is set to zero at times where the absolute value of the probe field is below one percent of its peak value.
program is able to deal with systems where a simultaneous coupling of the electric and magnetic field component of the same laser field has to be taken into account. The program employs an Adams predictor method together with a Lax–Wendroff scheme to perform the numerical time and space step. For the Lax–Wendroff scheme the source function has been suitably adapted to take into account a coupling to the polarization and the magnetization of the medium. In addition to a detailed description of the algorithm and the program structure we have provided three examples which illustrate the capabilities of our program. Acknowledgements We would like to thank T.N. Dey for helpful discussions. References [1] Z. Ficek, S. Swain, Quantum Interference and Coherence: Theory and Experiments, Springer, Berlin, 2005.
[2] M.O. Scully, M.S. Zubairy, Quantum Optics, Cambridge University Press, Cambridge, 1997. [3] S.L. McCall, E.L. Hahn, Phys. Rev. Lett. 18 (1967) 908; S.L. McCall, E.L. Hahn, Phys. Rev. 183 (1969) 457. [4] S.E. Harris, Phys. Today 50 (1997) 36. [5] M. Fleischhauer, A. Imamoglu, J.P. Marangos, Rev. Mod. Phys. 77 (2005) 633. [6] R. Grobe, F.T. Hioe, J.H. Eberly, Phys. Rev. Lett. 73 (1994) 3183. [7] J.R. Csesznegi, R. Grobe, Phys. Rev. Lett. 79 (1997) 3162. [8] M.Ö. Oktel, Ö.E. Müstecaplıo˘glu, Phys. Rev. A 70 (2004) 053806. [9] Q. Thommen, P. Mandel, Phys. Rev. Lett. 96 (2006) 053601. [10] J.Q. Shen, S. He, Ann. Phys. 16 (2007) 364. [11] J. Kästel, M. Fleischhauer, S.F. Yelin, R.L. Walsworth, Phys. Rev. Lett. 99 (2007) 073602. [12] R. Fleischhaker, J. Evers, Phys. Rev. A 80 (2009) 063816. [13] S.J. Buckle, S.M. Barnett, P.L. Knight, M.A. Lauder, D.T. Pegg, Opt. Acta 33 (1986) 1129. [14] D.V. Kosachiov, B.G. Matisov, Y.V. Rozhdestvensky, J. Phys. B 25 (1992) 2473. [15] M. Mahmoudi, J. Evers, Phys. Rev. A 74 (2006) 063827. [16] R. Fleischhaker, J. Evers, Phys. Rev. A 77 (2008) 043805. [17] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes, Cambridge University Press, Cambridge, 1992. [18] K.R. Hansen, Standing wave electromagnetically induced transparency, Master thesis, University of Aarhus, 2006.
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
A Maxwell–Schrödinger solver for quantum optical few-level systems ✩ Robert Fleischhaker, Jörg Evers ∗ Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, D-69117 Heidelberg, Germany
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 14 April 2010 Received in revised form 18 August 2010 Accepted 16 October 2010 Available online 28 October 2010 Keywords: Maxwell–Schrödinger equation Quantum optical few-level system Cross-coupling Magnetic field component Second order Lax–Wendroff algorithm Adams predictor
The msprop program presented in this work is capable of solving the Maxwell–Schrödinger equations for one or several laser fields propagating through a medium of quantum optical few-level systems in one spatial dimension and in time. In particular, it allows to numerically treat systems in which a laser field interacts with the medium with both its electric and magnetic component at the same time. The internal dynamics of the few-level system is modeled by a quantum optical master equation which includes coherent processes due to optical transitions driven by the laser fields as well as incoherent processes due to decay and dephasing. The propagation dynamics of the laser fields is treated in slowly varying envelope approximation resulting in a first order wave equation for each laser field envelope function. The program employs an Adams predictor formula second order in time to integrate the quantum optical master equation and a Lax–Wendroff scheme second order in space and time to evolve the wave equations for the fields. The source function in the Lax–Wendroff scheme is specifically adapted to allow taking into account the simultaneous coupling of a laser field to the polarization and the magnetization of the medium. To reduce execution time, a customized data structure is implemented and explained. In three examples the features of the program are demonstrated and the treatment of a system with a phase-dependent cross coupling of the electric and magnetic field component of a laser field is shown. Program summary Program title: msprop Catalogue identifier: AEHR_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEHR_v1_0.html Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 507 625 No. of bytes in distributed program, including test data, etc.: 10 698 552 Distribution format: tar.gz Programming language: C (C99 standard), Mathematica, bash script, gnuplot script Computer: Tested on x86 architecture Operating system: Unix/Linux environment RAM: Less than 30 MB Classification: 2.5 External routines: Standard C math library, accompanying bash script uses gnuplot, bc (basic calculator), and convert (ImageMagick) Nature of problem: We consider a system of quantum optical few-level atoms exposed to several near-resonant continuous-wave or pulsed laser fields. The complexity of the problem arises from the combination of the coherent and incoherent time evolution of the atoms and its dependence on the spatially varying fields. In systems with a coupling to the electric and magnetic field component the simultaneous treatment of both field components poses an additional challenge. Studying the system dynamics requires solving the quantum optical master equation coupled to the wave equations governing the spatio-temporal dynamics of the fields [1,2]. Solution method: We numerically integrate the equations of motion using a second order Adams predictor method for the time evolution of the atomic density matrix and a second order Lax–Wendroff scheme for iterating the fields in space [3]. For the Lax–Wendroff scheme, the source function is adapted such
✩ This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/ science/journal/00104655). Corresponding author. E-mail addresses: robert.fl[email protected] (R. Fleischhaker), [email protected] (J. Evers).
*
0010-4655/$ – see front matter doi:10.1016/j.cpc.2010.10.018
©
2010 Elsevier B.V. All rights reserved.
740
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
that a simultaneous coupling to the polarization and the magnetization of the medium can be taken into account. Restrictions: The evolution of the fields is treated in slowly varying envelope approximation [2] such that variations of the fields in space and time must be on a scale larger than the wavelength and the optical cycle. Propagation is restricted to the forward direction and to one dimension. Concerning the description of the atomic system, only a finite number of basis states can be treated and the laser-driven transitions have to be near-resonant such that the rotating-wave approximation can be applied [2]. Unusual features: The program allows the dipole interaction of both the electric and the magnetic component of a laser field to be taken into account at the same time. Thus, a system with a phasedependent cross coupling of electric and magnetic field component can be treated (see Section 4.2 and [4]). Concerning the implementation of the data structure, it has been optimized for faster memory access. Compared to using standard memory allocation methods, shorter run times are achieved (see Section 3.2). Additional comments: Three examples are given. They each include a readme file, a Mathematica notebook to generate the C-code form of the quantum optical master equation, a parameter file, a bash script which runs the program and converts the numerical data into a movie, two gnuplot scripts, and all files that are produced by running the bash script. Running time: For the first two examples the running time is less than a minute, the third example takes about 12 minutes. On a Pentium 4 (3 GHz) system, a rough estimate can be made with a value of 1 second per million grid points and per field variable. References: [1] Z. Ficek, S. Swain, Quantum Interference and Coherence: Theory and Experiments, Springer, Berlin, 2005. [2] M.O. Scully, M.S. Zubairy, Quantum Optics, Cambridge University Press, Cambridge, 1997. [3] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes, Cambridge University Press, Cambridge, 1992. [4] R. Fleischhaker, J. Evers, Phys. Rev. A 80 (2009) 063816.
© 2010 Elsevier B.V. All rights reserved.
1. Introduction In quantum optics one usually deals with a system of fewlevel atoms exposed to one or several driving laser fields. In a semi-classical description the dynamics of the atomic density matrix is then given by a quantum optical master equation [1]. It combines the coherent evolution due to the laser fields with the incoherent evolution due to spontaneous emission processes. The interaction with the laser fields is formulated in rotating-wave approximation and depends on the classical field amplitudes which are represented by the corresponding Rabi frequencies in the master equation [2]. In cases where the laser fields have an approximately constant intensity over the spatial region of the atomic sample, the dynamics of the system can be described by the master equation alone. But if spatial variations of the laser fields play a role, such as for pulsed laser fields, strong absorption or gain, in addition the wave equation governing the spatio-temporal dynamics of the laser fields has to be taken into account. For the description of typical quantum optical experiments, it is usually sufficient to treat the propagation in one dimension and to use the slowly varying envelope approximation [2]. This approximation reduces the general wave equation which is of second order in time and space to a set of first order differential equations, one for each laser field envelope. The master equation for the atomic density matrix together with the first order wave equations form a closed set of equations, the Maxwell–Schrödinger equations (MSE). When the initial conditions for the atomic density matrix and the boundary conditions for the laser fields are specified, the MSE define a unique solution for the spatio-temporal evolution of the system. In some simple cases such a solution can be expressed analytically. For example, in a medium of two-level atoms exists a group of soliton solutions that rely on self-induced transparency [3]. Similarly, in a three-level system under conditions of electromagnetically induced transparency (EIT) [4], the propagation of a weak, adiabatic probe pulse can be described analytically [5]. A more general class of so-
lutions was found also in the adiabatic regime of an EIT system, which holds for a strong probe field as well [6,7]. However, for systems where additional atomic states have to be considered and/or the driving fields exhibit a more complex behavior, an analytical treatment can easily be rendered impossible. This holds especially true for a new class of systems which received considerable interest in recent years [8–12]. These systems are characterized by a simultaneous interaction with the electric and magnetic component of the same laser field. Here, the additional coupling to the magnetic component adds a new degree of freedom to the propagation dynamics. For example, we have studied a system in which a combination of a closed interaction loop [13–16] and a simultaneous coupling of electric and magnetic field component leads to a dependence on the relative phase of the control fields. Changing this phase during propagation time gives rise to a pronounced spatio-temporal dependence of the fields [12]. To study the dynamics of such systems, one is dependent on a numerical solution of the MSE. Thus, we have implemented a numerical algorithm that integrates the MSE in one space dimension and in time. We use a Lax–Wendroff type [17] method for the space integration in combination with an Adams predictor scheme for the time evolution. We have specifically adapted the Lax–Wendroff scheme, such that a simultaneous coupling to the polarization and to the magnetization of the medium can be taken into account. The communication is organized as follows. In the next section (Section 2) we present the MSE and bring them into a dimensionless form suitable for the numerical treatment (Sections 2.1 and 2.2). We address the basic integration steps and give the explicit formulas to evolve the equations in space and time (Section 2.3). In Section 3 we present our program msprop. First, we discuss the structure of the program (Section 3.1) and the way a time efficient memory usage is realized (Section 3.2). Then, we shortly document the structure of the input and output file (Sections 3.3 and 3.4). In Section 4 we give three examples for a system in which the msprop program is used to calculate a numerical solution. The first example propagates a slow light pulse in an EIT medium. For
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
741
In addition to the coherent laser driving, we also consider spontaneous emission, which in Born–Markov approximation leads to the master equation [1]
∂t =
1 ih¯
[ H , ] + L ,
(3a)
L = L31 + L32 + L34 + L21 + L42
+ L54 + L51 , Fig. 1. (Color online.) Examples for quantum optical few-level systems discussed in this communication. The atomic eigenstates are represented by the black horizontal bars which are labeled with |1–|5. The arrows represent laser fields driving optical transitions between the eigenstates. Red solid arrows denote the probe fields, whereas blue dotted arrows indicate coupling fields. Our analysis focuses on setups in which both the electric (Rabi frequency Ω E ) and the magnetic (Ω B ) component of a single laser field couple to the atoms, as in example (b). Ω1 , Ω2 and ΩC are Rabi frequencies of auxiliary laser fields coupling only with their respective electric components. The discussed program can easily be adopted for other level schemes.
testing purpose, the numerically calculated data is compared to an analytical solution. The second example deals with an EIT system in which a switching of one laser field is used to store and retrieve a second pulsed laser field (Section 4.1). In the third example a system with a simultaneous coupling to the electric and magnetic component of a laser field is treated. Here, the relative phase of the additional laser fields is used to switch the refractive index during propagation time (Section 4.2). Finally, Section 5 summarizes our results.
The numerical method introduced in the following applies to general atomic few-level schemes and laser field couplings. Two examples are shown in Fig. 1. To illustrate the notation and the abstract set of equations solved in the numerical part, we start by considering the specific example in Fig. 1(b) in Section 2.1, and then proceed with an abstract description of the general case in Section 2.2. 2.1. Motivation of the problem The free atomic Hamiltonian for the system in Fig. 1(b) consisting of n = 5 states with energy h¯ ωl (l ∈ {1, . . . , n}) is given by [1, 12]
Ha =
5
h¯ ωl S ll ,
with S lm = |lm|. The atom–field interaction Hamiltonian in rotating-wave approximation can be written as
h¯
Ω B S 21 e −i ν B t + Ω E S 34 e −i ν E t + ΩC S 32 e −i νC t + Ω1 S 51 e −i ν1 t + Ω2 S 54 e −i ν2 t + H.c. . 2
(2)
Here, νl and Ωl denote the frequency and the Rabi frequency of field l ∈ {1, 2, B , E , C }. Note that in this example, all Rabi frequencies except for Ω B are electric field components coupling to the atoms in electric dipole approximation, and Ω E is the coupling of the electric component of the probe field. In addition, Ω B describes the coupling of the magnetic component of the probe field to a magnetic dipole transition of the atom. This generalization of the electric dipole approximation is motivated by the assumption that the atoms have suitable near-degenerate electric and magnetic dipole transitions such that both components of the probe field can couple [12]. In contrast, only the electric dipole coupling is considered for the other fields, for which there are typically only near-resonant electric dipole transitions.
(3c)
with total Hamiltonian H = H a + H I . Here, Llm is the standard Lindblad form for spontaneous emission from state |l to |m, and these equations can easily be generalized to also include dephasing of incoherent pump processes. The propagation of the electric (E) and magnetic (B) field component of a single laser field is described by wave equations following from Maxwell’s equation,
2 ∂ E = μ0 ∂t2 P + μ0 ∂t ∇ × M, c2 t 1 2 − 2 ∂t B = μ0 M − μ0 ∂t ∇ × P, −
1
(4a) (4b)
c
with polarization P and magnetization M acting as source terms. Specializing to a one-dimensional propagation in positive zdirection, we can write X ∈ E, B, P, M} as,
1 2
X 0 ( z, t )e X e i (k0 z−ω0 t ) + c.c.,
(5)
with polarization e X , and wave number of the carrier wave in vacuum k0 with ω0 = ck0 as the corresponding frequency. Assuming that variations in space and time of the envelope function are on a much larger scale than the wavelength λ0 = 2π /k0 and the oscillation period T 0 = 2π /ω0 of the carrier wave, the slowly varying envelope approximation can be applied. In this approximation derivatives of the envelope function in space and time are neglected compared to derivatives of the carrier wave. As a consequence, SVEA reduces the wave equations to first order equations for the envelope functions [2],
1 ik0 k0 ∂z + ∂t E 0 ( z, t ) = P 0 ( z, t ) ∓ M 0 ( z, t ), c 2ε0 2ε0 c 1 ik0 μ0 k0 ∂z + ∂t B 0 ( z, t ) = M 0 ( z, t ) ± P 0 ( z, t ). c
(1)
l =1
HI = −
Llm = γlm (2S ml S lm − S lm S ml − S lm S ml ),
X=
2. Theoretical background
(3b)
2
2ε0 c
(6a) (6b)
Throughout the light propagation, the atoms described by Eqs. (3) act as source terms for the wave propagation equations (6). Using the relations P 0 = Nd E ρ E and M 0 = Nd B ρ B , where N is the atom density, d E [d B ] is the dipole moment of the transition coupling to the electric [magnetic] probe field component, and ρ E [ρ B ] the density matrix element describing the coherence on the electric [magnetic] probe field transition, one can rewrite Eqs. (6) as
1 ∂z + ∂t Ω E = ξ(i ρ E ∓ αρ B ), c 1 ∂z + ∂t Ω B = ξ i α 2 ρ B ± αρ E ,
c
(7a) (7b)
with a constant ξ to be specified later. Here, we have assumed that the coupling strength of the magnetic probe field transition is suppressed by a factor of the fine structure constant α compared to the electric probe field transition. The sign of the magnetization contribution to the electric component and the polarization contribution to the magnetic component each depends on the helicity of the field such that the two different signs apply for leftor right-circular polarization [12].
742
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
It can be seen that the source terms on the right-hand side of Eqs. (6) are a linear combination of probe field coherences. The coupling of light and matter arises since these coherences are determined by Eqs. (3), which in turn depends on the Rabi frequencies determined by Eqs. (7). Together, Eqs. (3) and (7) constitute a closed set of equations that combine the time evolution of the atomic degrees of freedom with a one-dimensional propagation of the electromagnetic fields. As a result, we deal with a problem with two coordinates, the time t and the position z. For the numerical treatment, we define dimensionless time t and position z by
t = tγ ,
(8a)
z = zγ / v ,
(8b)
where 1/γ is a characteristic time scale for the atomic evolution and v is the expected propagation velocity. More specifically, for γ we used the fastest decay rate appearing in the respective level scheme and v can usually be estimated from an analytical formula for the group velocity (see examples in Section 4). On the basis of Eqs. (8) all other physical quantities are scaled to numerical units and Eqs. (3) and (7) become dimensionless equations. Since in what follows, we only deal with dimensionless quantities, we drop the primed notation for simplicity. Then, Eqs. (7) can be written as [12]
[∂z + u ∂t ]Ω E = g E (),
(9a)
[∂z + u ∂t ]Ω B = g B (),
(9b)
and
g B () = η i α 2 B ± α E .
(10)
Here, η = 3π L v /λγ is the dimensionless coupling constant, L = N λ3 /4π 2 , N is the density, and λ is the transition wave length. Note that in the system in Fig. 1(a) with only the electric component of the probe field coupling to the atoms, the corresponding propagation equation is simply
[∂z + u ∂t ]Ω E = i η E .
The numerical integration of the EOM [Eqs. (12) and (13)] is performed on a two-dimensional grid with M grid points in space direction and N grid points in time direction. The value of the density matrix and the fields at a grid point is denoted by
kj := ( j z, kt ) and Ω kj := Ω( j z, kt ),
(14)
where j and k are the indices for space and time direction, z and t are the respective stepsizes, and for notational clarity, we exemplary describe the descritization scheme for only one Rabi frequency Ω . In the same way as for Ω , we denote the space and time dependence of the RHS of the master equation by f kj and of the source function in the wave equation by g kj . The master equation (12) is evolved in time by a standard Adams predictor formula of second order accuracy. This linear multistep method has the advantage that it refers to known values at previous steps, rather than calculating intermediate steps such as in Runge–Kutta methods. This is of particular importance here, since the calculation of intermediate steps requires the values of the light fields at intermediate steps, which are priori are unknown. The Adams predictor formula can be written as,
kj +1 = kj +
t 2
3 f kj − f kj −1 + O t 3 .
(15)
For the first time step (k = 1), the values of kj −1 and Ω kj −1 are outside the numerical grid boundaries and we use a first order formula instead,
2j = 1j + t f j1 + O t 2 .
where the source functions take the form
g E () = η(i E ∓ α B )
2.3. Discretization scheme
(11)
(16)
To evolve the wave equation, we use a second order Lax– Wendroff scheme [17], modified to also include a source term [18]. First, Eq. (13) is discretized with a leap frog method centered around the grid point ( j + 1/2, k),
Ω kj+1 − Ω kj z
k+1/2
= −u
k−1/2
Ω j +1/2 − Ω j +1/2 t
k−1/2
(17)
k+1/2
Then, the unknown grid points Ω j +1/2 and Ω j +1/2 are calculated from a Lax scheme centered at ( j , k − 1/2) and ( j , k + 1/2), k−1/2
Ω j +1/2 − 12 (Ω kj −1 + Ω kj )
2.2. General equations of motion
+ g kj+1/2 .
z /2
= −u
k+1/2
Ω j +1/2 − 12 (Ω kj + Ω kj +1 )
Ω kj − Ω kj −1 t Ω kj +1 − Ω kj
k−1/2
+ gj
,
(18a)
.
(18b)
We now leave to specific examples of Fig. 1 and proceed with an abstract formulation of the modeled system to allow for an easy adaptation to other level schemes. As in Eqs. (3), in the general case, the time evolution of the atomic degrees of freedom is governed by a master equation, which in a compact dimensionless form can be written as,
The unknown grid points at half integer steps of the source function g in Eqs. (17) and (18) are replaced by the mean value of the nearest neighbors,
∂t ( z, t ) = f , {Ω} ,
g kj+1/2 =
(12)
where {Ω} denotes the set of dimensionless Rabi frequencies of the different applied fields. We use the function f to abbreviate the right-hand side (RHS) of the master equation. Its exact form depends on the respective level scheme and is not relevant for the algorithm. The wave equations for the respective dimensionless Ω become [see Eqs. (9)]
[∂z + u ∂t ]Ω = g (),
(13)
where u = v /c is the expected velocity divided by the velocity of light c.
z /2
k−1/2
gj
k+1/2
gj
= =
1 2 1 2 1 2
= −u
t
k+1/2
+ gj
g kj + g kj+1 ,
g kj −1 + g kj ,
g kj + g kj +1 .
(19a) (19b) (19c)
Note that this averaging introduces numerical dissipation. For example, the average in Eq. (19a) can be rewritten as,
1
1 k g kj + g kj+1 = g kj+1/2 + g j +1 − 2g kj+1/2 + g kj 2 2 k 1 k = g j +1/2 + z2 ∂z2 g j +1/2 + O z4 . 2
(20)
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
Here, (∂z2 g )kj+1/2 denotes the second derivative of g with respect to z. Thus, the replacement of g at intermediate values by averages in Eqs. (19) effectively adds a numerical dissipation contribution with diffusion constant z2 /2. The dissipation is restricted to small wavelengths, and can be neglected for the solution of the differential equation on scales much larger than z. On the other hand, the averaging has the advantage that it can lead to better numerical stability, see Chapter 19.1 in [17]. Using Eqs. (17), (18), and (19), we finally arrive at the formula for a space step,
Ω kj+1 = Ω kj + + +
uz
2t
Ω kj −1 − Ω kj +1
uz
Ω kj +1 − 2Ω kj + Ω kj −1
t
u z2
4t
g kj −1 − g kj +1 +
z 2
g kj + g kj+1 .
(21)
+1 At the final time k = N, Ω N and Nj +1 are not known. Here, j Eq. (21) cannot be applied and we use a nearly identical scheme centered around ( j + 1/2, N − 1/2) to derive
Ω Nj+1 =
1 1+
+
u z 2t
z 2
Ω Nj +
uz 2t
g Nj − g Nj+1
.
Ω Nj +−11 + Ω Nj −1 − Ω Nj
(22)
At the initial time k = 1, Ω kj −1 is outside the numerical grid boundaries. To still use Eq. (21) we replace it by the known vacuum value which, at the initial time, is usually zero. We conclude the discussion of the numerical implementation by noting that our numerical implementation is of second order precision in space and time [17]. 3. msprop program 3.1. Program structure The msprop program uses the formulas we have derived in the last section to iteratively calculate the numerical solution for the density matrix and the fields. It is written in C code and uses the C99 standard to include variables of complex type. In Fig. 2 we give an overview of the program structure. First, the program initializes all data structures. It reads the necessary parameters from the input file, allocates memory, opens an output file for each storage time, and sets the boundary conditions for the fields. Then, the density matrix kj and the fields Ω kj are calculated line-byline. The algorithm evolves along the time direction, taking time step k = 1 . . . N for a fixed space index j. This is what we refer to as a time line. After one time line is finished, the algorithm proceeds to the next space index. In the first time line, the values of 1k for k = 2 . . . N can be calculated from the initial condition 11
using only Eq. (15), since Ω1k for k = 1 . . . N is given by the boundary conditions of the fields. After the first time line is calculated, it is written to the output files. All following time lines are calculated in a loop over the space direction (see dark colored frame in Fig. 2). At time index k = 1, the density matrix 1j is given by the initial condition. Then, Ω 1j can be calculated from the last time line in a space step using Eq. (21). All following values of Ω kj are calculated in a second, nested loop (see light colored frame in Fig. 2). Since kj is necessary for the source function, it has to be calculated from the just gained value of Ω kj −1 and kj −1 using Eq. (15). In this manner, Eqs. (21) and (15) are used in alternation until the
743
last space step at k = N is performed by Eq. (22). After the calculation of a time line is finished, it is written to the output files if a storage interval is reached. Once the loop over the space direction is finished, the final time line is also written to an additional output file. 3.2. Memory structure The data structures for the density matrix and the fields each consist of a three-dimensional array. The first index represents the time step, the second index distinguishes between the last two space steps, and the third index runs over the elements of the density matrix or the field envelopes. Since the number of time steps depends on the values given in the input file, the actual size of the data structure is not known at compilation time such that it has to be allocated dynamically. Rather than allocating memory for each time step and both stored space steps separately, this is done in one block. This is advantageous because then, the data for one calculation step is close together in the physical memory which decreases memory access time. It requires, however, to set up the corresponding pointer structure by hand. The scheme is shown in Fig. 3. Each of the three layers of the data structure is allocated in one block and afterwards the connecting pointers are set. In the algorithm, further speed optimization is achieved by interchanging pointers rather than copying data once the new values for the elements of the density matrix and the field envelopes have been calculated. The new values are stored into the data structure addressed by the second index equal to 1 which overwrites the one before last values in space direction. Then, the pointers with second index 0 and 1 are interchanged, updating the data structure to the new values in the correct order for the next space step. 3.3. Input file The input file is a plain text file, used to set the physical and numerical parameters necessary for the numerical integration. Its name can be freely chosen and must be provided as the first command line argument when the msprop program is executed. The input file starts with an initial line describing the purpose of the file. Then, for each block of physical or numerical parameters a single comment line is followed by the actual data. First, the different decay rates, detunings, and field strength values are given. In the case of the system with a cross coupling of the electric and magnetic field component (see the third example, Section 4.3), there are additional values for the multi-photon phase, the helicity of the probe field, and the parameters of the dark state. Then, the peak and switch times, the respective width of pulsed fields, and the type of pulse envelope function are given. These parameters fix the initial conditions of the fields at the medium entry. It follows the value for L, necessary for the dimensionless coupling constant in the source function (see Section 2.2). Next, the parameters determining the numerical grid are specified, namely the maximum propagation time and distance as well as the numerical time and space step. The number of time and space steps are calculated in the program. But to reduce the amount of data, additional values for the number of time and space steps written to the output files are set. Finally, a frames-per-second value for the conversion to a movie file by the accompanying bash script is given. 3.4. Output files The output files store the calculated numerical data. They are plain text files where the name consists of a first part provided by the second command line argument and a second part of the form “_t = time.dat”. The program writes one output file for each
744
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
Fig. 2. (Color online.) Overview of the msprop program structure. The program flow is represented by the arrows and the boxes combine related processes in groups. The two larger frames represent the loop for the space direction nested with the loop for the time direction. Conditional statements are represented by diamond-shaped boxes.
Fig. 3. (Color online.) Memory structure for the density matrix and the fields. The structure consists of three arrays (light, medium and dark colored frame) which are connected by pointers (arrows). Data cells with one index contain a pointer to an array of pointers (R [i ]), with two indices a pointer to an array (R [i ][ j ]), and with three indices the actual data (R [i ][ j ][k]).
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
745
Fig. 4. (Color online.) First example, time line at the medium exit. The numerical probe pulse (solid red line) is compared to the analytic solution (dashed green line) for constant control field (dotted blue line). All fields are normalized to their initial value and the dimensionless time has been used.
storage time containing the numerical data for all storage points in space direction. These output files are written into a “data” subfolder. Another single output file with the name “timeline.dat” is written into the main folder. It contains the numerical data at the medium exit (last grid point in space direction) for all storage times. The structure of the output files is column based. The output files in the “data” subfolder contain the numerical value of the space grid point in the first column. Similarly the “timeline.dat” file contains the value of the time grid points in the first column. The other columns contain the values of the fields and the elements of the density matrix in the order of the corresponding data structures. Since this order depends on the name conventions of the respective quantum optical level scheme, a short specification of this data structures is given in the beginning of each source code file. 4. Examples With this manuscript, we provide three examples to illustrate the capabilities of the msprop program. The first example deals with a slow light pulse in an EIT system with constant control field. It can be solved analytically to a very good approximation and serves as a test case for the correct functioning of the program. The second example takes the same system to a different parameter range. Here, the numerical solution displays the full spatio-temporal dynamics of the light storage mechanism which in the general parameter regime is not accessible to analytical solutions. Finally, the third example deals with a system in which both the electric and magnetic components of the same probe field couple to the medium. It demonstrates how both components can be taken into account simultaneously in the numerical treatment. The numerical solution shows the specific phase-dependent dynamics which result from a cross-coupling of the electric and magnetic component. For all examples we show a plot of the time line at the medium exit, the numerical data from the last grid point in space direction and for all stored grid points in time direction. These plots can be generated directly from the numerical data with the gnuplot scrips “gp_timeline.sc”.
4.1. Slow light pulse The dynamics of a slow light pulse propagating in a three-level EIT medium can be solved analytically in the limit of a weak, adiabatic probe pulse [4,5]. We use the analytic formula for the group velocity to fix the numerically expected propagation velocity. This leads to the following expressions for the dimensionless propagation velocity u and the dimensionless coupling η [see Eqs. (13) and (11)]
u= 1+
3L γ p ω
Ωc2
−1 (23)
,
3
η = L ωu ,
(24)
2
where L = N λ3 /4π 2 is defined by the density N and the transition wavelength λ, and γ p , ω , Ωc are the probe field decay rate, probe field frequency, and the control field Rabi frequency. Here, and in the following, we express all numerical quantities in units of γ , the overall decay rate of the excited state. For a Gaussian probe pulse the dimensionless solution for the envelope function reads
Ω p ( z, t ) = Ω p
(t − z)2 σ , exp − σ˜ 2σ˜ 2
(25)
where Ω p and σ are the initial probe field amplitude and temporal width and σ˜ 2 = σ 2 + 4z/Ωc2 is the propagation dependent width of the pulse. To compare this solution to numerical data generated by the msprop program, we use a probe field amplitude of Ω p = 0.01, a control field amplitude of Ωc = 2, and a pulse length of σ = 50 to propagate the pulse over a length of z = 400. With these parameters the regime of a weak adiabatic probe field is well fulfilled. In Fig. 4 we show the data from the time line at the medium exit and compare it to the analytic solution. The exit time of the pulse is given by adding the time when the pulse enters the medium (t = 200) and the propagation time through the medium. Due to the scaling, the numerical velocity z/t = 1 and the propagation time through the medium (t = 400) is numerically
746
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
Fig. 5. (Color online.) Second example, time line at the medium exit. The probe field (solid red line) and the control field (dashed blue line) are shown. Both fields are scaled to their initial value and the dimensionless time has been used.
equivalent to the medium length (z = 400). Thus, the exit time of the pulse is given by t = 200 + 400 = 600 which is confirmed by the data shown in Fig. 4. The complete propagation dynamics can be studied with the “movie.avi” file which is generated by the accompanying bash script “run.bash”. 4.2. Storage of light The second example features the storage of a light pulse by switching the control field during the propagation time of a slow light pulse in an EIT medium. Since the EOM are the same as in the first example, the same expressions for u and η apply. To demonstrate the evolution of the spin coherence, we use a probe field strength of Ω p = 1.5 which is comparable to the control field strength of Ωc = 2 and a small ground state decoherence of γdec = 0.0025. The medium is now z = 150 long and the pulse with a width of σ = 20 enters the medium at t = 75. The control field is switched off at t = 150 and switched back on at t = 300 such that we have a storage time of t = 150. Since the exit time of the pulse is given by the sum of the entering time (t = 75), the propagation time (t = 150), and the storage time (t = 150), the retrieved pulse exits the medium at t = 375. In Fig. 5 we show the time line at the medium exit for the probe and control field. One can clearly see the signature of the energy transfer from the probe field to the control field at the moment where the pulse enters the medium (peak on top of control field at t = 75), the switch off time of the control field, and the adiabaton dip in the control field accompanying the pulse. The complete dynamics of the storage process can be studied with the movie file generated by the bash script. Here, in addition to the probe and control field, the build up and decay of the spin coherence are shown. 4.3. Phase controlled propagation in a system with a coupling of the electric and magnetic field component In the third example we deal with a system that exhibits a simultaneous coupling to the electric and magnetic component of a laser field. It serves to demonstrate the capabilities of the msprop
program to take into account both components and their coupling to the polarization and magnetization of the medium. In particular, it shows how to implement the modified Lax–Wendroff source function [Eq. (10)] into the numerical algorithm. A detailed discussion of the example system and its theoretical modeling are given in [12]. There, it was shown that the group velocity for a pulsed probe field is in principle given by an EIT-type expression and we find for the numerically expected propagation velocity,
u= 1+
3L γ34 ω44
Ωc2
−1 ,
(26)
where γ34 is the decay rate on the electric probe field transition and 44 is the population fraction in the lower state of the electric probe field transition. The scaling of all time units is done with the overall decay rate out of the upper state of the electric probe field transition. The physical parameters are closed to the example which is discussed in [12]. A Gaussian probe pulse of temporal width σ = 50 enters a z = 1200 long medium at t = 400. We assumed a leftcircular polarization of the probe field (helicity value +1) and an initial multi-photon phase of π /2. During the propagation time at t = 800 the multi-photon phase is switched to zero and at t = 1200 back to π /2. This changes the effective refractive index of the probe field. While in the beginning and at the end the probe field accumulates phase, it experiences gain during the intermediate time. In Fig. 6 we show the time line at the medium exit for the probe field and the multi-photon phase. One can clearly see the switching of the multi-photon phase as well as the change in phase and increase in amplitude of the probe pulse. The exit time (t = 1600) is given by the sum of the entering time and the propagation time (t = 1200). For the complete dynamics, we refer to the movie file provided with the example. 5. Summary We have implemented a numerical integration of the Maxwell– Schrödinger equations capable of solving the propagation dynamics of several laser fields through a medium of quantum optical fewlevel atoms in time and in one space dimension. In particular, our
R. Fleischhaker, J. Evers / Computer Physics Communications 182 (2011) 739–747
747
Fig. 6. (Color online.) Third example, time line at medium exit. The absolute value (thick solid black line), the real part (dashed blue line), and the imaginary part (solid red line) are shown. Their value is scaled to the initial value of the probe field on the left vertical axis and the dimensionless time has been used. The phase of the probe pulse (straight dashed green line) and the multi-photon phase (dash-dotted cyan line) are also shown and their value can be read off from the right vertical axis. For a clear view the probe field phase is set to zero at times where the absolute value of the probe field is below one percent of its peak value.
program is able to deal with systems where a simultaneous coupling of the electric and magnetic field component of the same laser field has to be taken into account. The program employs an Adams predictor method together with a Lax–Wendroff scheme to perform the numerical time and space step. For the Lax–Wendroff scheme the source function has been suitably adapted to take into account a coupling to the polarization and the magnetization of the medium. In addition to a detailed description of the algorithm and the program structure we have provided three examples which illustrate the capabilities of our program. Acknowledgements We would like to thank T.N. Dey for helpful discussions. References [1] Z. Ficek, S. Swain, Quantum Interference and Coherence: Theory and Experiments, Springer, Berlin, 2005.
[2] M.O. Scully, M.S. Zubairy, Quantum Optics, Cambridge University Press, Cambridge, 1997. [3] S.L. McCall, E.L. Hahn, Phys. Rev. Lett. 18 (1967) 908; S.L. McCall, E.L. Hahn, Phys. Rev. 183 (1969) 457. [4] S.E. Harris, Phys. Today 50 (1997) 36. [5] M. Fleischhauer, A. Imamoglu, J.P. Marangos, Rev. Mod. Phys. 77 (2005) 633. [6] R. Grobe, F.T. Hioe, J.H. Eberly, Phys. Rev. Lett. 73 (1994) 3183. [7] J.R. Csesznegi, R. Grobe, Phys. Rev. Lett. 79 (1997) 3162. [8] M.Ö. Oktel, Ö.E. Müstecaplıo˘glu, Phys. Rev. A 70 (2004) 053806. [9] Q. Thommen, P. Mandel, Phys. Rev. Lett. 96 (2006) 053601. [10] J.Q. Shen, S. He, Ann. Phys. 16 (2007) 364. [11] J. Kästel, M. Fleischhauer, S.F. Yelin, R.L. Walsworth, Phys. Rev. Lett. 99 (2007) 073602. [12] R. Fleischhaker, J. Evers, Phys. Rev. A 80 (2009) 063816. [13] S.J. Buckle, S.M. Barnett, P.L. Knight, M.A. Lauder, D.T. Pegg, Opt. Acta 33 (1986) 1129. [14] D.V. Kosachiov, B.G. Matisov, Y.V. Rozhdestvensky, J. Phys. B 25 (1992) 2473. [15] M. Mahmoudi, J. Evers, Phys. Rev. A 74 (2006) 063827. [16] R. Fleischhaker, J. Evers, Phys. Rev. A 77 (2008) 043805. [17] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes, Cambridge University Press, Cambridge, 1992. [18] K.R. Hansen, Standing wave electromagnetically induced transparency, Master thesis, University of Aarhus, 2006.