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

H^=ωi(v)|ii|E(r,t)μij(θ,ϕ)|ij|,
with implicit summation of i and j over the N available energy levels. ħωi is the energy of the energy eigenstate |i〉, and μ ij = 〈i|μ^|j〉 is the ij element of the effective dipole moment operator μ^. The spherical angles θ, ϕ and the velocity v are used to classify absorbers that have a characteristic orientation in the θ, ϕ 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 ωi are not functions of 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

Σ=n(v,θ,ϕ)dvdΩ,
where 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 v]. The total number density of the medium is denoted by Na.

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 ρ^. Therefore, we seek to replace the von Neumann equation dtρ^ =()−1 [Ĥ, ρ^] by the vectorized version

dt|ρ=Lˇ|ρ.
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

|O=Ojk|jk,
where a sum over jk is implicit. Ojk 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
O1|O2=Tr(O^1O^2).
This implies the orthogonality 〈〈jk|mn〉〉 = δjmδkn, and a completeness relation Σjk|jk〉〉〈〈jk| = 1. Every Liouville space operator can therefore be expanded as Ǒ = Σjk,mn Ojk,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

Lˇ=iH^(H^).
The Kronecker sum ⊕ is defined AB = 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),

LˇLˇ+Rˇ.
The only fundamental restriction on Ř is that the equation dt|ρ〉〉 = Ř|ρ〉〉 must preserve the physical properties of ρ^. It is convenient to derive Ř directly from the master equation
dtρ^=(i)1[H^,ρ^]+i,j(γ^ijρ^γ^ij12ρ^γ^ijγ^ij12γ^ijγ^ijρ^),
where the operators γ^ij are of the form γ^ij=γij|ij| with γij real. Equation (8) preserves positiveness, hermiticity and trace. Generally, operators γ^ij with ij represent damping that cause population transfer from state j to state i at a rate γij(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
dtρmn=(i)1[H^,ρ^]mn12j(γjm+γjn)ρmn,mn
dtρmm=(i)1[H^,ρ^]mn+jmγmjρjjjmγjmρmm
The total off-diagonal decay rates are sums of both pure decoherence rates (γii’s) as well as amplitude decay rates (γij’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
Rˇ=i,j[γ^ijγ^ij12(γ^ijγ^ij)(γ^ijγ^ij)],
where we have made use of the fact that all γij’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/kBT 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 = γjm exp[−β(ωmωj)] where β = ħ/kBT. Independent of forcing the thermal restriction, the equilibrium state is the solution to the rate balance equation Σj mγmjρjj − Σ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γmjρjj − Σ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 = 〈Tr(μ^ ρ^)〉Σ and the time-derivative of the polarization, dt P, can be written by the use of Eq. (3), Eq. (5) and the hermiticity of μ^ as

P=Naμ|ρΣ,
dtP=Naμ|Lˇ|ρΣ.

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

|ρ(t)=Uˇ0(t,t0)exp[n=1Ωˇ(t,t0)]|ρ(t0),
where the free propagator Ǔ 0(t, t 0) is Ǔ 0(t, t 0) = exp[Ľ 0(tt 0)]. Σn=1Ωˇ(t,t0) is the Magnus series [21] where the first term reads
Ωˇ1(t,t0)=t0tdτLˇI(τ),
where LˇI(t)=Uˇ01(t,t0)Lˇ1(t)Uˇ0(t,t0). 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 tt 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 ρ^.

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 tc = Ω−1 = (Ecμ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 H is scaled according to H(Ecɛ/μ)H, r(tc/μɛ)r, ĽĽ/tc 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

tH=×E,
tE=×Hζμ|Lˇ0+ELˇ1)|ρΣ,
dt|ρ=(Lˇ0+ELˇ1)|ρ,
where ζ=1ɛNaμcEc and the Liouville operators (that contain only dimensionless frequencies and dipole moments) are
Lˇ0=iH^0(H^0)+Rˇ,
Lˇ1=iμ(μ^).

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(tt 0)] and UˇI(t,t0)=exp[Σn=1Ωˇn(t,t0)]. 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

|ρn+1/2=Uˇ0UˇIn|ρn1/2,
Uˇ0=exp(Lˇ0Δt),
UˇIn=exp(LˇInΔt),
where LˇIn=Uˇ01EnLˇ1Uˇ0. We have truncated the Magnus series at the first term and approximated the integral (n1/2)Δt(n+1/2)ΔtdτLˇI(τ)LˇInΔt. Equation (16) preserves all the physical properties of ρ^ as long as the matrix exponentials can be calculated exactly. Note that |ρ〉〉 and 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(LˇInΔ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 LˇInEnLˇ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 LˇInEnLˇ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 UˇIn can be written

UˇIn=exp(iEnμ^Δt)*exp(iEnμ^Δt),
where we have made use of the property exp(AB) = exp(A) ⊗ exp(B) and the hermiticity of μ^. The point is that instead of calculating the matrix exponential exp (–i E n · [μ^ ⊕ (−μ^)]Δ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 N 2. The approximations made so far essentially leads to the method introduced in [14], and the solution is summarized as
|ρn+1/2=Uˇ0[exp(iEnμ^Δt)*exp(iEnμ^Δt)]|ρn1/2.

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

|ρn+1/2=Uˇ01/2UˇInUˇ01/2|ρn1/2.
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.

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 〈〈μ| 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, 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 d by

dF=d1(2πikdd(F))
where kd is the spatial frequency in the d-direction, F is a component of either E or H, and d and d1 denote the Fourier and inverse Fourier transforms in the d-direction. Taking the numerical approximation tF ≈ (𝔻tF)n, Maxwell’s equations can be summarized in their discrete forms as
(𝔻tHx)i,j,kn=(yEzzEy)i,j,kn,
(𝔻tHy)i,j,kn=(zExxEz)i,j,kn,
(𝔻tHz)i,j,kn=(xEyyEx)i,j,kn,
(𝔻tEx)i,j,kn+1/2=(yHzzHy)i,j,kn+1/2(ζμx|(Lˇ0+ELˇ1|ρΣ)i,j,kn+1/2,
(𝔻tEy)i,j,kn+1/2=(zHxxHz)i,j,kn+1/2(ζμy|(Lˇ0+ELˇ1|ρΣ)i,j,kn+1/2,
(𝔻tEz)i,j,kn+1/2=(xHyyHx)i,j,kn+1/2(ζμz|(Lˇ0+ELˇ1|ρΣ)i,j,kn+1/2,
where |μx〉〉 = ·|μ〉〉, |μy〉〉 = ŷ ·|μ〉〉, |μz〉〉 = ·|μ〉〉. (𝔻tF)n is traditionally taken as the second order Yee-type leapfrogging (𝔻tF)n=Fn+1/2Fn1/2Δ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 (𝔻tF)n, which will allow us to push closer to the Nyquist limit. In one dimension, E = Ex, H = Hy, x = y = 0, the fourth order numerical approximation for t F can be written (𝔻tF)n=Fn+1/2Fn1/2ΔtΔt2241[(2πik)3(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)/( 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.

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 (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 Na = 1024m−3, μc = 10−29Cm and ω 0 = 4π · 1014rad/s, which reflects the choice of parameters made in [7].

The injected pulse is E(t) = (t)sin(ω 0 t) = Ecsech [(t – 6τ)/τ]sin(ω 0 t) where Ec = 4.2186·109V/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Δzcπ0.318fs. The chosen values for Ec, μc and τ reflect that the pulse area is A=μcε(t)dt=2π, which is the necessary condition for the lossless propagation of the 2π hyperbolic secant soliton shown in Fig. 1.

 figure: 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 |(t, z)| of the electric field E(t,z)= (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 counter-rotating terms in the matter-field interaction [7].

Download Full Size | PPT Slide | PDF

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

g(v)=1πvp2exp(v2vp2).
Inserting the detuning ω = ωo(1 + v/c) = ω 0 + Δω into Eq. (23) and integrating the free induction decay P(t)=d(Δω)g(Δω)cos[(ω0+Δω)t] we can extract a Doppler lifetime T2*=2cvp1ω0. In order to demonstrate and validate the role of the mean free velocity vp, which acts as a control parameter of the inhomogeneous bandwidth, we choose to simulate three materials with vp/c = 1/(50π), vp/c = 1/(100π) and vp/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 T2*=250fs, T2*=500fs and T2*=1ps. The inhomogeneous bandwidth ( 1/T2*) lies in the THz range.

Numerically, we represent the broadening function by sampling nl = 1001 absorption lines in the velocity interval vc[−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)=Ec4(sech[(t50τ)/τ]+sech[(t800τ)/(2π)])sin(ω0t), 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 T2*’s agree very well with the theoretically predicted values of 250fs, 500fs and 1ps, which serves to validate our method.

 figure: 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) = 𝒫(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).

Download Full Size | PPT Slide | PDF

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π · 1013rad/s. The injected pulse is E(t) = Ec exp[−(tt 0)2/τ 2]sin(ω 0 t) where τ = 100fs, t 0 = 3τ and Ec = 5 · 108 V/m. This pulse has a spectral content that covers all the resonances ωn +1ωn, and the peak of the spectrum is located at the ω 4ω 3 resonance. We include the Doppler broadening function in Eq. (23) sampled at 301 Doppler lines with vp/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 γ121=1ps, γ231=1.1ps, γ341=1.2ps, γ451=1.3ps, γ561=1.4ps and γ111=γ221==γ661=1ps. We force the thermal restriction γmj = γjm exp[−β (ωmωj)] with 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 Na = 1025m−3 and take μ 12 = μ 23 = ...μ 56 = 10−29Cm. The length of our simulation region is L = 1mm with spatial and temporal time steps Δz = L/500 = 2μm and Δt = 2Δz/() ≈ 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) Ω¯=A/τ=(Ecμcπ)/ 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: Fig. 3

Fig. 3 Few-cycle pulse propagation shown in a 6-level ladder configuration. The left figure shows the rapid absorption of the pulse in the medium. Even after only one pulselength the peak amplitude has been reduced to half its initial amplitude. The right figure shows the pulse spectrum during propagation.

Download Full Size | PPT Slide | PDF

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.

 figure: 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.

Download Full Size | PPT Slide | PDF

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 ρ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.

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).

References

  • View by:
  • |
  • |
  • |

  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).

2010 (1)

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]

2008 (2)

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]

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]

2006 (1)

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

2005 (3)

2004 (1)

2003 (2)

2000 (2)

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

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

1998 (3)

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]

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

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

1996 (1)

1995 (1)

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]

1993 (1)

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]

1966 (1)

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]

1926 (1)

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]

Allegrini, M.

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]

Allen, L.

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

Araujo, L. E. E.

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

Arnold, J. M.

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]

Bartel, T.

Bidegaray, B.

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

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

Blanes, S.

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]

Bouloufa, N.

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]

Bourgeade, A.

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

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

Bucksbaum, P. H.

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]

Casas, F.

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]

Chotia, A.

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]

Comparat, D.

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]

Cui, N.

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]

Dulieu, O.

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]

Eberly, J. H.

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

Elsaesser, T.

Gaal, P.

Gogny, D. M.

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]

Gong, S.

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]

Gruetzmacher, J. A.

Hagness, S. C.

Herrmann, J.

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]

Hughes, S.

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

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

Jones, R. R.

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]

Kalosha, V. P.

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]

Lee, T.-W.

Li, R.

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]

Liu, Q. H.

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]

Mukamel, S.

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

Niu, Y.

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]

Oteo, J. A.

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]

Piket-May, M.

Pillet, P.

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]

Reignier, D.

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

Reimann, K.

Ros, J.

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]

Rosseland, S.

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]

Saut, O.

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

Scherer, N. F.

Schlottau, F.

Sidje, R. B.

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

Sweetser, J. N.

Viteau, V.

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]

Wagner, K.

Walmsley, I. A.

Woerner, M.

Xia, K.

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]

Yee, K.

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]

You, D.

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]

Ziolkowski, R. W.

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]

ACM Trans. Math. Softw. (1)

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

Astrophys. J. (1)

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]

Eur. J. Phys. (1)

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]

IEEE Trans. Antennas Propag. (1)

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]

J. Comp. Phys. (2)

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

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

J. Opt. Soc. Am. B (2)

Numer. Methods Partial Differ. Equ. (1)

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

Opt. Express (1)

Opt. Lett. (2)

Phys. Rev. A (5)

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

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

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

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]

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]

Phys. Rev. Lett. (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]

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]

Science (1)

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]

Other (3)

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

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]

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

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (4)

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 |(t, z)| of the electric field E(t,z)= (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 counter-rotating 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) = 𝒫(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).
Fig. 3
Fig. 3 Few-cycle pulse propagation shown in a 6-level ladder configuration. The left figure shows the rapid absorption of the pulse in the medium. Even after only one pulselength the peak amplitude has been reduced to half its initial amplitude. The right figure shows the pulse spectrum during propagation.
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.

Equations (34)

Equations on this page are rendered with MathJax. Learn more.

H ^ = ω i ( v ) | i i | E ( r , t ) μ i j ( θ , ϕ ) | i j | ,
Σ = n ( v , θ , ϕ ) d v d Ω ,
d t | ρ = L ˇ | ρ .
| O = O j k | j k ,
O 1 | O 2 = Tr ( O ^ 1 O ^ 2 ) .
L ˇ = i H ^ ( H ^ ) .
L ˇ L ˇ + R ˇ .
d t ρ ^ = ( i ) 1 [ H ^ , ρ ^ ] + i , j ( γ ^ i j ρ ^ γ ^ i j 1 2 ρ ^ γ ^ i j γ ^ i j 1 2 γ ^ i j γ ^ i j ρ ^ ) ,
d t ρ m n = ( i ) 1 [ H ^ , ρ ^ ] m n 1 2 j ( γ j m + γ j n ) ρ m n , m n
d t ρ m m = ( i ) 1 [ H ^ , ρ ^ ] m n + j m γ m j ρ j j j m γ j m ρ m m
R ˇ = i , j [ γ ^ i j γ ^ i j 1 2 ( γ ^ i j γ ^ i j ) ( γ ^ i j γ ^ i j ) ] ,
P = N a μ | ρ Σ ,
d t P = N a μ | L ˇ | ρ Σ .
| ρ ( t ) = U ˇ 0 ( t , t 0 ) exp [ n = 1 Ω ˇ ( t , t 0 ) ] | ρ ( t 0 ) ,
Ω ˇ 1 ( t , t 0 ) = t 0 t d τ L ˇ I ( τ ) ,
t H = × E ,
t E = × H ζ μ | L ˇ 0 + E L ˇ 1 ) | ρ Σ ,
d t | ρ = ( L ˇ 0 + E L ˇ 1 ) | ρ ,
L ˇ 0 = i H ^ 0 ( H ^ 0 ) + R ˇ ,
L ˇ 1 = i μ ( μ ^ ) .
| ρ n + 1 / 2 = U ˇ 0 U ˇ I n | ρ n 1 / 2 ,
U ˇ 0 = exp ( L ˇ 0 Δ t ) ,
U ˇ I n = exp ( L ˇ I n Δ t ) ,
U ˇ I n = exp ( i E n μ ^ Δ t ) * exp ( i E n μ ^ Δ t ) ,
| ρ n + 1 / 2 = U ˇ 0 [ exp ( i E n μ ^ Δ t ) * exp ( i E n μ ^ Δ t ) ] | ρ n 1 / 2 .
| ρ n + 1 / 2 = U ˇ 0 1 / 2 U ˇ I n U ˇ 0 1 / 2 | ρ n 1 / 2 .
d F = d 1 ( 2 π i k d d ( F ) )
( 𝔻 t H x ) i , j , k n = ( y E z z E y ) i , j , k n ,
( 𝔻 t H y ) i , j , k n = ( z E x x E z ) i , j , k n ,
( 𝔻 t H z ) i , j , k n = ( x E y y E x ) i , j , k n ,
( 𝔻 t E x ) i , j , k n + 1 / 2 = ( y H z z H y ) i , j , k n + 1 / 2 ( ζ μ x | ( L ˇ 0 + E L ˇ 1 | ρ Σ ) i , j , k n + 1 / 2 ,
( 𝔻 t E y ) i , j , k n + 1 / 2 = ( z H x x H z ) i , j , k n + 1 / 2 ( ζ μ y | ( L ˇ 0 + E L ˇ 1 | ρ Σ ) i , j , k n + 1 / 2 ,
( 𝔻 t E z ) i , j , k n + 1 / 2 = ( x H y y H x ) i , j , k n + 1 / 2 ( ζ μ z | ( L ˇ 0 + E L ˇ 1 | ρ Σ ) i , j , k n + 1 / 2 ,
g ( v ) = 1 π v p 2 exp ( v 2 v p 2 ) .

Metrics