Multilevel Maxwell-Bloch simulations in inhomogeneously broadened media

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 OCIS codes: (030.1670) Coherent optical effects; (160.2710) Inhomogeneous optical media; (320.7110) Ultrafast nonlinear optics. 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). 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). 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). 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). 5. L. E. E. Araujo, “Ultrashort pulse propagation in multilevel systems,” Phys. Rev. A 72, 053802 (2005). 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). 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). 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). 10. S. Hughes, “Subfemtosecond soft-x-ray generation from a two-level atom: extreme carrier-wave Rabi flopping,” Phys. Rev. A 62, 055401 (2000). 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). 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). 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). 14. B. Bidegaray, A. Bourgeade, and D. Reignier, “Introducing physical relaxations terms in Bloch equations,” J. Comp. Phys. 170, 603–613 (2000). 15. B. Bidegaray, “Time discretizations for Maxwell-Bloch equations,” Numer. Methods Partial Differ. Equ. 19, 284–300 (2003). 16. A. Bourgeade and O. Saut, “Numerical methods for the bidimensional Maxwell-Bloch equations in nonlinear crystals,” J. Comp. Phys. 213, 823–843 (2006). #147292 $15.00 USD Received 10 May 2011; revised 23 Jul 2011; accepted 27 Jul 2011; published 15 Aug 2011 (C) 2011 OSA 29 August 2011 / Vol. 19, No. 18 / OPTICS EXPRESS 16784 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. 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). 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). 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). 22. R. B. Sidje, “Expokit: a software package for computing matrix exponentials,” ACM Trans. Math. Softw. 24, 130–156 (1998). 23. L. Allen and J. H. Eberly, Optical Resonance and Two-Level Atoms (Dover, 1987).


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.

Theoretical model
The Hamiltonian Ĥ = Ĥ(r r r,t,v v v,θ ,ϕ) for an absorber interacting with a classical electromagnetic field E E E = E E E(r r r,t) under the dipole approximation and expanded in the energy basis is with implicit summation of i and j over the N available energy levels.̵ hω i is the energy of the energy eigenstate |i⟩, and μ μ μ i j = ⟨i| μ μ μ| j⟩ is the i j element of the effective dipole moment operator μ μ μ.The spherical angles θ ,ϕ and the velocity v v v are used to classify absorbers that have a characteristic orientation in the θ ,ϕ direction and a spatial velocity v v v. Ĥ is implicitly understood to belong to one of these classes in the unit volume at (r r r,t).It is natural to assume that μ μ μ and ω i are not functions of r r r and 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 where n(v v v,θ ,ϕ)dv v vdΩ is the fraction of absorbers with orientation (θ ,ϕ) within the solid angle dΩ, and a spatial velocity in the interval The total number density of the medium is denoted by N a .

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 ρ.Therefore, we seek to replace the von Neumann equation ] by the vectorized version In a Hilbert space of dimension 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 ρ, 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 where a sum over jk is implicit.O jk is the 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 This implies the orthogonality ⟪ jk |mn⟫ = δ jm δ kn , and a completeness relation 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 The Kronecker sum ⊕ is defined where ⊗ is the Kronecker product.

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), The only fundamental restriction on Ř is that the equation d t |ρ⟫ = Ř|ρ⟫ must preserve the physical properties of ρ.It is convenient to derive Ř directly from the master equation where the operators γij are of the form γij = √ γ i j |i⟩⟨ j| with γ i j real.Equation ( 8) preserves positiveness, hermiticity and trace.Generally, operators γij with i ≠ j represent damping that cause population transfer from state j to state i at a rate γ i j (for example due to inelastic collisions, spontaneous emission or a thermal background).Operators γii cause decoherence only (for example due to elastic collisions).The off-diagonal and diagonal terms decay according to The total off-diagonal decay rates are sums of both pure decoherence rates (γ ii 's) as well as amplitude decay rates (γ i j 's).The effective decay rates for the populations ρ mm 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 where we have made use of the fact that all γ i j 's are real.
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 γ m j = γ jm exp[−β (ω m − ω j )] where β = ̵ h/k B T .Independent of forcing the thermal restriction, the equilibrium state is the solution to the rate balance equation ∑ j≠m γ m j ρ j j − ∑ j≠m γ jm ρ mm = 0, and all off-diagonal elements equal to zero.This can be written as the linear system Γ|diag( ρ)⟩ = 0 where |diag( ρ)⟩ is the column vector corresponding to the diagonal elements of ρ and Γ is the matrix corresponding to the linear system ∑ j≠m γ m j ρ j j − ∑ j≠m γ jm ρ mm = 0.The solution is known [20], and the populations are ρ mm = Γ jm ∑ j Γ jm where Γ jm is the cofactor of the element Γ jm .Γ has the special shape that all cofactors in every column are identical (Γ 1m = Γ 2m = ... etc.), so the cofactor Γ jm is independent of j.
With Ř now included implicitly in Ľ, the polarization P P P = ⟨Tr( μ μ μ ρ)⟩ ∑ and the time-derivative of the polarization, d t P P P, can be written by the use of Eq. ( 3), Eq. ( 5) and the hermiticity of μ μ μ as

Time evolution
By splitting the Liouville operator into time dependent and time independent parts, Ľ = Ľ0 + Ľ1 (t), the exact solution to Eq. ( 3) is where the free propagator Ǔ0 (t, ) is the Magnus series [21] where the first term reads where ĽI (t) = Ǔ−1 0 (t,t 0 ) Ľ1 (t) Ǔ0 (t,t 0 ).We shall assume that the higher order terms, which contains the time integral of the commutator [ Ľ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 ρ.

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 = ω i0 (1 + v/c) where ω i0 is the energy divided by ̵ h 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.

Maxwell's equations and dimensionless units
To ensure maximal precision of numerical simulations, we introduce characteristic scales such that E E E is of order one, and a time-scale t c = Ω −1 = (E c μ c / ̵ h) −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 where E c and the Liouville operators (that contain only dimensionless frequencies and dipole moments) are

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 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 where Ľn . We have truncated the Magnus series at the first term and approximated the integral ∫ (n+1/2)Δt (n−1/2)Δt dτ ĽI (τ) ≈ Ľn I Δt.Equation ( 16) preserves all the physical properties of ρ as long as the matrix exponentials can be calculated exactly.Note that |ρ⟫ and E E E are stored on staggered temporal grids, which is the general trend throughout this paper.
There are several approximations and simplifications that can be made in Eq. ( 16).Firstly, the computation of the matrix exponential exp( Ľn I Δ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 Ľn 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 + O(Δt).

If inhomogeneous broadening is included, taking Ľn
I ≈ E E E n ⋅ Ľ L 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 Ǔn I can be written where we have made use of the property exp ) and the hermiticity of μ μ μ.The point is that instead of calculating the matrix exponential exp(−iE where the argument is of order N 4 , we have reduced it to calculating the exponential exp(iE E E n ⋅ μ μ μ) where the argument is of order 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), Equation ( 19) is an example of a second order approximation to the solution for |ρ 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.

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 E 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 |ρ⟫ = ETr( μ,[ μ, ρ]) = 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 E E, H H H and |ρ⟫, 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 where k d is the spatial frequency in the d-direction, F is a component of either E E E or H H H, and F d and F −1 d denote the Fourier and inverse Fourier transforms in the d-direction.Taking the numerical approximation ∂ t F ≈ (D t F) n , Maxwell's equations can be summarized in their discrete forms as where ) n is traditionally taken as the second order Yee-type leapfrogging

Δ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 (D t F) n , which will allow us to push closer to the Nyquist limit.In one dimension, E E E = E x , H H H = H y , ∂ x = ∂ y = 0, the fourth order numerical approximation for ∂ t F can be written [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)/( √ Dπ) (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.

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 .

Two-level systems
Our first numerical simulation shows how an energetic electromagnetic pulse can travel losslessly 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 (2n + 1)π pulse inverts the medium through n inversions and a (2n + 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 , μ c = 10 −29 Cm and ω 0 = 4π ⋅ 10 14 rad/s, which reflects the choice of parameters made in [7].(t, z) of the medium.The "flattenings" on the inversion bump are manifestations of counterrotating terms in the matter-field interaction [7].
The injected pulse is E(t) = E(t)sin(ω 0 t) = E c sech[(t − 6τ)/τ]sin(ω 0 t) where E c = 4.2186 ⋅ 10 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 Δz = L/250 = λ 10 = 0.15μm and Δt = 2Δz cπ ≈ 0.318fs.The chosen values for E c , μ c and τ reflect that the pulse area is A 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 Inserting the detuning ω = ω o (1+v/c) = ω 0 +Δω into Eq.( 23) and integrating the free induction decay . In order to demonstrate and validate the role of the mean free velocity v p , which acts as a control parameter of the inhomogeneous bandwidth, we choose to simulate three materials with 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.approximation, the SVEA and RWA.

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 inhomoge-neous 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 ρ i j and ρ ji ).At any point, the transition from a mixed-state formalism to a pure state formalism can be made by replacing |ρ⟫ 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.
and |μ μ μ⟫ → μ 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

ρFig. 1 .
Fig. 1.Propagation of a 2π hyperbolic secant soliton shown in a two-level absorber of length L = 37.5μm.The left figure shows the envelope |E(t, z)| of the electric field E(t, z) = E(t, z)cos φ (t, z).The envelope has been extracted from the full-field data by use of a Hilbert transform.The right figure shows the traveling inversion bump ρ 22 (t, z) − ρ 11(t, z) of the medium.The "flattenings" on the inversion bump are manifestations of counterrotating terms in the matter-field interaction[7].

Fig. 2 .ω
Fig. 2. Photon echo(es) shown in a two-level absorber for a simulation with Δt = 0.318fs for a total duration of 11ps.The top figure shows the injected pulse envelopes obtained through a Hilbert transform of the full-field, while the inset shows the population inversion ⟨ρ 22 − ρ 11 ⟩ ∑ during the passage of the first pulse.The second figure (below) shows the envelope(s) of the polarization(s) P(t) = P(t) cos ϕ(t) obtained through a Hilbert transform for vp/c = 1/(50π)(dotted), vp/c = 1/(100π)(dashed) and vp/c = 1/(200π)(solid).Numerically, we represent the broadening function by sampling n l = 1001 absorption lines in

Fig. 4 .
Fig. 4. The figure on the left shows the normalized energy during the passage of the pulse at the entrance of the medium, while the figure on the right shows the level populations.
The theoretically derived Doppler lifetimes are T *