## Abstract

A compact numerical method for simulating ultrafast pulse interaction with inhomogeneously broadened multi-level media is reported. We use a low-dispersion pseudospectral scheme with fourth order time stepping for Maxwell’s equations, and a weakly coupled operator splitting method for the Bloch equations where inhomogeneous broadening and relaxations are also taken into account. The underlying physics is briefly discussed with emphasis on the formalism used.

©2011 Optical Society of America

## 1. Introduction

Today, laser technology makes it possible to generate electromagnetic pulses with durations of only a few cycles per pulse. It is, for example, possible to obtain 5fs pulses from titanium:sapphire lasers at 800nm by using double chirped mirrors for pulse compression. The interaction of the fundamental 800nm wavelength from a titanium:sapphire laser system and its second harmonic can also be used to excite third order nonlinearities in nitrogen plasmas, which acts as a source for single-cycle THz pulses with electric field amplitudes as high as 40MV/m [1]. As examples of applications of ultrashort pulse technology, femtosecond pulses are currently being used by spectroscopists in the characterization of the electronic properties of matter, sub-cycle THz pulses have found useage in the manipulation of atomic wavefunctions in Rydberg atoms [2], and vibrational cooling of molecules has recently been achieved by optical pumping with broadband pulses [3]. The coherent matter-field interaction lies at the heart of many of these applications. However, for short and/or strong pulses many of the standard approximations start to break down. The slowly varying envelope approximation(SVEA), for example, requires the pulse’s bandwidth to be small compared to its central frequency. For femtosecond pulses at optical wavelengths there are only a few cycles per pulse and the SVEA is not very well satisfied. For the single- and sub-cycle regimes the SVEA has lost all validity. Also, a two- or three-level approximation is often made for the medium, on the grounds that the pulse is spectrally narrow and effectively interacts with only one or two transitions, and all other transitions are spectrally far removed. In practice, the energy levels of atoms and molecules are more evenly spaced and a few-level approximation becomes questionable when the bandwidth of the pulse increases. An example of such a situation is the propagation of broadband THz pulses in the vibration-rotation bands of molecular gases. Some works, related to ultrashort pulse propagation in multilevel systems within SVEA, can be found in e.g. [4, 5] and references therein. Another standard approximation is the rotating wave approximation(RWA), which requires near-resonant pulses and Rabi periods long compared to an optical cycle. The RWA is a basic matter-field approximation, and it also breaks down in the short and strong pulse regimes due to the broad spectrum of the pulse, as well as the possible generation of high-frequency content when Rabi flops take place on the timescale of an optical cycle [6]. With the breakdown of these approximations it is very difficult to develop analytical tools that are appropriate for investigating pulse propagation effects, and one must therefore often look to approximate numerical methods.

Traditionally, simulations of nonlinear optical interactions deal with phenomenological models that are based on macroscopic properties of the medium, such as split-step finite difference or spectral methods based on Debye or Drude dispersion models. However, a phenomenological approach fails spectacularly in the strong and short pulse regimes when coherent quantum effects first begin to manifest. This necessitates simulation models that take into account the inherent microscopic structure of matter. An approach that has proven successful is to model the system by use of slowly varying envelopes for the electric field and the density matrix elements, which is often supplemented by a rotating wave approximation for the matter-field interaction. However, such models hold no validity for few-cycle pulses. The first numerical simulations of the full two-level Maxwell-Bloch equation set was initiated by Ziolkowski in [7] for two-level atoms. The approach in [7] is based on the Yee scheme [8] for Maxwell’s equations together with a Crank-Nicholson scheme for the Bloch equations, with corrections later given in [9]. Since then, the Ziolkowski approach has been used to predict attosecond radiation from femtosecond carrier-wave Rabi flopping [10], the formation of optical subcycle pulses in dense media [11] and breakup and self-focusing of 2*π* pulses [12]. The method has also been extended to inhomogeneously broadened media [13]. Indeed, the Ziolkowski approach works quite well for two-level systems but does not generalize to more levels as it would lose the positiveness property for the density matrix [14]. In [14] and [15], a generally valid method that decouples the Maxwell and Bloch equations is introduced. The results in [14] and [15] are quite good, and this method has since been extended to propagation in anisotropic media [16].

In this paper, we extend the work in [14] by working in a Liouville [17] space. In particular we focus on the implementation of an arbitrary number of levels and show how inhomogeneous broadening can be included. Our simulations are based on weakly coupled methods for the Bloch equations together with pseudospectral time domain methods [18, 19] for Maxwell’s equations. Although our model is formulated generally, we focus on one-dimensional simulations only.

## 2. Theoretical model

The Hamiltonian *Ĥ* = *Ĥ*(** r**,

*t*,

**,**

*v**θ*,

*ϕ*) for an absorber interacting with a classical electromagnetic field

**=**

*E***(**

*E***,**

*r**t*) under the dipole approximation and expanded in the energy basis is

*i*and

*j*over the

*N*available energy levels.

*ħω*is the energy of the energy eigenstate |

_{i}*i*〉, and

*μ*_{ij}= 〈

*i*|

**|**

*$\widehat{\mathit{\mu}}$**j*〉 is the

*ij*element of the effective dipole moment operator

**. The spherical angles**

*$\widehat{\mathit{\mu}}$**θ*,

*ϕ*and the velocity

**are used to classify absorbers that have a characteristic orientation in the**

*v**θ*,

*ϕ*direction and a spatial velocity

**.**

*v**Ĥ*is implicitly understood to belong to one of these classes in the unit volume at (

**,**

*r**t*). It is natural to assume that

**and**

*$\widehat{\mathit{\mu}}$**ω*are not functions of

_{i}**and**

*r**t*which means that the medium is isotropic at all times, although the extension to spatially varying media is straightforward.

For notational simplicity, we introduce a notation that denotes the fractional sum over all the velocities and orientations as

*n*(

**,**

*v**θ*,

*ϕ*)

*d*

*v**d*Ω is the fraction of absorbers with orientation (

*θ*,

*ϕ*) within the solid angle

*d*Ω, and a spatial velocity in the interval [

**,**

*v***+**

*v**d*

**]. The total number density of the medium is denoted by**

*v**N*.

_{a}#### 2.1. Liouville space

From a numerical point of view it is convenient to modify the von Neumann equation by working with a vectorized form of the density operator *$\widehat{\rho}$*. Therefore, we seek to replace the von Neumann equation *d _{t}$\widehat{\rho}$* =(

*iħ*)

^{−1}[

*Ĥ*,

*$\widehat{\rho}$*] by the vectorized version

*N*, the density operator is typically represented as an

*N*×

*N*self-adjoint, positive definite matrix with a unity trace(if the system is closed). Here, we have replaced that matrix by a column vector |

*ρ*〉〉 of length

*N*

^{2}containing all the entries of

*$\widehat{\rho}$*, and replaced the commutator [

*Ĥ*,] of dimension

*N*×

*N*by another operator

*Ľ*of dimension

*N*

^{2}×

*N*

^{2}. To stay consistent with spectroscopic literature [17], we shall call the space in which |

*ρ*〉〉 exists for a Liouville space. We shall denote Liouville space vectors by double kets and Liouville space operators by inverted hats.

For every operator *Ô* in a Hilbert space, there exists an equivalent vector |*O*〉〉 in a Liouville space. These are generally written as

*jk*is implicit.

*O*is the

_{jk}*jk*element of the Hilbert space operator

*Ô*and the

*jk*entry in the Liouville vector |

*O*〉〉. We define a bra-vector 〈〈

*O*| that corresponds to

*Ô*

^{†}and a scalar product of two vectors by

*jk*|

*mn*〉〉 =

*δ*, and a completeness relation Σ

_{jm}δ_{kn}*|*

_{jk}*jk*〉〉〈〈

*jk*| = 1. Every Liouville space operator can therefore be expanded as

*Ǒ*= Σ

_{jk,mn}*O*|

_{jk,mn}*jk*〉〉〈〈

*mn*|.

The order of the entries in the column vector |*ρ*〉〉 is not fundamentally important, and the operator *Ľ* can be formed in a general fashion once the order is chosen (the order of every other operator follows the same order as |*ρ*〉〉). However, we observe that the matrix equation *C* = *AXB* can be written as |*C*〉〉= (*B*
^{⊤} ⊗ *A*)|*X*〉〉 if *A* and *B* are square matrices and |*C*〉〉 and |*X*〉〉 are column major ordered (a similar expression exists for row major ordering). *Ľ* can therefore be written in the form

*A*⊕

*B*=

*A*⊗

*Î*

_{dim(}

_{B}_{)}+

*Î*

_{dim(}

_{A}_{)}⊗

*B*, where ⊗ is the Kronecker product.

#### 2.2. Relaxations and polarization

Physical effects not accounted for in the Hamiltonian in Eq. (1) can included phenomenologically by the addition of a relaxation operator *Ř* on the right hand side of Eq. (3),

*Ř*is that the equation

*d*|

_{t}*ρ*〉〉 =

*Ř*|

*ρ*〉〉 must preserve the physical properties of

*$\widehat{\rho}$*. It is convenient to derive

*Ř*directly from the master equation

*$\widehat{\gamma}$*are of the form ${\widehat{\gamma}}_{ij}\hspace{0.17em}=\hspace{0.17em}\sqrt{{\gamma}_{ij}}|i\u3009\u3008j|$ with

_{ij}*γ*real. Equation (8) preserves positiveness, hermiticity and trace. Generally, operators

_{ij}*$\widehat{\gamma}$*with

_{ij}*i*≠

*j*represent damping that cause population transfer from state

*j*to state

*i*at a rate

*γ*(for example due to inelastic collisions, spontaneous emission or a thermal background). Operators

_{ij}*γ*

*̂*cause decoherence only (for example due to elastic collisions). The off-diagonal and diagonal terms decay according to

_{ii}*γ*’s) as well as amplitude decay rates (

_{ii}*γ*’s). The effective decay rates for the populations

_{ij}*ρ*consist of population transfer from all states to

_{mm}*ρ*minus the population transfer to all states from

_{mm}*ρ*. The relaxation operator can be written in Liouville space as

_{mm}*γ*’s are real.

_{ij}In the literature, the initial condition for the medium is often taken as a coherent ensemble, or the ground state if the transitions are highly energetic or the medium is cold(all factors Δ*E*/*k _{B}T* being very large). Other times, the initial state of the medium must be taken as an incoherent ensemble with thermally distributed populations. An example of such a situation is the propagation of strong terahertz and microwave pulses in the vibrational and rotational levels in molecular gases. Relaxation to a thermal equilibrium is forced by taking

*γ*=

_{mj}*γ*exp[−

_{jm}*β*(

*ω*−

_{m}*ω*)] where

_{j}*β*=

*ħ*/

*k*. Independent of forcing the thermal restriction, the equilibrium state is the solution to the rate balance equation Σ

_{B}T

_{j}_{≠}

*− Σ*

_{m}γ_{mj}ρ_{jj}

_{j}_{≠}

_{m}*γ*= 0, and all off-diagonal elements equal to zero. This can be written as the linear system Γ|diag(

_{jm}ρ_{mm}*$\widehat{\rho}$*)〉= 0 where |diag(

*$\widehat{\rho}$*)〉 is the column vector corresponding to the diagonal elements of

*$\widehat{\rho}$*and Γ is the matrix corresponding to the linear system Σ

_{j}_{≠}

*− Σ*

_{m}γ_{mj}ρ_{jj}

_{j}_{≠}

_{m}*γ*= 0. The solution is known [20], and the populations are ${\rho}_{mm}\hspace{0.17em}=\hspace{0.17em}\frac{{\mathrm{\Gamma}}^{jm}}{{\mathrm{\Sigma}}_{j}{\mathrm{\Gamma}}^{jm}}$ where Γ

_{jm}ρ_{mm}*is the cofactor of the element Γ*

^{jm}*. Γ has the special shape that all cofactors in every column are identical (Γ*

_{jm}^{1m}= Γ

^{2m}= ... etc.), so the cofactor Γ

*is independent of*

^{jm}*j*.

With *Ř* now included implicitly in *Ľ*, the polarization ** P** = 〈Tr(

*$\widehat{\mathit{\mu}}$**$\widehat{\rho}$*)〉

_{Σ}and the time-derivative of the polarization,

*d*

_{t}**, can be written by the use of Eq. (3), Eq. (5) and the hermiticity of**

*P***as**

*$\widehat{\mathit{\mu}}$*#### 2.3. Time evolution

By splitting the Liouville operator into time dependent and time independent parts, *Ľ* = *Ľ*
_{0} + *Ľ*
_{1}(*t*), the exact solution to Eq. (3) is

*Ǔ*

_{0}(

*t*,

*t*

_{0}) is

*Ǔ*

_{0}(

*t*,

*t*

_{0}) = exp[

*Ľ*

_{0}(

*t*−

*t*

_{0})]. ${\mathrm{\Sigma}}_{n=1}^{\infty}\stackrel{\u02c7}{\mathrm{\Omega}}(t,\hspace{0.17em}{t}_{0})$ is the Magnus series [21] where the first term reads

*Ľ*(

_{I}*τ*

_{2}),

*Ľ*(

_{I}*τ*

_{1})], is negligible over the time

*t*–

*t*

_{0}= Δ

*t*. The transition in Eq. (7) does not preserve the anti-hermiticity in Eq. (6), so the free propagator

*Ǔ*

_{0}is not necessarily unitary, although it still preserves the trace, hermiticity and positiveness of

*$\widehat{\rho}$*.

#### 2.4. Inhomogeneous broadening

For non-zero temperatures the thermal motion of absorbers will give rise to collision broadening and Doppler shifts on optical transitions. Consequently, instead of sharp-line absorption, there will be a spectrum of absorption frequencies centered around every resonance. In one-dimensional unidirectional propagation, Doppler broadening can be incorporated into our model by letting *ω _{i}* =

*ω*

_{i}_{0}(1 +

*v*/

*c*) where

*ω*

_{i}_{0}is the energy divided by

*ħ*of the level

*i*for an absorber at rest, and

*v*is its velocity. Since Doppler broadening is a statistical time-independent phenomena in our model, we shall see how it, in practice, only affects

*Ǔ*

_{0}, which leads to a very efficient numerical evaluation.

The line profile that is most commonly used in spectroscopic literature is the Voigt profile, which results from convolving two line broadening mechanisms, one that gives a Lorentzian profile(e.g. collision broadening) and one that gives a Gaussian profile(e.g. Doppler broadening). For gases it is usually a good approximation to only use a Gaussian profile due to the broad Doppler distribution, and later we give a numerical example of a photon echo simulation where Doppler broadening plays a key role.

#### 2.5. Maxwell’s equations and dimensionless units

To ensure maximal precision of numerical simulations, we introduce characteristic scales such that ** E** is of order one, and a time-scale

*t*= Ω

_{c}^{−1}= (

*E*/

_{c}**μ**_{c}*ħ*)

^{−1}which is equal to the inverse Rabi frequency Ω, a quantity that is often on the order of an optical transition for very strong matter-field interactions. Scaling of the rest of the model follows from these quantities: the magnetic field

**is scaled according to $\mathit{\text{H}}\hspace{0.17em}\to \hspace{0.17em}({E}_{c}\sqrt{\varepsilon /\mu})\mathit{\text{H}}$, $\mathit{\text{r}}\hspace{0.17em}\to \hspace{0.17em}({t}_{c}/\sqrt{\mu \varepsilon})\mathit{\text{r}}$,**

*H**Ľ*→

*Ľ*/

*t*and |

_{c}**〉〉 →**

*μ**μ*|

_{c}**〉〉. The density operator is already of order one and dimensionless, and therefore remains untouched. The entire Maxwell-Bloch model can now be summarized by the dimensionless equations**

*μ*## 3. Numerical model

#### 3.1. Method for Bloch equations

A generally valid method for the Bloch equations is what is known as an *operator splitting method*, and it was first introduced in [14] for the Maxwell-Bloch system. Here we show how the splitting method can be recovered from the exact solution |*ρ*(*t*)〉〉= *Ǔ*
_{0}(*t*, *t*
_{0})*Ǔ _{I}*(

*t*,

*t*

_{0})|

*ρ*(

*t*

_{0})〉〉, where

*Ǔ*

_{0}(

*t*,

*t*

_{0}) = exp[

*Ľ*

_{0}(

*t*–

*t*

_{0})] and ${\stackrel{\u02c7}{U}}_{I}(t,\hspace{0.17em}{t}_{0})\hspace{0.17em}=\hspace{0.17em}\text{exp}[{\mathrm{\Sigma}}_{n=1}^{\infty}{\stackrel{\u02c7}{\mathrm{\Omega}}}_{n}(t,\hspace{0.17em}{t}_{0})]$. Evaluating this on

*t*

_{0}= (

*n*− 1/2)Δ

*t*,

*t*= (

*n*+ 1/2)Δ

*t*and keeping only the first order term in the Magnus series gives

*$\widehat{\rho}$*as long as the matrix exponentials can be calculated exactly. Note that |

*ρ*〉〉 and

**are stored on staggered temporal grids, which is the general trend throughout this paper.**

*E*There are several approximations and simplifications that can be made in Eq. (16)
. Firstly, the computation of the matrix exponential
$\text{exp}({\stackrel{\u02c7}{L}}_{I}^{n}\mathrm{\Delta}t)$ must be performed at every time step and therefore has the potential of being the bottleneck in a numerical simulation. It may therefore be tempting to approximate *Ǔ _{I}*, but in doing so one must be careful because truncating the exponential has the potential to ultimately lead to violation of the trace, hermiticity and positiveness property of

*ρ*

*̂*. On the other hand,

*Ǔ*

_{0}is time independent and only needs to be calculated once for each class of absorbers.

Secondly, we may take
${\stackrel{\u02c7}{L}}_{I}^{n}\hspace{0.17em}\approx \hspace{0.17em}{\mathit{\text{E}}}^{n}\hspace{0.17em}\cdot \hspace{0.17em}{\stackrel{\u02c7}{\mathit{\text{L}}}}_{1}$. This reduces the number of matrix multiplications that are necessary at each time step and also has the potential of increasing the rate of convergence when splitting the operators. This approximation is also more consistent with the previous approximation of keeping only the first order term in the Magnus series since *Ľ _{I}* =

*Ľ*

_{1}+

*𝒪*(Δ

*t*). If inhomogeneous broadening is included, taking ${\stackrel{\u02c7}{L}}_{I}^{n}\hspace{0.17em}\approx \hspace{0.17em}{\mathit{\text{E}}}^{n}\hspace{0.17em}\cdot \hspace{0.17em}{\stackrel{\u02c7}{\mathit{\text{L}}}}_{1}$ means that only one interaction propagator needs to be calculated at each spatial point at each time step. This drastically reduces the necessary computations with introducing negligible errors only.

Next, with these approximations the interaction propagator ${\stackrel{\u02c7}{U}}_{I}^{n}$ can be written

**⊕**

*A***) = exp(**

*B***) ⊗ exp(**

*A***) and the hermiticity of**

*B***. The point is that instead of calculating the matrix exponential exp (–**

*$\widehat{\mathit{\mu}}$**i*

*E**· [*

^{n}

*$\widehat{\mathit{\mu}}$*^{⊤}⊕ (−

**)]Δ**

*$\widehat{\mathit{\mu}}$**t*) where the argument is of order

*N*

^{4}, we have reduced it to calculating the exponential exp(

*i*

*E**·*

^{n}**) where the argument is of order**

*$\widehat{\mathit{\mu}}$**N*

^{2}. The approximations made so far essentially leads to the method introduced in [14], and the solution is summarized as

Additionally, higher order approximations can be achieved by rearranging the propagators in Eq. (18),

*ρ*

^{n+1/2}〉〉 and converges faster than Eq. (18). We shall refer to rearrangements of the operators (such as Eq. (19)) as

*the*splitting method.

Finally, the propagators are very sparse and a large computational gain can be obtained by storing the propagators in sparse formats and using optimized matrix multiplication. Besides, for inhomogeneously broadened media the update of |*ρ*〉〉 rather than the calculation of the interaction propagator is the computational bottleneck, and the simulation time is greatly reduced by using sparse matrix multiplication. As an example, in our 6-level simulation only 66 of the 1296 elements in each of the 301 free propagators were non-zero. The evolution of |*ρ*〉〉 is also local, which leads to a very efficient parallelization.

#### 3.2. Pseudospectral time domain

The use of the Yee scheme comes with two large drawbacks that make it inefficient for modeling coherent pulse propagation. Firstly, the Yee scheme requires a very fine sampling in space and time in order to accurately model wave propagation. Since the largest computational burden in our model is the updating of the Bloch equations, reducing the grid resolution can lead to drastic improvements in efficiency. Secondly, the Yee grid leads to a large and complicated linear system that needs to be solved at every time step because ** E** enters on both sides of Eq. (14b). The problem is rooted in the staggeredness of the Yee grid, but does not enter in the one-dimensional model because 〈〈

*μ*|

*EĽ*

_{1}|

*ρ*〉〉 =

*E*Tr(

*$\widehat{\mu}$*, [

*$\widehat{\mu}$*,

*$\widehat{\rho}$*]) = 0. Both drawbacks of the Yee scheme can be remedied by the use of pseudospectral methods [18] for Maxwell’s equations. In this case, one can use a non-staggered spatial grid for

**,**

*E***and |**

*H**ρ*〉〉, but a temporally staggered grid is still necessary in order to achieve at least second order accuracy in time.

We define the Fourier derivative operator ${\partial}_{d}^{\mathcal{F}}$ by

*k*is the spatial frequency in the

_{d}*d*-direction,

*F*is a component of either

**or**

*E***, and**

*H**ℱ*and ${\mathcal{F}}_{d}^{-1}$ denote the Fourier and inverse Fourier transforms in the

_{d}*d*-direction. Taking the numerical approximation

*∂*≈ (𝔻

_{t}F*)*

^{t}F*, Maxwell’s equations can be summarized in their discrete forms as*

^{n}*μ*〉〉 =

_{x}**·|**

*x̂***〉〉, |**

*μ**μ*〉〉 =

_{y}**·|**

*ŷ***〉〉, |**

*μ**μ*〉〉 =

_{z}**·|**

*ẑ***〉〉. (𝔻**

*μ**)*

^{t}F*is traditionally taken as the second order Yee-type leapfrogging ${({\text{\mathbb{D}}}^{t}F)}^{n}\hspace{0.17em}=\hspace{0.17em}\frac{{F}^{n+1/2}-{F}^{n-1/2}}{\mathrm{\Delta}t}$. However, while the use of Fourier transforms means that the spatial derivatives are calculated exactly as long as the sampling is above the Nyquist limit, the temporal Yee-type leapfrogging scheme is a source for numerical dispersion. This can lead to artificial chirps in the pulse, and the spatial resolution must therefore be taken significantly higher than the Nyquist limit in order to achieve sufficient accuracy. We seek to partially overcome the limitations of the second order pseudospectral scheme(PSTD-2) by using a 4th order time stepping for (𝔻*

^{n}*)*

^{t}F*, which will allow us to push closer to the Nyquist limit. In one dimension,*

^{n}**=**

*E**E*,

_{x}**=**

*H**H*,

_{y}*∂*=

_{x}*∂*= 0, the fourth order numerical approximation for

_{y}*∂*

_{t}*F*can be written ${({\text{\mathbb{D}}}^{t}F)}^{n}\hspace{0.17em}=\hspace{0.17em}\frac{{F}^{n+1/2}-{F}^{n-1/2}}{\mathrm{\Delta}t}-\frac{\mathrm{\Delta}{t}^{2}}{24}{\mathcal{F}}^{-1}\hspace{0.17em}{[{(2\pi ik)}^{3}\hspace{0.17em}\mathcal{F}(F)]}^{n}$ [19]. We have neglected the cross terms

*ζ*Δ

*t*

^{3}because the dominant contribution from the polarization is already captured by the terms on the right hand side of Eq. (22) . The equivalent expressions in two and three dimensions can be found in [19].

With the same grid resolutions, the fourth order scheme is not significantly slower than a second order scheme because the largest computational burden lies in the update of |*ρ*〉〉 and the polarization terms on the right hand side of Eq. (22)
, and not the calculation of the Fourier transforms. Since the fourth order PSTD scheme(PSTD-4) allows a rougher grid resolution than PSTD-2 to achieve the same accuracy, the PSTD-4 method generally performs faster than an equivalent PSTD-2 method.

The PSTD method also has a more stringent condition on the time step than the Yee scheme does, and must be chosen according to Δ*t* ≤ 2min(Δ*x*, Δ*y*, Δ*z*)/(
$\sqrt{D}\pi $) (in dimensionless units), where *D* is the dimensionality. It is always necessary to choose spatial and temporal time steps small enough to resolve the generated spectrum, and not just the initial spectrum. The periodicity introduced through the Fourier transforms also means that absorbing layers must be added to the boundaries of the computational domain to counter wraparound effects.

As a final comment, the linear system that enters in Eq. (22) is block-diagonal and can be inverted by direct methods.

## 4. Numerical examples

We now show several scenarios that can be modeled by full-field Maxwell-Bloch equations, although our numbers do not necessarily describe a particular physical system. While we solve our equations in dimensionless form, all simulations are presented in SI units only. In all our simulations we take |** μ**〉〉 = |

*μ*〉〉 and use the PSTD-4 method for Maxwell’s equations and Eq. (19) for updating the Bloch equations. Our propagators are all calculated in full using the expokit [22] software library, and our computational domains are padded with absorbing layers with reflection coefficients less than 10

^{−6}.

#### 4.1. Two-level systems

Our first numerical simulation shows how an energetic electromagnetic pulse can travel loss-lessly through a two-level sharp line absorber. This pulse, the 2*π* pulse, owes its transparency to its ability to locally excite the medium and then coherently return the medium to its initial state in one inversion cycle. Likewise, a 2*πn* pulse will return the material to its initial state in *n* inversion cycles, a (2*n* + 1)*π* pulse inverts the medium through *n* inversions and a (2*n* + 1)*π*/2 pulse leaves the material in a state with *ρ*
_{11} = *ρ*
_{22} = 0.5 after inverting the medium *n* times. The material that we choose here is taken as *N _{a}* = 10

^{24}m

^{−3},

*μ*= 10

_{c}^{−29}Cm and

*ω*

_{0}= 4

*π*· 10

^{14}rad/s, which reflects the choice of parameters made in [7].

The injected pulse is *E*(*t*) = *ℰ*(*t*)sin(*ω*
_{0}
*t*) = *E _{c}*sech [(

*t*– 6

*τ*)/

*τ*]sin(

*ω*

_{0}

*t*) where

*E*= 4.2186·10

_{c}^{9}V/m and

*τ*= 5fs. This pulse has a FWHM ≈ 13fs. The length of our medium is

*L*= 37.5

**m with spatial and temporal time steps $\mathrm{\Delta}z\hspace{0.17em}=\hspace{0.17em}L/250\hspace{0.17em}=\hspace{0.17em}\frac{\lambda}{10}\hspace{0.17em}=\hspace{0.17em}0.15\hspace{0.17em}\mu m$ and $\mathrm{\Delta}t\hspace{0.17em}=\hspace{0.17em}\frac{2\mathrm{\Delta}z}{c\pi}\hspace{0.17em}\approx \hspace{0.17em}0.318\text{fs}$. The chosen values for**

*μ**E*,

_{c}*μ*and

_{c}*τ*reflect that the pulse area is $A\hspace{0.17em}=\hspace{0.17em}\frac{{\mu}_{c}}{\hslash}\hspace{0.17em}{\int}_{-\infty}^{\infty}\epsilon (t\prime )dt\prime \hspace{0.17em}=\hspace{0.17em}2\pi $, which is the necessary condition for the lossless propagation of the 2

*π*hyperbolic secant soliton shown in Fig. 1.

Secondly, we show how a photon echo [23] can form locally in an inhomogeneously broadened medium. Coupling to a Maxwell code is straightforward but does not illustrate the role of Doppler broadening further. We use the same parameters as in the previous example but include a Doppler distribution function

*ω*=

*ω*(1 +

_{o}*v/c*) =

*ω*

_{0}+ Δ

*ω*into Eq. (23) and integrating the free induction decay $P(t)\hspace{0.17em}=\hspace{0.17em}{\int}_{-\infty}^{\infty}d(\mathrm{\Delta}\omega )g(\mathrm{\Delta}\omega )\hspace{0.17em}\text{cos}[({\omega}_{0}\hspace{0.17em}+\hspace{0.17em}\mathrm{\Delta}\omega )t]$ we can extract a Doppler lifetime ${T}_{2}^{*}\hspace{0.17em}=\hspace{0.17em}\frac{2c}{{v}_{p}}\frac{1}{{\omega}_{0}}$. In order to demonstrate and validate the role of the mean free velocity

*v*, which acts as a control parameter of the inhomogeneous bandwidth, we choose to simulate three materials with

_{p}*v*/

_{p}*c*= 1/(50

*π*),

*v*/

_{p}*c*= 1/(100

*π*) and

*v*/

_{p}*c*= 1/(200

*π*). Note that these velocities may be unrealistically high for some materials, but they nevertheless serve to demonstrate the validity of our method. The theoretically derived Doppler lifetimes are ${T}_{2}^{*}\hspace{0.17em}=\hspace{0.17em}250\text{fs}$, ${T}_{2}^{*}\hspace{0.17em}=\hspace{0.17em}500\text{fs}$ and ${T}_{2}^{*}\hspace{0.17em}=\hspace{0.17em}1\text{ps}$. The inhomogeneous bandwidth ( $1/{T}_{2}^{*}$) lies in the THz range.

Numerically, we represent the broadening function by sampling *n _{l}* = 1001 absorption lines in the velocity interval

*v*∈

*c*[−0.025, 0.025] distributed according to Eq. (23). A large number of absorption lines is usually necessary in order to avoid artificial echoes due to finite sampling. The photon echo shown in Fig. 2 is generated by the injection of a

*π*/2 pulse that lifts the absorbers onto the equator plane of the Bloch sphere [23]. When the pulse has passed, the dipoles dephase due to their different resonance frequencies which causes a free induction decay of the macroscopic polarization. This decay is not irreversible since the microscopic polarization of each absorber is unaffected. An echo can therefore be reproduced from the dipoles by injecting a

*π*-pulse a time

*T*/2 after the first pulse has passed, which rotates the dipoles

*π*radians around the Bloch sphere. When again in the equator plane the dipoles continue their rotation around the Bloch sphere but the phase-shift of

*π*causes a rephasing after an additional time

*T*/2, which results in an echo of the original polarization. Our injected pulses are $E(t)\hspace{0.17em}=\hspace{0.17em}\frac{{E}_{c}}{4}\hspace{0.17em}(\text{sech}[\text{(}t\hspace{0.17em}-\hspace{0.17em}50\tau )/\tau ]\hspace{0.17em}+\hspace{0.17em}\text{sech}[(t\hspace{0.17em}-\hspace{0.17em}800\tau )/(2\pi )])\hspace{0.17em}\text{sin}({\omega}_{0}t)$, so the

*π*/2 and

*π*pulses have the same peak amplitude, but different pulselengths. We have verified, although not shown here explicitly, that the numerically calculated ${T}_{2}^{*}$’s agree very well with the theoretically predicted values of 250fs, 500fs and 1ps, which serves to validate our method.

#### 4.2. Multilevel systems

As a final example we show the absorption of a strong few-cycle pulse in a 6-level anharmonic ladder. We take the transitions of the medium at the spectral lines *ω _{n}*

_{+1}–

*ω*=

_{n}*ω*

_{0}[1 – 0.1(

*n*–3)], where

*n*= 1,2,3,4,5, and

*ω*

_{0}is the angular carrier frequency

*ω*

_{0}= 2

*π*· 10

^{13}rad/s. The injected pulse is

*E*(

*t*) =

*E*exp[−(

_{c}*t*–

*t*

_{0})

^{2}/

*τ*

^{2}]sin(

*ω*

_{0}

*t*) where

*τ*= 100fs,

*t*

_{0}= 3

*τ*and

*E*= 5 · 10

_{c}^{8}V/m. This pulse has a spectral content that covers all the resonances

*ω*

_{n}_{+1}–

*ω*, and the peak of the spectrum is located at the

_{n}*ω*

_{4}–

*ω*

_{3}resonance. We include the Doppler broadening function in Eq. (23) sampled at 301 Doppler lines with

*v*/

_{p}*c*≈ 10

^{−4}, which indicates an inhomogeneous bandwidth in the GHz range. Note that the inhomogeneous bandwidth is much smaller than spectral width of the pulse, and that the presence of inhomogeneous broadening therefore gives rise to minor modifications only. We take relaxation rates ${\gamma}_{12}^{-1}\hspace{0.17em}=\hspace{0.17em}1\text{ps}$, ${\gamma}_{23}^{-1}\hspace{0.17em}=\hspace{0.17em}1.1\text{ps}$, ${\gamma}_{34}^{-1}\hspace{0.17em}=\hspace{0.17em}1.2\text{ps}$, ${\gamma}_{45}^{-1}\hspace{0.17em}=\hspace{0.17em}1.3\text{ps}$, ${\gamma}_{56}^{-1}\hspace{0.17em}=\hspace{0.17em}1.4\text{ps}$ and ${\gamma}_{11}^{-1}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}{\gamma}_{22}^{-1}\hspace{0.17em}=\hspace{0.17em}\dots \hspace{0.17em}=\hspace{0.17em}{\gamma}_{66}^{-1}\hspace{0.17em}=\hspace{0.17em}1\text{ps}$. We force the thermal restriction

*γ*=

_{mj}*γ*exp[−

_{jm}*β*(

*ω*–

_{m}*ω*)] with

_{j}*T*= 600K, which means that the initial state of our medium is

*ρ*

_{11}≈ 0.6,

*ρ*

_{22}≈ 0.23,

*ρ*

_{33}≈ 0.096,

*ρ*

_{44}≈ 0.044,

*ρ*

_{55}≈ 0.02,

*ρ*

_{66}≈ 0.01 and all coherences equal to zero. We take the density as

*N*= 10

_{a}^{25}m

^{−3}and take

*μ*

_{12}=

*μ*

_{23}= ...

*μ*

_{56}= 10

^{−29}Cm. The length of our simulation region is

*L*= 1mm with spatial and temporal time steps Δ

*z*=

*L*/500 = 2

**m and Δ**

*μ**t*= 2Δ

*z*/(

*cπ*) ≈ 4.2fs.

In the following, we make no attempt to explain the physics but rather comment on the general behaviour of the system. Our pulse is strong enough to appreciably affect the populations since the normalized pulse area (the average Rabi frequency)
$\overline{\mathrm{\Omega}}\hspace{0.17em}=\hspace{0.17em}A/\tau \hspace{0.17em}=\hspace{0.17em}({E}_{c}{\mu}_{c}\sqrt{\pi )}/\hslash $ is higher than all decoherence rates. Consequently, we expect to observe population fluctuations during the passage of the pulse. Figure 3 shows the time and frequency domain propagation of the pulse through the medium. The pulse is quickly absorbed at the entrance of the medium and it is easy to observe flattenings of the spectrum on every resonance. The most conspicuous absorption is found at the *ω*
_{2} – *ω*
_{1} and *ω*
_{3} – *ω*
_{2} resonances where the population is highest. The pulse propagates with lower absorption after *z* ≈ 0.5mm due to the spectral hole at the *ω*
_{3} – *ω*
_{2} and *ω*
_{2} – *ω*
_{1} resonances.

Figure 4 shows the normalized energy (〈*Ĥ*
_{0}〉_{Σ} – *ω*
_{10})/(*ω*
_{60} – *ω*
_{10}), the electric field and averaged populations 〈*ρ*
_{11}〉_{Σ}, 〈*ρ*
_{22}〉_{Σ},..., 〈*ρ*
_{66}〉_{Σ} at the entrance, *z* = 0, of the medium. It is quite obvious that the strong pulse excites the medium, and that the medium relaxes toward its thermal equilibrium after the pulse has passed. During the passage, the medium is excited to a high energy level, and the population of the highest energy level, 〈*ρ*
_{66}〉_{Σ}, reaches 0.33 during propagation. This final example also illustrates the simultaneous breakdown of the two- or three-level approximation, the SVEA and RWA.

## 5. Conclusion

The usefulness of the method presented here lies in the compactness of the formalism. Indeed, the development of an *N*-level code should be of the same complexity as 2-level code, apart from the additional degrees of freedom of the medium (it will also require longer simulation times). The use of a Liouville space formalism together with pseudospectral methods not only allows gain in computational time, but also allows us to reduce the number of polarization terms in Maxwell’s equations which leads to simpler formulations of multi-dimensional systems. We include homogeneous relaxations via the master equation, which ensures that the positive-definite property of the density matrix is not violated, and we also include inhomogeneous broadening via the free propagators. The most expensive computational part for inhomogeneously broadened systems is the update of the density operators and the polarization terms, which is made more efficient by use of sparse matrix multiplication and storage. The method can be made even more efficient by use of numerical approximations and simplifications for the interaction propagator, parallelization of the code and more optimized algorithms for updating the density operator (for example, it is strictly not necessary to store and calculate both *ρ _{ij}* and

*ρ*). At any point, the transition from a mixed-state formalism to a pure state formalism can be made by replacing |

_{ji}*ρ*〉〉 by the state vector |Ψ〉, and replacing the Liouville operator by the Hamiltonian, although such systems do not allow one to consider the effects of homogeneous relaxations.

## Acknowledgments

This research was partially supported by the Norwegian University of Science and Technology. We also wish to extend our gratitude to Prof. J.H. Eberly and the University of Rochester for allowing one of us (R.M.) to work in the Physics department for the last 6 months.

## References and links

**1. **T. Bartel, P. Gaal, K. Reimann, M. Woerner, and T. Elsaesser, “Generation of single-cycle THz transients with high electric-field amplitudes,” Opt. Lett. **30**, 2805–2807 (2005). [CrossRef] [PubMed]

**2. **R. R. Jones, D. You, and P. H. Bucksbaum, “Ionization of Rydberg atoms by subpicosecond half-cycle electromagnetic pulses,” Phys. Rev. Lett. **70**, 1236–1239 (1993). [CrossRef] [PubMed]

**3. **V. Viteau, A. Chotia, M. Allegrini, N. Bouloufa, O. Dulieu, D. Comparat, and P. Pillet, “Optical pumping and vibrational cooling of molecules,” Science **321**, 232–234 (2008). [CrossRef] [PubMed]

**4. **J. N. Sweetser and I. A. Walmsley, “Linear pulse propagation in stationary and nonstationary multilevel media in the transient regime,” J. Opt. Soc. Am. B **13**, 601–612 (1996). [CrossRef]

**5. **L. E. E. Araujo, “Ultrashort pulse propagation in multilevel systems,” Phys. Rev. A **72**, 053802 (2005). [CrossRef]

**6. **S. Hughes, “Breakdown of the Area theorem: carrier-wave Rabi flopping of femtosecond optical pulses,” Phys. Rev. A **81**, 3363–3366 (1998).

**7. **R. W. Ziolkowski, J. M. Arnold, and D. M. Gogny, “Ultrafast pulse interactions with two-level atoms,” Phys. Rev. A **52**, 3082–3094 (1995). [CrossRef] [PubMed]

**8. **K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. **14**, 302–307 (1966). [CrossRef]

**9. **J. A. Gruetzmacher and N. F. Scherer, “Finite-difference time-domain simulations of ultrashort pulse propagation incorporating quantum-mechanical response functions,” Opt. Lett. **28**, 573–575 (2003). [CrossRef] [PubMed]

**10. **S. Hughes, “Subfemtosecond soft-x-ray generation from a two-level atom: extreme carrier-wave Rabi flopping,” Phys. Rev. A **62**, 055401 (2000). [CrossRef]

**11. **V. P. Kalosha and J. Herrmann, “Formation of optical subcycle pulses and full Maxwell-Bloch solitary waves by coherent propagation effects,” Phys. Rev. Lett. **83**, 544–547 (1998). [CrossRef]

**12. **Y. Niu, K. Xia, N. Cui, S. Gong, and R. Li, “Spatiotemporal evolution and multiple self-focusing of ultrashort pulses in a resonant two-level medium,” Phys. Rev. A **78**, 063835 (2008). [CrossRef]

**13. **F. Schlottau, M. Piket-May, and K. Wagner, “Modeling of femtosecond pulse interaction with inhomogeneously broadened media using an iterative predictor-corrector FDTD method,” Opt. Express **13**, 182–194 (2005). [CrossRef] [PubMed]

**14. **B. Bidegaray, A. Bourgeade, and D. Reignier, “Introducing physical relaxations terms in Bloch equations,” J. Comp. Phys. **170**, 603–613 (2000). [CrossRef]

**15. **B. Bidegaray, “Time discretizations for Maxwell-Bloch equations,” Numer. Methods Partial Differ. Equ. **19**, 284–300 (2003). [CrossRef]

**16. **A. Bourgeade and O. Saut, “Numerical methods for the bidimensional Maxwell-Bloch equations in nonlinear crystals,” J. Comp. Phys. **213**, 823–843 (2006). [CrossRef]

**17. **S. Mukamel, *Principles of Nonlinear Optical Spectroscopy* (Oxford University Press, 1995).

**18. **Q. H. Liu, “The pseudospectral time-domain (PSTD) method: A new algorithm for solutions of Maxwell’s equations,” in *Antennas and Propagation Society International Symposium Digest* (IEEE, 1997), Vol. 1, pp. 122–125. [CrossRef]

**19. **T.-W. Lee and S. C. Hagness, “Pseudospectral time-domain methods for modeling optical wave propagation in second-order nonlinear materials,” J. Opt. Soc. Am. B **21**, 330–342 (2004). [CrossRef]

**20. **S. Rosseland, “On the transmission of radiation through an absorbing medium in motion, with applications to the theory of sun-spots and solar rotation,” Astrophys. J. **63**, 342–367 (1926). [CrossRef]

**21. **S. Blanes, F. Casas, J. A. Oteo, and J. Ros, “A pedagogical approach to the Magnus expansion,” Eur. J. Phys. **31**, 907–918 (2010). [CrossRef]

**22. **R. B. Sidje, “Expokit: a software package for computing matrix exponentials,” ACM Trans. Math. Softw. **24**, 130–156 (1998). [CrossRef]

**23. **L. Allen and J. H. Eberly, *Optical Resonance and Two-Level Atoms* (Dover, 1987).