An efficient quantum memory based on two-level atoms

Light is an ideal carrier of information both in the classical and in the quantum regime. To harness the full potential of light for quantum information processing it is, however, a major advantage to have access to quantum memories capable of storing the information from the light into an atomic medium and later releasing the information. Ideally such quantum memories should introduce as little disturbance as possible to the information encoded in the light field, i.e. the memory should be fully coherent and efficient. Furthermore it should be as simple as possible to operate.

A large number of proposals have been developed for how one can construct quantum memories based on atomic ensembles [1]. As opposed to the most natural approach of storing one photon in one atom, the idea behind the ensemble-based approach is to store photons into the collective states of an ensemble of atoms. Thereby, one avoids the technical challenges associated with efficiently coupling single atoms and single photons. As a consequence the ensemble-based approach considerably simplifies the experimental realization of quantum memories. A large class of quantum memory approaches based on atomic ensembles uses classical laser fields to control the memory process in such a way that the incoming field is mapped into the ground-state coherence of the atoms. Based on this a number of key experimental advances have been achieved (see for instance [2–8]). Here we shall pursue a different approach based on controlled reversible inhomogeneous broadening (CRIB) [9]. The original (transverse) CRIB is developed for impurities embedded in solid-state systems. Due to the random orientations of the (static) dipole moments of the atoms or a transverse field gradient, an external applied field will lead to different shifts in the transition frequency of different atoms. Hence an external applied field will dephase the polarization of the atoms. This will in essence turn off the reemission of the light absorbed into the ensemble which is therefore stored. The central idea in the CRIB approach is to later reverse the direction of the external field. If the atomic dipoles are fixed, the reversal of the field generates a shift in the opposite direction (see figure  1(c)). After a certain time the polarization will rephase, causing an echo, where the light is reemitted. In the original approach an additional ground level was used and it was shown that by driving the transition to this level with laser fields one can obtain a read-out in the backward direction which can, in principle, reach 100% efficiency (within the theoretical models we shall consider here, the efficiency is the only parameter characterizing the performance of the memory for a single incident mode [1, 10], and we shall therefore focus on this parameter below). From a practical perspective it is highly desirable to avoid the use of a third level and the associated lasers to drive that transition.

Zoom In Zoom Out Reset image size

Since the original CRIB proposal [9] (see also [11]) a number of modifications have appeared [12–14] and much experimental progress has been reported [12, 15–20]. In particular, a different longitudinal CRIB approach has been developed where a gradient of an external field is applied such that the shift depends on the position of the atoms along the propagation direction [12]. For longitudinal CRIB, it has been shown that one can construct an efficient memory based on two-level systems which, in principle, can reach an efficiency of 100% [13] if the atoms have a sufficiently long coherence time. For transverse CRIB, it was also shown that one can construct a quantum memory based on two-level systems, but in this case the efficiency was found to be limited to 54% [21]. Here we develop the theory for an efficient quantum memory based on two-level atoms subject to transverse CRIB. We consider varying the broadening in time so that it is not turned on initially but only applied during a short time interval as shown in figure 1(d). We show that with this procedure the efficiency of the memory can, in principle, reach 100% if the optical depth of the ensemble is very high, the atoms have a very long coherence time, and pulses with a large bandwidth are incident. A highly related proposal was considered in [22], where the broadening was also applied after a certain time as in our case. Since that work only considered narrow bandwidth pulses, similar to [21], the efficiency was still limited to 54%. We believe that the use of wide bandwidth pulses, as we consider here, may help overcome that limitation.

Our system consists of N two-level atoms that are coupled to a quantized electric field. We treat everything in one dimension and thus the electric field is only dependent on the spatial coordinate z. The atoms are assumed to be confined to a length L along the z-axis and are distributed with uniform density. The transverse extent of the ensemble and the cross-section of the electric field mode are assumed to be the same and have area A. See figure 1(a).

The nth atom has two internal levels denoted by |g〉_n and |e〉_n for the ground and excited state, respectively. The two levels are connected by a transition with a matrix element $\wp$. We assume that the atomic transition frequencies are affected by both an intrinsic inhomogeneous broadening that we cannot control and a controlled inhomogeneous broadening that we can rapidly turn on and off and reverse in a time much shorter than the duration of the incoming pulse. These two types of broadening give rise to a detuning of the atomic transition from the central frequency ω_L of the incoming light pulse. The two detunings for the nth atom are denoted by Δ_0n for the intrinsic broadening and Δ_n for the controlled broadening. Since we apply a reversed controlled broadening in stage 4 compared to stage 2, the transition frequency is given by ω_eg,n = ω_L + Δ_0n ± Δ_n during these two stages. See figure 1(b).

We are interested in storing light in two-level systems. Such storage is only possible if the lifetime of the excited state is very long. Therefore, we shall assume that the excited state has a sufficiently long lifetime so that spontaneous emission can be neglected. The only source of decoherence of the atomic transition which we consider here comes from having an intrinsic inhomogeneous broadening of the atomic transition. Note that this assumption of negligible spontaneous emission from the excited state does not contradict the assumption that light couples to the transition. The quantum light we are interested in couples to collective states that can have much higher emission rates due to collective enhancement effects. The memory can thus have a large bandwidth even if the individual atoms have a negligible decay rate. Indeed as discussed below, the bandwidth of the memory is characterized by the quantity

$$\begin{eqnarray*} &&\mu = \frac{Ng_0^2\wp^2}{\hbar^2}, \end{eqnarray*}$$

where g₀ can be written in terms of speed of light in vacuum c and vacuum permittivity ₀ as $g_0=\sqrt {\hbar \omega _{\mathrm { L}}/2c\epsilon _0 A}$. The bandwidth μ can be related to more experimentally accessible variables if we consider a specific model for the broadening. Throughout this work we shall assume a Gaussian distribution of the intrinsic broadening (as defined by (4)). In this case the decay of the collective atomic polarization due to dephasing is described by a Gaussian in time and we define the coherence time T₂ as the 1/e decay time of the atomic polarization. The precise definition and a more detailed discussion can be found in appendix A. The memory bandwidth can then be related to the optical depth d₀ that will be observed when the controlled broadening is turned off, and the coherence time through $\mu =d_0/(T_2\sqrt {\pi })$.

As was shown in [23, 24], 1/μ characterizes the fastest time scale for variations in the pulse shape that can be allowed if the memory is to store the light with high efficiency. For a memory based on three-level Λ-type atoms controlled by a classical drive, any pulse shape with a pulse length τ_p ≫ 1/μ can be stored efficiently in the ensemble. Here we are, however, interested in two-level systems. It was also shown in [24] that there exist certain incoming pulse shapes of duration τ_p ∼ 1/μ that can be efficiently absorbed and retrieved by having the light interacting only with two-level atoms. To control the storage and retrieval process, it was proposed in [24] to rapidly apply a π-pulse to transfer the population from the excited state to an auxiliary state for long-term storage after the pulse was absorbed in the ensemble. Later another rapid π pulse can be applied to initiate the retrieval process. It was shown that for the optimal pulse shape of the incoming field the efficiency of this memory protocol is identical to the highest obtainable efficiency for the memory. Recently this storage mechanism was also considered in the context of a concrete implementation in rare-earth solid state systems [25]. Here we follow a similar approach as the one suggested in [24] but without the use of a third auxiliary level. Instead we follow a CRIB-like approach [9, 21] and apply a broadening of the atomic levels. In our scheme the light pulse is first absorbed in the unbroadened ensemble. Afterwards an externally applied inhomogeneous broadening rapidly dephases the polarization of the atomic transition. This dephasing effectively turns off the reemission of the stored light allowing for its long-term storage if there is negligible decoherence of the atomic state. At a later stage we assume that we can reverse the applied inhomogeneous broadening, resulting in a rephasing of the atoms that effectively turns on the reemission of the pulse. As we shall show below, this process allows for an efficient quantum memory if the ensemble has a large optical depth.

The memory protocol that we propose can be divided into five stages shown in figure 1(d). In stage 1 the incoming pulse with duration τ_p is incident on the atoms and gets absorbed without any controlled broadening present. In stage 2 with a (short) duration τ_d the controlled broadening is applied, which turns off the reemission of the excitation. We find that considerably higher efficiencies are obtained if we allow the incoming pulse to have a small tail into stage 2 and we thus allow this in the detailed simulations below. Stage 3 of duration τ_s constitutes the main storage period. In this stage, we assume that the external broadening is turned off. The storage duration τ_s is only limited by the coherence time T₂ of the atomic polarization. Stage 4 marks the beginning of the read-out process. Here the external broadening is reversed and applied for a time τ_d equal to the duration of the initial broadening in stage 2. The read-out of the light pulse continues in stage 5 with duration τ_p where the external broadening is again turned off.

In our protocol most of the interaction with the light happens during stages 1 and 5. These stages work according to a different principle than considered in most studies of CRIB, which consider pulses that are long (have a narrow bandwidth) compared to the width of the absorption profile. For such protocols, the important parameter is the observable optical depth that characterizes the steady state absorption of the medium. In contrast, our protocol works in a transient regime identified in [24] where the pulse duration is short compared to the width of the absorption line. In this regime the optical depth and the absorption profile do not directly account for the storage of the light which rather depends on the dispersive part of the atomic line (since the dispersive part is related to the absorption profile through the Kramers–Kronig relations, we can, however, still characterize the ensemble through the optical depth). This dispersive part slows down and completely reshapes the incoming pulse, as illustrated in figure 1(a). For each instance in time, one can identify certain incoming mode shapes that are localized inside the ensemble. These are the modes stored by our scheme when we turn on the broadening that shuts off the propagation of the stored excitation. Due to the dispersive nature of the protocol, the duration of these pulses is set by the requirement that all frequency components should get a phase shift of order unity by propagating through the ensemble. Exploiting the dispersive part of the atomic line rather than the absorptive part is key to overcoming the 54% limitation for forward CRIB identified in [21]: because the dispersive part falls off more slowly than the absorption, we can work with pulses that are wider in frequency than the absorption line. Since most of the light is not inside the absorption peak, it is not absorbed by the ensemble allowing the efficiency to approach 100% in the ideal limit.

Furthermore most approaches to CRIB-based memories [12, 21] do not include the stage 3 introduced here but instead use longer stages 2 and 4. The idea behind this stage is to increase the flexibility of the read-out. If the broadening is applied during the whole storage period as in figure 1(a), the time it takes from when we reverse the broadening to when the pulse is retrieved is equal to the time from the storage until the reversal. Thus a long storage time gives a long response time of the memory. On the other hand, as soon as the atoms are sufficiently dephased in stage 2, the reemission of the light is turned off and any further dephasing will not affect the performance of the memory. Hence we might as well turn off the broadening after the dephasing time τ_d, which can be much smaller than the storage time τ_s. Stage 3 cannot be longer than T₂ if the memory is to be efficient but it can be arbitrarily short. If we want to retrieve the excitation stored in the atoms at any point in stage 3, we can turn on the reversed broadening (stage 4) and retrieve the excitations after a time τ_d that is independent of the storage time τ_s. Additionally, having a small τ_d makes it easier to do the numerical simulations of the memory as discussed below.

In this section we present the model that we shall use to evaluate the performance of the memory. Here we shall be rather brief in describing the formalism. A more thorough derivation of the equations of motion and the principles used is given in [1, 21]. To describe the light we consider the electric field operator given in the Schrödinger picture by

$$\begin{eqnarray*} &&\hat{\mathcal{E}}(z) =g_0\hat{\mathcal{E}}_{\rm slow}(z){\rm e}^{{\rm i}\omega_{\rm L} z/c} +{\rm h.c.} \end{eqnarray*}$$

in terms of its slowly varying components that satisfy

$$\begin{eqnarray*} &&[\hat{\mathcal{E}}_{\rm slow}(z), \hat{\mathcal{E}}_{\rm slow}^\dagger(z')]=c\delta(z-z'). \end{eqnarray*}$$

To be able to describe a realistic situation where the atomic states may suffer from decoherence, we introduce a finite coherence time of the atoms by including a distribution of the intrinsic detunings. Often such decoherence is introduced by adding a decay term to the coherence. Such a model can, however, only account for Markovian (exponential) decay of the coherence, whereas a finite line width from inhomogenous broadening may lead to different kinds of decay. Such different decoherence models can have rather different results for the memory performance [26]. As transverse CRIB protocols are most relevant for state systems that typically show strong inhomogeneous broadening, we chose to include this type of decoherence. The atomic transition is thus subject to two types of broadening, intrinsic broadening and controlled reversible broadening. The probability distributions of the detunings are described by the functions $\mathcal {G}_0$ for the intrinsic broadening and $\mathcal {G}$ for the controlled broadening. The distributions are assumed to be even functions and are normalized according to $\int _{-\infty }^\infty \mathcal {G}_0(\Delta _0){\mathrm { d}}\Delta _0=1$ and $\int _{-\infty }^\infty \mathcal {G}(\Delta ){\mathrm { d}}\Delta =1$.

To derive the equations of motion, we split the position and detunings of the intrinsic and controlled broadening into intervals of size z_δ, Δ_0δ and Δ_δ, respectively, and define the atomic polarization operators

$$\begin{eqnarray*} &&\hat{\mathcal{S}}_{klm,ij}=\frac{\sqrt{N}}{N_{klm}} \sum_{n=1}^{N_{klm}} |i\rangle_n \langle j|_n. \end{eqnarray*}$$

Here the sum is over atoms having a position in the interval [kz_δ,(k + 1)z_δ], an intrinsic detuning in the interval [lΔ_0δ,(l + 1)Δ_0δ] and a controlled detuning in the interval [mΔ_0δ,(m + 1)Δ_0δ]. The number of atoms that have positions and detunings in the mentioned intervals is $N_{klm} = (N/L)z_\delta \mathcal {G}_0\left (l\Delta _{0\delta }\right ) \Delta _{0\delta }\mathcal {G}\left (m\Delta _{\delta }\right ) \Delta _{\delta }$. In the limit z_δ,Δ_0δ,Δ_δ → 0, the operator $\hat {\mathcal {S}}_{klm,ij}$ gets replaced by $\hat {\mathcal {S}}_{ij}(z,\Delta _0,\Delta )$. The commutation relation is then

$$\begin{eqnarray*}&&\fl [\hat{\mathcal{S}}_{ij}(z,\Delta_0,\Delta), \hat{\mathcal{S}}_{i'j'}(z',\Delta_0',\Delta')] =\frac{L} {\sqrt{N}\mathcal{G}_0(\Delta_0)\mathcal{G}(\Delta)} \delta(z-z')\delta(\Delta_0-\Delta_0')\delta(\Delta-\Delta')\\ &&\times(\delta_{ji'}\hat{\mathcal{S}}_{ij'}(z,\Delta_0,\Delta) -\delta_{j'i}\hat{\mathcal{S}}_{i'j}(z,\Delta_0,\Delta)). \end{eqnarray*}$$

We shall assume that all atoms are initially prepared in the ground state |g〉 and that only weak fields are incident such that most of the population remains in the ground state. When finding commutators we therefore approximate $\hat {\mathcal {S}}_{gg}(z,\Delta _0,\Delta )\approx \sqrt {N}$ and $\hat {\mathcal {S}}_{ee}(z,\Delta _0,\Delta )\approx 0$. Hence the commutator of $\hat {\mathcal {S}}_{ge}$ and $\hat {\mathcal {S}}_{{eg}}$ is a constant, and $\hat {\mathcal {S}}_{ge}$ can be regarded as proportional to an effective harmonic oscillator annihilation operator.

If we introduce the slowly varying operator $\hat {\sigma }(z,\Delta _0,\Delta )={\mathrm { e}}^{-{\mathrm { i}}\omega _{\mathrm { L}} z/c}\hat {\mathcal {S}}_{ge}(z,\Delta _0,\Delta )$, the Hamiltonian in the interaction picture and rotating wave approximation can be written as

$$\begin{eqnarray*} &&\skew3\hat{H}=\frac{\sqrt{N}}{L} \int_{-\infty}^\infty\int_{-\infty}^\infty\int_{0}^{L} \mathcal{G}_0(\Delta_0)\mathcal{G}(\Delta)\{ \hbar\left(\Delta_0+\Delta\right) \hat{\mathcal{S}}_{ee}(z,\Delta_0,\Delta)\\ &&-g_0 \wp[\hat{\mathcal{E}}_{\rm slow}(z)\hat{\sigma}^\dagger(z,\Delta_0,\Delta) +\hat{\mathcal{E}}_{\rm slow}^\dagger(z)\hat{\sigma}(z,\Delta_0,\Delta)]\} {\rm d} z\, {\rm d}\Delta_0 \,{\rm d}\Delta. \end{eqnarray*}$$

This Hamiltonian describes the evolution during stage 2 where the controlled broadening is turned on. In stages 1, 3 and 5 where there is no controlled broadening, we can find the Hamiltonian by removing the term with $\Delta \hat {\mathcal {S}}_{ee}(z,\Delta _0,\Delta )$ in the integrand. In stage 4, where the detuning is reversed we change the sign of this term. From the Hamiltonian, we can find the equations of motion for the relevant operators. Assuming most of the population to be in the ground state, the equations of motion are given by

$$\begin{eqnarray*} &&\left(\frac{\partial}{\partial t}+c\frac{\partial}{\partial z}\right) \hat{\mathcal{E}}_{\rm slow}(z,t) ={\rm i}\sqrt{\mu}\frac{c}{L} \int_{-\infty}^\infty \int_{-\infty}^\infty \mathcal{G}_0(\Delta_0)\mathcal{G}(\Delta) \hat{\sigma}(z,t,\Delta_0, \Delta){\rm d} \Delta_0 \,{\rm d}\Delta,\\ &&\frac{\partial}{\partial t}\hat{\sigma}(z,t,\Delta_0,\Delta) =-{\rm i}(\Delta_0+\Delta)\hat{\sigma}(z,t,\Delta_0,\Delta) +{\rm i}\sqrt{\mu}\hat{\mathcal{E}}_{\rm slow}(z,t). \end{eqnarray*}$$

Since these equations are linear and only couple effective annihilation operators we can ignore the hats on the operators and consider them to be equations of functions instead of operators. This allows us to calculate any normally ordered product of operators, and evaluate the efficiency of the memory. If we only consider a single incoming mode as we shall do throughout this paper, this efficiency is the only important parameter for characterizing the performance of the memory [1, 24].

Using the bandwidth of the memory μ defined above, we can change to dimensionless units. We introduce the spatial coordinate ${\tilde {z}}=z/L$ and the time ${\tilde {t}}=\mu (t-z/c)$. We also introduce the dimensionless detunings $\tilde {\Delta }_0 = \Delta _0/\mu$ and $\tilde {\Delta } = \Delta /\mu$. Then we can define the dimensionless broadening distributions G₀ and G and the dimensionless electric field E. These are related to the old quantities by $G_0(\tilde {\Delta }_0)=\mu \mathcal {G}_0(\mu \tilde {\Delta }_0)$, $G(\tilde {\Delta })=\mu \mathcal {G}(\mu \tilde {\Delta })$ and $E(z,t)=\mathcal {E}_{\mathrm { slow}}(z,t)/\sqrt {\mu }$. Using these dimensionless variables, the equations of motion become

$$\begin{equation} \frac{\partial}{\partial {\tilde{z}}}E({\tilde{z}},{\tilde{t}}) ={\rm i}\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G_0(\tilde{\Delta}_0)G(\tilde{\Delta}) \sigma({\tilde{z}},{\tilde{t}},\tilde{\Delta}_0,\tilde{\Delta}) {\rm d}\tilde{\Delta}_0 \,{\rm d}\tilde{\Delta}, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial}{\partial {\tilde{t}}}\sigma({\tilde{z}},{\tilde{t}},\tilde{\Delta}_0,\tilde{\Delta}) =-{\rm i}(\tilde{\Delta}_0+\tilde{\Delta})\sigma({\tilde{z}},{\tilde{t}},\tilde{\Delta}_0,\tilde{\Delta}) +{\rm i} E({\tilde{z}},{\tilde{t}}). \end{equation} \tag{ 2 }$$

To describe reversal of the controlled broadening, we can replace $\tilde {\Delta }$ by $-\tilde {\Delta }$ in (1) and (2). In the stages without any broadening (1, 3 and 5), we omit the term with $\tilde {\Delta }$.

To simplify the discussions below, it will be convenient to introduce the polarization P for each value of the intrinsic detuning. This is defined by

$$\begin{equation} P({\tilde{z}},{\tilde{t}},\tilde{\Delta}_0) =\int_{-\infty}^\infty G(\tilde{\Delta})\sigma({\tilde{z}},{\tilde{t}},\tilde{\Delta}_0,\tilde{\Delta}){\rm d}\tilde{\Delta}. \end{equation} \tag{ 3 }$$

In stages 1 and 5 the equations of motion (1) and (2) can be written in terms of P instead of σ since the controlled broadening is not present (see (B.1) and (B.2) in appendix B).

We also note that in the detailed calculations we shall assume that the broadening distributions are Gaussian, i.e.

$$\begin{equation} G_0(\tilde{\Delta}_0)=\frac{1}{\sqrt{2\pi\tilde{\gamma}_0^2}} \exp\left(-\frac{\tilde{\Delta}_0^2}{2\tilde{\gamma}_0^2}\right), \quad G(\tilde{\Delta})=\frac{1}{\sqrt{2\pi\tilde{\gamma}^2}} \exp\left(-\frac{\tilde{\Delta}^2}{2\tilde{\gamma}^2}\right). \end{equation} \tag{ 4 }$$

Here the dimensionless widths are $\tilde {\gamma }_0 = \gamma _0/\mu$ and $\tilde {\gamma } = \gamma /\mu$.

Before looking at the results of the numerical simulations which provide a full assessment of the efficiency of the memory, we can gain some intuition about the important aspects of the proposed protocol by doing simple perturbative calculations. Stages 1 and 5 have been described in [24]. Here we therefore investigate the process in stages 2 and 4 where the broadening first rapidly turns off and later turns on the absorption–reemission process in such a way that a large fraction of the excitation is left in the atoms during the storage interval (stage 3). To focus on this part of the dynamics, we completely ignore the evolution in stages 1, 3 and 5 for now.

For simplicity, we shall neglect the intrinsic broadening so that our system is described by the equations

$$\begin{equation} \frac{\partial}{\partial {\tilde{z}}}E({\tilde{z}},{\tilde{t}})={\rm i}\int_{-\infty}^{\infty} G(\tilde{\Delta}) \sigma({\tilde{z}},{\tilde{t}},\tilde{\Delta}) {\rm d}\tilde{\Delta}, \end{equation} \tag{ 5 }$$

$$\begin{equation} \frac{\partial}{\partial {\tilde{t}}}\sigma({\tilde{z}},{\tilde{t}},\tilde{\Delta}) =-{\rm i}\tilde{\Delta}\sigma({\tilde{z}},{\tilde{t}},\tilde{\Delta}) +{\rm i} E({\tilde{z}},{\tilde{t}}). \end{equation} \tag{ 6 }$$

We assume that all of the input pulse was mapped to P in stage 1 so that $E({\tilde {z}}=0,{\tilde {t}})=0$. At the beginning of stage 2 the initial condition is then $\sigma ({\tilde {z}},{\tilde {t}}=0,\Delta )=P^{(1)}({\tilde {z}})$, where P⁽¹⁾ denotes the polarization at the end of stage 1. To proceed, we Laplace transform the equations (${\tilde {z}}\rightarrow u$) and combine them into

$$\begin{eqnarray*} &&\frac{\partial}{\partial {\tilde{t}}}\bar{\sigma}(u,{\tilde{t}},\tilde{\Delta}) =-{\rm i}\tilde{\Delta}\bar{\sigma}(u,{\tilde{t}},\tilde{\Delta}) -\frac{1}{u} \int_{-\infty}^{\infty} G(\tilde{\Delta}') \bar{\sigma}(u,{\tilde{t}},\tilde{\Delta}') {\rm d} \tilde{\Delta}'. \end{eqnarray*}$$

We want to solve this equation perturbatively. To do this we introduce $\sigma _S({\tilde {z}},{\tilde {t}},\tilde {\Delta })={\mathrm { e}}^{{\mathrm { i}}\tilde {\Delta }{\tilde {t}}}\sigma ({\tilde {z}},{\tilde {t}},\tilde {\Delta })$, which satisfies the differential equation

$$\begin{equation} \frac{\partial}{\partial {\tilde{t}}}\bar{\sigma}_S(u,{\tilde{t}},\tilde{\Delta}) =-\frac{1}{u}{\rm e}^{{\rm i}\tilde{\Delta}{\tilde{t}}}\int_{-\infty}^\infty G(\tilde{\Delta}')\bar{\sigma}_S(u,{\tilde{t}},\tilde{\Delta}')\,{\rm e}^{-{\rm i}\tilde{\Delta}'{\tilde{t}}}\,{\rm d}\tilde{\Delta}'. \end{equation} \tag{ 7 }$$

Under the assumption that σ_S changes slowly as a function of ${\tilde {t}}$, we can replace $\bar {\sigma }_S(u,{\tilde {t}},\Delta ')$ in the integrand on the right hand side of (7) by $\bar {\sigma }_S(u,{\tilde {t}}=0,\Delta ')=\bar {P}^{(1)}(u)$. Defining the Fourier transform of G by $\tilde {G}({\tilde {t}})=\int _{-\infty }^\infty G(\tilde {\Delta }')\,{\mathrm { e}}^{-{\mathrm { i}}\tilde {\Delta }'{\tilde {t}}} \,{\mathrm { d}}\tilde {\Delta }'$, the approximate solution for σ at the end of stage 2 is

$$\begin{equation} \bar{\sigma}^{(2)}(u,\tilde{\Delta}) ={\rm e}^{-{\rm i}\tilde{\Delta}\mu\tau_{\rm d}} \left(1-\frac{1}{u}\int_0^{\mu\tau_{\rm d}} {\rm e}^{{\rm i}\tilde{\Delta} t'}\tilde{G}(t')\,{\rm d} t'\right) \bar{P}^{(1)}(u). \end{equation} \tag{ 8 }$$

For the case of Gaussian G, we have $\tilde {G}({\tilde {t}})=\exp \left (-\tilde {\gamma }^2{\tilde {t}}^2/2\right )$ describing the dephasing in time.

Expression (8) gives the perturbative approximation for the evolution during stage 2 in the absence of an incoming pulse. In stage 4, we make the same approximation. The equations are the same as above but with $\tilde {\Delta }$ replaced by $-\tilde {\Delta }$ throughout to describe reversal of the broadening. Hence we look at the evolution of the polarization with the reversed sign of $\tilde {\Delta }$ that evolves from the initial condition (8). At the end of stage 4, the solution for σ to first order in the perturbation is then

$$\begin{eqnarray*}&&\fl \bar{\sigma}^{(4)}(u,-\tilde{\Delta}) =\left(1 -\frac{1}{u}\int_0^{\mu\tau_{\rm d}} {\rm e}^{{\rm i}\tilde{\Delta}t'}\tilde{G}(t')\,{\rm d}t' -\frac{1}{u}\int_0^{\mu\tau_{\rm d}} {\rm e}^{{\rm i}\tilde{\Delta}(t'-\mu\tau_{\rm d})} \tilde{G}(t'-\mu\tau_{\rm d})\,{\rm d}t'\right) \bar{P}^{(1)}(u). \end{eqnarray*}$$

From this expression and definition (3), we can find P. Using the assumption that G and hence $\tilde {G}$ is an even function of its argument, the total polarization at the end of stage 4 can be written as

$$\begin{equation} \bar{P}^{(4)}(u) =\left(1-\frac{2}{u}\int_0^{\mu\tau_{\rm d}} (\tilde{G}(t'))^2 {\rm d}t'\right) \bar{P}^{(1)}(u). \end{equation} \tag{ 9 }$$

The main idea behind the memory protocol we propose here is that excitations mapped into the memory during stage 1 are released during stage 5. The expression above describes how an initial excitation stored in the polarization of the atoms described by P⁽¹⁾ at the end of stage 1 is mapped to the end of stage 4. To get an idea about the performance of the broadening mechanism, we consider the efficiency with which this mapping is achieved. To do this we calculate the efficiency that can be expressed by $\eta =\int _0^1 |P^{(4)}({\tilde {z}})|^2 \,{\mathrm { d}}{\tilde {z}}$ if we assume that P⁽¹⁾ is normalized to unity ($\int |P^{(1)}({\tilde {z}})|^2 {\mathrm { d}} {\tilde {z}}=1$). Now we take the inverse Laplace transform ($u\rightarrow {\tilde {z}}$) of (9). Assuming G to be Gaussian we can carry out the integration explicitly. Keeping only terms to lowest order in the perturbation we find

$$\begin{equation} \eta(\tau_{\rm d}) =1-2\sqrt{\pi}\frac{\mu}{\gamma}{\rm erf}(\gamma\tau_{\rm d}) \int_0^1 P^{(1)}({\tilde{z}})\int_0^{\tilde{z}} P^{(1)}(z')\,{\rm d} z'\,{\rm d}{\tilde{z}}. \end{equation} \tag{ 10 }$$

This expression is the main result of this section. Ideally, we would like to have an efficiency of unity. From the expression we see that an efficient memory operation can be achieved if γ ≫ μ, i.e. if the width of the broadening is much larger than the bandwidth of the memory. In this limit the broadening is sufficiently strong to dephase the atoms and turn off the reemission of the light that happens on a time scale ∼1/μ. Alternatively this can be expressed in terms of the optical depth: as derived in appendix A, the measurable optical depth in the presence of broadening is $d\approx \sqrt {2\pi }\mu /\gamma$. Hence the limit where the broadening is efficient at storing the excitation is equivalent to the limit where the optical depth is much smaller than unity after the broadening has been turned on.

The argument of the error function in (10) reflects the remaining reemission that has not been completely dephased at the end of stage 2 (the emission rate at the end is ∼dη/dτ_d). Hence if τ_d ≫ 1/γ the reemission is completely turned off and additional dephasing does not improve the efficiency considerably. Therefore, there is no need to have the controlled broadening turned on during stage 3 as mentioned earlier.

To have an idea about the validity of (10), we can compare it to a numerical solution of (5) and (6) using the method described in appendix B. In figure 2, we plot the efficiencies calculated by the numerical and the perturbative approach. The two approaches are seen to agree rather well for γ ≫ μ.

Zoom In Zoom Out Reset image size

The results of [24] show that ideal performance of stages 1 and 5 is achievable for sufficiently dense samples. The results obtained in this section demonstrate that also stages 2 and 4 can be made to work, thus demonstrating that an efficient memory is achievable if we can apply a strong reversible broadening of the atomic transition γ ≫ μ. To reach this result, we have for simplicity assumed that the incoming field is only incident in stage 1. As mentioned above we, however, allow for an incoming field also in stage 2 of the protocol in our numerical evaluation of the efficiency of the protocol. As we shall see below, allowing for this small tail of the incoming pulse to leak into stage 2 improves the efficiency of the protocol and allows for a much more rapid convergence with the broadening γ than predicted by the results of this section.

Having verified that the memory can, in principle, work with near 100% efficiency in the ideal limit of dense samples and large broadening, we now turn to a detailed numerical investigation of the performance of the protocol for real parameters. Specifically, we shall investigate the performance of the memory protocol for finite optical depth and finite broadening of the atomic transition.

The details of the numerical procedure used to solve the equations of motion are given in appendix B. We want to find the relation between the incoming and outgoing light fields. Since the underlying equations are all linear 'beam-splitter equations', which couple the annihilation operators of the effective harmonic oscillators, the relation between the incoming fields is also a beam splitter relation. As a result, the connection between the incoming fields in stages 1 and 2 and the outgoing fields in stages 4 and 5 can be written in the form

$$\begin{equation} E_{\rm out}({\tilde{t}}) =\int_0^{\mu\tau_{\rm R}} K_E({\tilde{t}},t') E_{\rm in}(\mu\tau_{\rm R}-t')\,{\rm d}t', \end{equation} \tag{ 11 }$$

where τ_R = τ_p + τ_d is the total duration of the read-in sequence and also the duration of the read-out sequence. The detailed expression for the matrix kernel K_E can be found in appendix B. In the end we are interested in the efficiency of the whole process, which is given by

$$\begin{equation} \eta=\int_0^{\mu\tau_{\rm R}} |E_{\rm out}({\tilde{t}})|^2 \,{\rm d}{\tilde{t}}, \end{equation} \tag{ 12 }$$

if we assume that the incoming pulse is normalized, i.e. $\int _0^{\mu \tau _{\mathrm { R}}} |E_{\mathrm { in}}({\tilde {t}})|^2\, \mathrm {d}{\tilde {t}} = 1$. Using expression (11) we can write the efficiency as

$$\begin{equation} \eta =\int_0^{\mu\tau_{\rm R}}\int_0^{\mu\tau_{\rm R}} E^*_{\rm in}(\mu\tau_{\rm R}-t') K_{\rm eff}(t',t'') E_{\rm in}(\mu\tau_{\rm R}-t'') \,{\rm d}t' \,{\rm d}t'', \end{equation} \tag{ 13 }$$

where

$$\begin{equation} K_{\rm eff}(t',t'') =\int_0^{\mu\tau_{\rm R}}K_E^*({\tilde{t}},t')K_E({\tilde{t}},t''){\rm d} {\tilde{t}}. \end{equation} \tag{ 14 }$$

In the numerical simulations, we discretize the time and use a quadrature rule [27] to set up a matrix K_eff(t,t') for the discrete times (nodes of the quadrature rule). With this matrix it is thus possible to obtain the memory efficiency for any vector containing the input fields at the discrete times.

For evaluating the performance of the memory, there are various approaches that one can take. From (13) we see that the memory efficiency depends on the shape of the incoming pulse. One approach to evaluating the performance is to look for the incoming mode that has the highest efficiency. Since the final expression for the efficiency (13) can be written in the form of a simple vector and matrix product, where the kernel matrix K_eff is self-adjoint, the maximal efficiency can be shown to be given by the largest eigenvalue of the matrix K_eff. Alternatively, the performance of the memory can be assessed by investigating how well the memory operates with a specially chosen mode that one may be interested in. Below we shall consider both these approaches. We emphasize, however, that there are many more methods of characterizing the performance. For instance, the approaches that we take here do not characterize the ability of the memory to store multiple modes [28] and we shall only briefly look at how input pulses get reshaped after storage. These properties may be important for practical applications.

To simulate the performance of the memory we consider a realistic situation where we have a collection of atoms with an intrinsic broadening of the optical transition described by the Gaussian distribution G₀ defined in (4) with width $\tilde {\gamma }_0=\gamma _0/\mu$. Since we are interested in a quantum memory operating only on two-level atoms, the width of this distribution inherently leads to a decay of the atomic polarization by an amount $\exp (-(\tau _{\mathrm { s}}/T_2)^2)$ with $T_2=d_0/(\sqrt {\pi }\mu )$ as shown in appendix A. At the same time, the width of the atomic line also determines the measurable optical depth d₀ of the ensemble before the controlled broadening is turned on which is given by $d_0=\sqrt {2\pi }\mu /\gamma _0$. In the investigations below, we are mainly interested in how the performance scales with the optical depth of the atomic ensembles. We therefore fix the storage time τ_s to be equal to $1/(2\gamma _0)=T_2/\sqrt {8}$ so that the maximal attainable storage efficiency is equal to $\exp (-\tau _{\mathrm { s}}^2\gamma _0^2)=\exp (-1/4)\approx 0.78$. With the storage time fixed relative to the dephasing time T₂, the investigations of the dependence of the optical depth below essentially correspond to the scaling one obtains when varying the number of atoms in the ensemble while keeping all the other parameters fixed. The allowed duration of the pulse to be stored is primarily determined by the duration τ_p of stage 1. For a memory to make sense, the duration of the pulse must be shorter than the storage time τ_p < τ_s. We therefore restrict ourselves to a duration τ_p = τ_s/4. Finally, we chose a constant duration τ_d of stages 2 and 4. As shown above, the duration of this period is not so important as long as it is long enough that after stage 2, all the different polarizations are dephased sufficiently so that negligible light gets out during stage 3. In the numerical simulations we are, however, constrained by having only a finite number of discrete frequencies. This means that there exists a finite time when the polarizations rephase again. Hence τ_d has to be chosen such that it is well below this rephasing time and we fix it at τ_d = 1/μ.

In figure 3, we show the maximal efficiency obtainable for a given optical depth and width of the broadening. In the figure we see that the efficiency rapidly increases when we apply the broadening and saturates when the width of the applied broadening reaches a value around γ ≳ 3μ. Once the applied broadening reaches this value it is sufficiently broad to rapidly dephase the polarization and thus rapidly turn off the reemission of the absorbed light after the broadening has been applied. Furthermore, we see that as we increase the optical depth of the ensemble d₀ the efficiency approaches the maximally allowed efficiency of η ≈ 0.78. The procedure proposed here thus allows for efficient quantum memory operation using only two-level atoms. In particular, by allowing the broadening to be turned on after the pulse is incident we are able to surpass the limit of η = 0.54 identified in [21] for a two-level quantum memory based on the standard transverse CRIB approach where the broadening is turned on before the pulse enters the medium (figure 1(a)). Furthermore for a sample with a long coherence time T₂ ≫ τ_s and a high optical depth d₀ ≫ 1, we can, in principle, come arbitrarily close to an efficiency of 100%.

Zoom In Zoom Out Reset image size

Some examples of the optimal modes leading to the maximal efficiency in figure 3 are shown in figure 4. For low values of the broadening, where the memory is inefficient, the optimal shape is rather flat. As we begin to increase the width of the applied broadening, the optimal shape changes character. After a certain value of γ which in this case is γ ≈ 1.4 μ, the shape of the modes for ${\tilde {t}} \leqslant \mu \tau _{\mathrm { p}}$ begins to resemble a Bessel-like function, which is the shape of the optimal modes identified in [24] (plotted in [26]) for the so-called fast memory regime, which corresponds to our stage 1. This resemblance reflects that the mechanisms in these approaches are highly similar. The only difference is that our dephasing mechanism which shuts off the reemission of the stored field is replaced in [24] by a π-pulse taking the excitation from |e〉 to an auxiliary state. After the broadening is applied at ${\tilde {t}}=\mu \tau _{\mathrm { p}}$, the optimal mode shape rapidly drops to zero on a time scale set by the width of the broadening. This reflects that 1/γ is the time scale needed to dephase the excitation and thus terminate the read-in sequence. For the retrieval we find that in the efficient regime (γ ≳ 1.4 μ with the parameters above) when the Bessel-like optimal modes are stored and then retrieved, the outgoing mode shapes are the same as the incoming mode shapes only time reversed and multiplied by minus the square root of efficiency (figure 1(a)).

Zoom In Zoom Out Reset image size

It is interesting to note that the efficiency in figure 3 approaches the stationary value much more rapidly when we increase the width γ/μ than shown in figure 2. The reason for this difference is that we allow for a small tail of the pulse to leak into stage 2, which was not included in section 4. In stage 2 the effective optical depth is given by $d=\sqrt {2\pi }\mu /\sqrt {\gamma ^2+\gamma _0^2}\approx \sqrt {2\pi }\mu /\gamma$ for $\gamma \gtrsim \gamma _0=\sqrt {2\pi }\mu / d_0$. Hence for the regime where the memory is efficient γ > μ the optical depth is below unity during stage 2. Regardless of this, the inclusion of a small optical field during this stage is still sufficient to alter the memory efficiency. This result emphasizes that the optical depth is not necessarily the correct physical parameter for characterizing the transient absorption of a pulse. The optical depth characterizes the fraction of the incident energy that is absorbed and not reemitted if the parameters of the memory are stationary. In our dynamic situation where the parameters are varying in time, the optical depth does not correctly characterize the transient absorption of the pulse. In fact, in our protocol most of the energy of the pulse is absorbed during stage 1 where the atomic line width is very narrow. Hence the pulse that we store is much broader in frequency than the atomic line width and most frequency components of the field see an optical depth much less than unity. This, however, only means that there is no stationary absorption and does not exclude the fact that the field is absorbed and reemitted several times during the passage through the memory (e.g. the reduced group velocity of a pulse traveling through a transparent medium, such as glass, can be understood as a consequence of constant absorption and reemission events for the light; in this case the resulting pulse delay can be substantial even if the optical depth is very low). In this sense the idea behind the current memory protocol is to interrupt the frequent absorption and reemission events by applying the broadening which stops and later resumes this absorption and reemission process.

Above we have focused on the maximal efficiency obtainable for the optimal mode shape. Since these optimal modes may not be available experimentally, we can also consider how the memory performs for a given predefined mode shape. We therefore also calculate the efficiency for Gaussian input modes of the form

$$\begin{equation} E_{\rm in}({\tilde{t}}, {\tilde{t}}_{\rm c}, {\tilde{t}}_{\rm w}) =\frac{1}{(2\pi{\tilde{t}}_{\rm w}^2)^{1/4}} \exp\left(-\frac{({\tilde{t}}-{\tilde{t}}_{\rm c})^2}{4{\tilde{t}}_{\rm w}^2}\right). \end{equation} \tag{ 15 }$$

This mode is characterized by two parameters, the center ${\tilde {t}}_{\mathrm { c}}$ and the width ${\tilde {t}}_{\mathrm { w}}$, and in figure 5 we show the efficiency as a function of these parameters. As we can see in the figure, there is always a well-defined maximum that we can find via numerical optimization. Furthermore in the limit where we expect the memory to work γ ≫ μ, the best performance is achieved by sending in a rather narrow pulse of duration ${\tilde {t}}_{\mathrm { w}}\sim 1/\mu$ right before broadening is applied. This result indicates that the best performance is achieved when the incident pulse resembles the sharply peaked optimal mode functions in figure 4.

Zoom In Zoom Out Reset image size

If we optimize ${\tilde {t}}_{\mathrm { w}}$, and ${\tilde {t}}_{\mathrm { c}}$, we can find the optimal Gaussian pulse for storage into the memory. The results of this optimization are shown in figure 6. Here we see that for the optimal Gaussian mode the efficiency also rapidly increases with the applied broadening and reaches a maximum for γ ≈ 3 μ. The efficiency then starts to decrease, which is a consequence of the Gaussian mode shape not resembling the optimal mode shape. The efficiency also increases with increasing optical depth of the ensemble d₀ although much slower than for the optimal mode in figure 3. Unfortunately limited numerical resources prevent us from increasing the optical depth even further, but we believe that it will eventually converge toward the maximal attainable efficiency of $\eta \leqslant \exp (-2 (\tau _{\mathrm { s}}/T_2)^2)\approx 0.78$. Nevertheless our results are still able to surpass the limit of η = 54%·0.78 = 42% for the two-level memory protocol presented in [21] when we take into account the dephasing during the memory time.

Zoom In Zoom Out Reset image size

Before concluding, we briefly comment on the experimental requirements for realizing our scheme. The main requirements are that the ensemble should have a large optical depth and a long coherence time. Taking for instance a line width γ₀ = (2π)140 kHz [18], we find that the setting used in the numerical simulations corresponds to a storage time τ_s = 1/(2γ₀) ∼ 0.6 μs. In addition, we should be able to apply a sufficiently strong broadening to rapidly dephase the atomic transition. Taking an optimistic value of d₀ = 100, this corresponds to an induced width $\gamma \gtrsim 3\,\mu =3d_0\gamma _0/\sqrt {2\pi } \sim (2\pi )$17 MHz, which is about an order of magnitude larger than achieved in [18]. Note that the bandwidth μ ∼ (2π)6 MHz of the incoming pulse may be problematic for implementations based on rare-earth dopants with a small hyperfine splitting, e.g. europium and praseodymium, since the width is comparable to the achievable spectral hole burning width, but it is less of a problem for dopants with a larger hyperfine splitting, e.g. terbium and holmium [29]. In addition, the scheme introduces the requirement that the field should be switched on a time scale less than 1/μ ∼ 30 ns. For comparison [18] finds that the switching times in the experiment are less than 150 ns. Achieving the efficient memory operation described here thus seems challenging but not out of reach. In particular, the requirements for the switching time and width of the applied broadening can be relaxed if more narrow line shapes can be achieved.

We have proposed a method to make an efficient quantum memory for light based on two-level atoms. In our proposal the memory operation is controlled by applying and later reversing an external field that broadens the atomic transition. For sufficiently strong broadening and sufficiently high optical depth before the broadening is applied, we have shown that the memory operation can be limited only by the intrinsic decoherence (the T₂ time) of the atomic transition. Most importantly, in contrast to most protocols the proposed memory can be efficient even without employing any optical control field. Since there is no need for an additional laser field, we hope that this memory protocol may be simpler to implement in practice.

Our protocol is not the first proposal for a quantum memory based on two-level atoms. In [12, 13] an efficient memory was proposed for the so-called longitudinal CRIB where a field gradient is applied such that the shift of the atomic transition depends on the position along the propagation direction. In contrast, our protocol is based on transverse CRIB where the shift does not have this dependence. For this setting the only previous protocol [21] based on two-level systems had a maximal efficiency of 54%. For our protocol there is no such limit and the efficiency can approach 100% for sufficiently high optical depth. The efficiency is, however, dependent on the mode shape and converges to the ideal limit much more rapidly for the optimal mode shape than if we constrain it to be e.g. a Gaussian. It would be interesting to study in more detail how the protocol developed here compares to the other memory protocols based on two-level systems. The previous work, however, had a different scope than what we consider here and only studied the dependence on the optical depth after the broadening was applied, and the effects of the intrinsic broadening were not included in the final results for the efficiency. Within our framework these previous studies thus correspond to the limit d₀ → ∞ and cannot directly be compared to our results.

It has previously been argued that adding and reversing (transverse) broadening during the storage process generally reduces the memory efficiency for the optimal modes as compared to schemes that do not have this broadening [26]. In our scheme the broadening is absent during the period where most of the light is read into the memory. We therefore believe that for a given ensemble our proposed protocol will have the highest possible efficiency for the optimal mode. On the other hand, it has also been shown that the previous CRIB protocols [12, 21] have a much larger multimode capacity than protocols such as this one, where the broadening is not present during the storage [28]. This limited multimode capacity of our scheme is reflected in our results for Gaussian pulse shapes, where the efficiency becomes limited by the ability of the Gaussian to resemble the relatively few modes that are stored with high efficiency. In fact for the read-out process, we find that also for the case of an incoming Gaussian pulse, the outgoing mode shape resembles the outgoing mode shape for the optimal mode. Again this is most likely a consequence of the limited multimode efficiency: only a few storage modes are excited by the Gaussian, and these modes are later retrieved into outgoing modes that contain little information about the incoming pulse shape. This reshaping of the mode may be detrimental for some applications of the quantum memory, but for other applications the mode shape may be less important. For instance, for quantum repeaters [30, 31] we are interested in interfering with the output from two different quantum memories. Since both memories in this case will emit similar mode shapes, the reshaping of the mode is of minor consequence. For concrete applications of quantum memories, a more detailed study will have to be performed to determine whether the increased efficiency of the present protocol outweighs the drawback of the limited multimode capacity.

We thank A Grodecka-Grad and N Sangouard for helpful discussions. This work was supported by the Danish National Research Foundation and by the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC grant agreement no. 306576.

In this section we relate the parameters of the model to physically measurable quantities. Specifically we shall express everything in terms of the optical depth in the absence of controlled broadening d₀ and the coherence time of the polarization of the atoms T₂. To do this we consider light propagating through the atomic ensemble without the controlled broadening being present. This situation can be described by (B.1) and (B.2). If we formally integrate (B.2) under the initial condition $P({\tilde {z}}, {\tilde {t}}=-\infty ,\tilde {\Delta }_0)=0$, we obtain

$$\begin{eqnarray*} &&P({\tilde{z}},{\tilde{t}},\tilde{\Delta}_0) ={\rm i}\int_{-\infty}^{\tilde{t}}\exp(-\mathrm{i}\Delta_0({\tilde{t}}-t'))E({\tilde{z}}, t'){\rm d} t'. \end{eqnarray*}$$

Taking the Fourier transform (${\tilde {t}}\rightarrow \omega$) this equation becomes

$$\begin{equation} \tilde{P}({\tilde{z}},\omega,\tilde{\Delta}_0) ={\rm i}\tilde{E}({\tilde{z}}, \omega) \left(\pi\delta(\omega-\tilde{\Delta}_0)-{\rm PV}\frac{{\rm i}}{\omega-\tilde{\Delta}_0}\right), \end{equation} \tag{ A.1 }$$

where PV reminds us that we need to take the Cauchy principal value when integrating. Inserting (A.1) into the Fourier transform of (B.1) results in

$$\begin{eqnarray*} &&\frac{\partial}{\partial {\tilde{z}}}\skew3\tilde{E}({\tilde{z}},\omega) =-\skew3\tilde{E}({\tilde{z}},\omega)\left(\pi G_0(\omega) -{\rm PV}\int_{-\infty}^{\infty} G_0(\tilde{\Delta}_0')\frac{{\rm i}}{\omega-\tilde{\Delta}_0'}{\rm d} \tilde{\Delta}_0'\right). \end{eqnarray*}$$

From this we see that

$$\begin{eqnarray*} &&|\skew3\tilde{E}({\tilde{z}}=1,\omega)|^2 =\exp\left(-2\pi G_0(\omega))\right |\skew3\tilde{E}({\tilde{z}}=0,\omega)|^2. \end{eqnarray*}$$

From this expression we define the optical depth on resonance as d₀ = 2πG₀(0). If the controlled broadening is present, this amounts to having a distribution of the sum of the controlled and intrinsic broadening G_sum instead of G₀. Hence we have d = 2πG_sum(0). If both G₀ and the controlled broadening distribution G are Gaussian with widths γ₀/μ and γ/μ, respectively, then G_sum is also Gaussian with width $\gamma _{\mathrm { sum}}/\mu =\sqrt {\gamma _0^2+\gamma ^2}/\mu$. In the limit γ ≫ γ₀ we simply have γ_sum ≈ γ, so that

$$\begin{eqnarray*} &&d_0 =\sqrt{2\pi}\frac{\mu}{\gamma_0} \quad {\rm and} \quad d =\sqrt{2\pi}\frac{\mu}{\sqrt{\gamma_0^2+\gamma^2}} \approx\sqrt{2\pi}\frac{\mu}{\gamma}. \end{eqnarray*}$$

To derive the temporal decay of the polarization, we can assume that P in (3) is independent of $\tilde {\Delta }_0$ initially, i.e. $P({\tilde {z}},{\tilde {t}}=0,\Delta _0)=P_0({\tilde {z}})$. If no electric field is present, the evolution is given by (B.2) and the total polarization becomes

$$\begin{eqnarray*} &&P_{\rm total}({\tilde{z}},{\tilde{t}}) =N\int_{-\infty}^\infty G_0(\tilde{\Delta}) P({\tilde{z}}, {\tilde{t}},\Delta_0) {\rm d}\tilde{\Delta} =N P_0({\tilde{z}})\int_{-\infty}^\infty G_0(\tilde{\Delta}) {\rm e}^{-{\rm i}\tilde{\Delta}_0{\tilde{t}}} {\rm d}\tilde{\Delta}. \end{eqnarray*}$$

Hence for a Gaussian G₀ given by (4) we have

$$\begin{eqnarray*} &&P_{\rm total}({\tilde{z}},{\tilde{t}}) =P_{\rm total}({\tilde{z}},{\tilde{t}}=0)\exp\left(-\frac{{\tilde{t}}^2}{\mu^2 T_2^2}\right) \end{eqnarray*}$$

with $T_2=\sqrt {2}/\gamma _0=d_0/(\sqrt {\pi }\mu )$. This expression allows us to relate the bandwidth of the memory μ to the measurable quantities d₀ and T₂.

In this appendix we give the details of how we do the numerical simulations of the memory protocol. In short we Laplace transform the equations of motion in space (${\tilde {z}}\rightarrow u$) and rewrite them so that we obtain an equation only involving the polarization. Then we discretize the broadening variables $\tilde {\Delta }_0$ and $\tilde {\Delta }$ such that the integrals become sums and convolutions become matrix products. The resulting vector equations have simple solutions expressed in terms of the matrix exponential. In the end we can find the Laplace transform of the electric field using the discretized and Laplace transformed versions of (1) and (B.1). Then the electric field in real space can be evaluated by applying the numerical inverse Laplace transform that amounts to transforming the integration kernels. Detailed calculations for each stage are given below.

B.1. Stage 1

In stage 1 the controlled broadening is not present so that σ does not depend on $\tilde {\Delta }$. Hence using the definition (3), we see that $\sigma ({\tilde {z}},{\tilde {t}},\tilde {\Delta }_0,\tilde {\Delta })=P({\tilde {z}},{\tilde {t}},\tilde {\Delta }_0)$. Also the term $-{\mathrm { i}}\tilde {\Delta }\sigma ({\tilde {z}},{\tilde {t}},\tilde {\Delta }_0,\tilde {\Delta })$ is not present on the right-hand side of (2) so that we can write it together with (1) as

$$\begin{equation} \frac{\partial}{\partial {\tilde{z}}}E({\tilde{z}},{\tilde{t}})={\rm i} \int_{-\infty}^{\infty} G_0(\tilde{\Delta}_0) P({\tilde{z}},{\tilde{t}},\tilde{\Delta}_0) \,{\rm d}\tilde{\Delta}_0, \end{equation} \tag{ B.1 }$$

$$\begin{equation} \frac{\partial}{\partial {\tilde{t}}}P({\tilde{z}},{\tilde{t}},\tilde{\Delta}_0) =-{\rm i}\tilde{\Delta}_0 P({\tilde{z}},{\tilde{t}},\tilde{\Delta}_0) +{\rm i} E({\tilde{z}},{\tilde{t}}). \end{equation} \tag{ B.2 }$$

Now we Laplace transform these equations in space (${\tilde {z}}\rightarrow u$). Combining the resulting equations and using the definition $E_{\mathrm { in}}({\tilde {t}})=E({\tilde {z}}=0,{\tilde {t}})$ we obtain

$$\begin{equation}\fl \frac{\partial}{\partial {\tilde{t}}}\skew3\bar{P}(u,{\tilde{t}},\tilde{\Delta}_0) =-{\rm i}\tilde{\Delta}_0 \skew3\bar{P}(u,{\tilde{t}},\tilde{\Delta}_0) -\frac{1}{u}\int_{-\infty}^\infty G_0(\tilde{\Delta}_0') \skew3\bar{P}(u,{\tilde{t}},\tilde{\Delta}_0'){\rm d}\tilde{\Delta}_0' +\frac{{\rm i}}{u}E_{\rm in}({\tilde{t}}). \end{equation} \tag{ B.3 }$$

We choose K discrete values of $\tilde {\Delta }_0$ that we call $\tilde {\Delta }_{01}, \ldots , \tilde {\Delta }_{0K}$. For simplicity, we choose them so that they have a constant step $\tilde {\Delta }_{0\delta }$. After discretization (B.3) becomes

$$\begin{equation}\fl \frac{\partial}{\partial {\tilde{t}}}\skew3\bar{P}(u,{\tilde{t}},\tilde{\Delta}_{0j}) =-{\rm i}\tilde{\Delta}_{0j} \skew3\bar{P}(u,{\tilde{t}},\tilde{\Delta}_{0j}) -\frac{1}{u}\sum_{k=1}^K G_0(\tilde{\Delta}_{0k})\tilde{\Delta}_{0\delta}\skew3\bar{P}(u,{\tilde{t}},\tilde{\Delta}_{0k}) +\frac{{\rm i}}{u}E_{\rm in}({\tilde{t}}). \end{equation} \tag{ B.4 }$$

To simplify the notation we wish to write the discretized equations in vector form. To do this we define vectors Δ₀, $\mathbf {P}(u,{\tilde {t}})$ and g₀ with elements given by $({{\boldsymbol {\Delta}} }_0)_j=\tilde {\Delta }_{0j}$, $(\mathbf {P}(u,{\tilde {t}}))_j=\skew3\bar {P}(u,{\tilde {t}},\tilde {\Delta }_{0j})$ and $(\mathbf {g}_0)_j=\tilde {\Delta }_{0\delta }G_0(\tilde {\Delta }_{0j})$ for 1 ⩽ j ⩽ K. Also for any vector v we define a matrix D(v) with elements D(v)_jj = v_j for all j and D(v)_jk = 0 for j ≠ k, i.e. a diagonal matrix with v as its diagonal. If we further let h^(K) be a K-dimensional vector with constant elements h_j = 1 for all j, we can write (B.4) as a vector equation

$$\begin{equation} \frac{\partial}{\partial {\tilde{t}}}\mathbf{P}(u,{\tilde{t}}) =-{\rm i} D({{\boldsymbol {\Delta}}}_0) \mathbf{P}(u,{\tilde{t}}) -\frac{1}{u}\mathbf{h}^{(K)} \mathbf{g}_0^T\mathbf{P}(u,{\tilde{t}}) +\frac{{\rm i}}{u}E_{\rm in}({\tilde{t}})\mathbf{h}^{(K)}. \end{equation} \tag{ B.5 }$$

Note that h^(K)g^T₀ is a matrix with each row equal to the vector g₀. Defining

$$\begin{eqnarray*} &&M_1(u)=-{\rm i} D({{\boldsymbol {\Delta}}}_0)-\frac{1}{u}\mathbf{h}^{(K)} \mathbf{g}_0^T, \end{eqnarray*}$$

we can write (B.5) as

$$\begin{equation} \frac{\partial}{\partial {\tilde{t}}}\mathbf{P}(u,{\tilde{t}}) =M_1(u)\mathbf{P}(u,{\tilde{t}}) +\frac{{\rm i}}{u}E_{\rm in}({\tilde{t}})\mathbf{h}^{(K)}. \end{equation} \tag{ B.6 }$$

In stage 1 the initial condition is that the memory is empty, i.e $P({\tilde {z}},{\tilde {t}}=0)=0$ or equivalently $\mathbf {P}(u,{\tilde {t}}=0)=0$. The solution to (B.6) can then be expressed using the matrix exponential

$$\begin{equation} \mathbf{P}(u,{\tilde{t}})={\rm i}\int_0^{\tilde{t}} \frac{1}{u}\exp(M_1(u)t') \mathbf{h}^{(K)} E_{\rm in}({\tilde{t}}-t')\,{\rm d} t'. \end{equation} \tag{ B.7 }$$

B.2. Stage 2

The memory protocol works by having most of the light absorbed during stage 1. Nevertheless we still allow for some input light also during stage 2. Hence stage 2 is described by (1) and (2) with the initial conditions

$$\begin{eqnarray*} & &E({\tilde{z}}=0,{\tilde{t}})=E_{\rm in}({\tilde{t}}+\mu\tau_{\rm p}),\\ & &\sigma({\tilde{z}},{\tilde{t}}=0,\tilde{\Delta}_0,\tilde{\Delta}) =\sigma^{(1)}({\tilde{z}},\tilde{\Delta}_0,\tilde{\Delta}) =P^{(1)}({\tilde{z}},\tilde{\Delta}_0) \end{eqnarray*}$$

where σ⁽¹⁾ and P⁽¹⁾ are two ways of denoting the polarization at the end of stage 1; and E_in is the same function as defined for stage 1. Here and in the beginning of all the subsequent stages, the time ${\tilde {t}}$ represents the time since the beginning of the current stage. This resetting of time is convenient since the equations of motion change from one stage to another. Since we want to consider E_in as being an input pulse for both stage 1 and stage 2, we shift its argument by the duration of stage 1, when we solve the equations of motion in stage 2.

Taking the Laplace transform of (1) and (2) and combining the resulting equations as for stage 1, we obtain

$$\begin{eqnarray}&&\fl \frac{\partial}{\partial {\tilde{t}}}\bar{\sigma}(u,{\tilde{t}},\tilde{\Delta}_0,\tilde{\Delta}) \!=\!-{\rm i}(\tilde{\Delta}_0\!+\!\tilde{\Delta}) \bar{\sigma}(u,{\tilde{t}},\tilde{\Delta}_0,\tilde{\Delta})\!-\!\frac{1}{u}\!\int_{-\infty}^\infty \!\int_{-\infty}^\infty\! G_0(\tilde{\Delta}_0') G(\tilde{\Delta}') \bar{\sigma}(u,{\tilde{t}},\tilde{\Delta}_0',\tilde{\Delta})\,{\rm d}\tilde{\Delta}_0'\,{\rm d}\tilde{\Delta}'\nonumber\\ &&+\frac{{\rm i}}{u}E_{\rm in}({\tilde{t}}+\mu\tau_{\rm p}). \end{eqnarray} \tag{ B.8 }$$

Now we have to split both detunings $\tilde {\Delta }_0$ and $\tilde {\Delta }$ into respectively K discrete values $\tilde {\Delta }_{01}, \ldots , \tilde {\Delta }_{0K}$ and N discrete values $\tilde {\Delta }_{1}, \ldots , \tilde {\Delta }_{N}$. Again we assume that they have constant steps $\tilde {\Delta }_{0\delta }$ and $\tilde {\Delta }_\delta$, respectively. Defining $\bar {\sigma }_{j,k}(u,{\tilde {t}})=\bar {\sigma }(u,{\tilde {t}},\tilde {\Delta }_{0j},\tilde {\Delta }_k)$, the discretized version of (B.8) can be written as

$$\begin{eqnarray}&&\fl \frac{\partial}{\partial {\tilde{t}}}\bar{\sigma}_{j,k}(u,{\tilde{t}}) =-{\rm i}(\tilde{\Delta}_{0j}+\tilde{\Delta}_k) \bar{\sigma}_{j,k}(u,{\tilde{t}})-\frac{1}{u}\sum_{j'=1}^K \sum_{k'=1}^N G_0(\tilde{\Delta}_{0j'}) G(\tilde{\Delta}_{k'})\tilde{\Delta}_{0\delta}\tilde{\Delta}_\delta \bar{\sigma}_{j',k'}(u,{\tilde{t}})\nonumber\\ &&+\frac{{\rm i}}{u}E_{\rm in}({\tilde{t}}+\mu\tau_{\rm p}). \end{eqnarray} \tag{ B.9 }$$

We see that we can write this expression in vector notation as before if we combine the indices j and k into one. Hence we define vectors Δ_±, ${\boldsymbol {\sigma }}(u,{\tilde {t}})$ and g with elements $({{\boldsymbol {\Delta}}}_\pm )_{(j-1)N+k}=\tilde {\Delta }_{0j}\pm \tilde {\Delta }_k$, $({\boldsymbol {\sigma }}(u,{\tilde {t}}))_{(j-1)N+k}=\bar {\sigma }_{j,k}(u,{\tilde {t}})$ and $(\mathbf {g})_{(j-1)N+k}\!=\!\tilde {\Delta }_{0\delta }\tilde {\Delta }_\delta G _0(\tilde {\Delta }_{0j})G(\tilde {\Delta }_k)$ for 1 ⩽ j ⩽ K, 1 ⩽ k ⩽ N. With these definitions (B.9) becomes

$$\begin{equation}\fl \frac{\partial}{\partial {\tilde{t}}}{\boldsymbol{\sigma}}(u,{\tilde{t}}) =-{\rm i} D({{\boldsymbol {\Delta}}}_+) {\boldsymbol{\sigma}}(u,{\tilde{t}}) -\frac{1}{u}\mathbf{h}^{(KN)}\mathbf{g}^T {\boldsymbol{\sigma}}(u,{\tilde{t}}) +\frac{{\rm i}}{u}E_{\rm in}({\tilde{t}}+\mu\tau_{\rm p})\mathbf{h}^{(KN)}. \end{equation} \tag{ B.10 }$$

Defining M₂(u) = −iD(Δ₊) − h^(KN)g^T/u, we can write (B.10) as

$$\begin{equation} \frac{\partial}{\partial {\tilde{t}}}{\boldsymbol{\sigma}}(u,{\tilde{t}}) =M_2(u){\boldsymbol{\sigma}}(u,{\tilde{t}}) +\frac{{\rm i}}{u}E_{\rm in}({\tilde{t}}+\mu\tau_{\rm p})\mathbf{h}^{(KN)} \end{equation} \tag{ B.11 }$$

and the solution at the end of stage 2 is given by

$$\begin{equation}\fl {\boldsymbol{\sigma}}^{(2)}(u) =\frac{{\rm i}}{u}\int_0^{\mu\tau_{\rm d}} \exp(M_2(u)t)\mathbf{h}^{(KN)} E_{\rm in}(\mu\tau_{\rm R}-t)\,{\rm d}t +\exp(M_2(u)\mu\tau_{\rm d}){\boldsymbol{\sigma}}^{(1)}(u) \end{equation} \tag{ B.12 }$$

with τ_R = τ_p + τ_d as defined earlier. Here we have used the initial condition σ⁽¹⁾(u) that in the vector notation has elements (σ⁽¹⁾(u))_(j−1)N+k = (P(u,μτ_p))_j for 1 ⩽ j ⩽ K and 1 ⩽ k ⩽ N with $\mathbf {P}(u,{\tilde {t}})$ given by (B.7).

B.3. Stage 3

To describe stage 3 we could in principle have used the same equations as for stage 1 but then we would have lost the information about how the individual frequency components of the controlled broadening evolved. Hence doing this would not have allowed us to describe the rephasing in stage 4. Instead we use an equation similar to (B.9) and set $\tilde {\Delta }_k=0$. Here we do not have an incoming field so that the evolution is given by

$$\begin{equation}\fl \frac{\partial}{\partial {\tilde{t}}}\bar{\sigma}_{j,k}(u,{\tilde{t}}) =-{\rm i}\tilde{\Delta}_{0j} \bar{\sigma}_{j,k}(u,{\tilde{t}}) -\frac{1}{u}\sum_{j'=1}^K \sum_{k'=1}^N G_0(\tilde{\Delta}_{0j'}) G(\tilde{\Delta}_{k'})\tilde{\Delta}_{0\delta}\tilde{\Delta}_\delta \bar{\sigma}_{j',k'}(u,{\tilde{t}}). \end{equation} \tag{ B.13 }$$

We define a vector Δ with elements $({\boldsymbol{\Delta }})_{(j-1)N+k}=\tilde {\Delta }_{0j}$ for 1 ⩽ j ⩽ K, 1 ⩽ k ⩽ N and the corresponding matrix M₃(u) = −iD(Δ) − h^(KN)g^T/u. Then we can write (B.13) as

$$\begin{equation} \frac{\partial}{\partial {\tilde{t}}}{\boldsymbol{\sigma}}(u,{\tilde{t}}) =M_3(u){\boldsymbol{\sigma}}(u,{\tilde{t}}) \end{equation} \tag{ B.14 }$$

so that the solution at the end of stage 3 is given by

$$\begin{equation} {\boldsymbol{\sigma}}^{(3)}(u) =\exp(M_3(u)\mu\tau_{\rm d}){\boldsymbol{\sigma}}^{(2)}(u). \end{equation} \tag{ B.15 }$$

B.4. Stage 4

For stage 4 we define M₄(u) = −iD(Δ₋) − h^(KN)g^T/u and the evolution of σ is then described by an equation of the same form as (B.14) with the solution for σ given by

$$\begin{equation} {\boldsymbol{\sigma}}(u,{\tilde{t}})=\exp(M_4(u){\tilde{t}}){\boldsymbol{\sigma}}^{(3)}(u). \end{equation} \tag{ B.16 }$$

From this expression we can find the electric field using the Laplace transformed and discretized version of (1) which can be written as

$$\begin{equation} \skew3\bar{E}(u,{\tilde{t}}) =\frac{{\rm i}}{u}\mathbf{g}^T{\boldsymbol{\sigma}}(u,{\tilde{t}}). \end{equation} \tag{ B.17 }$$

We define the matrix

$$\begin{eqnarray*} &&J^{(KN\times KN)}(u)=\exp(M_3(u)\mu\tau_{\rm s})\exp(M_2(u)\mu\tau_{\rm d}) \end{eqnarray*}$$

and the kernel

$$\begin{eqnarray*} &&\bar{k}_1(u,{\tilde{t}},t')=-\frac{1}{u^2}\mathbf{g}^T \exp(M_4(u){\tilde{t}})\exp(M_3(u)\mu\tau_{\rm s})\exp(M_2(u)t') \mathbf{h}^{(KN)}. \end{eqnarray*}$$

Then we can use (B.12), (B.15)–(B.17) to write the field as

$$\begin{equation}\fl \skew3\bar{E}(u,{\tilde{t}}) =\int_0^{\mu\tau_{\rm d}} \bar{k}_1(u,{\tilde{t}},t') E_{\rm in}(\mu\tau_{\rm R}-t')\,{\rm d}t' -\displaystyle\frac{1}{u^2}\mathbf{g}^T\exp(M_4(u){\tilde{t}})J^{(KN\times KN)}(u) {\boldsymbol{\sigma}}^{(1)}(u). \end{equation} \tag{ B.18 }$$

Here the first term describes how the field incident in stage 2 is read-out during stage 4 and the last term describes the read-out of the field incident during stage 1.

Before we can write the second term in this expression as something involving the input field, we need to express σ⁽¹⁾(u) in terms of P⁽¹⁾(u) so that we can use (B.7). In vector notation, P⁽¹⁾(u) is a K-dimensional vector while σ⁽¹⁾(u) is a KN-dimensional vector as defined in the end of (B.2). This definition means that every element in block number j of σ⁽¹⁾(u) with length N has the same value (P⁽¹⁾(u))_j. Hence when the matrix J^(KN×KN)(u) is multiplied with σ⁽¹⁾(u), each element of the resulting vector is given by

$$\begin{eqnarray*} (J^{(KN\times KN)}(u){\boldsymbol{\sigma}}^{(1)}(u))_j &=&\sum_{k=1}^{KN}(J^{(KN\times KN)}(u))_{j,k}({\boldsymbol{\sigma}}^{(1)}(u))_k\\ &=&\sum_{j'=1}^{K}\sum_{k'=1}^{N}(J^{(KN\times KN)}(u))_{j,(j'-1)N+k'} ({\boldsymbol{\sigma}}^{(1)}(u))_{(j'-1)N+k'}. \end{eqnarray*}$$

In the last line, (σ⁽¹⁾(u))_(j'−1)N+k' = (P⁽¹⁾(u))_j' so that it can be taken out of the k' summation. Hence if we define a matrix J^(KN×K)(u) with elements

$$\begin{eqnarray*} &&(J^{(KN\times K)}(u))_{j,k} =\sum_{k'=1}^N (J^{(KN\times KN)}(u))_{j,(k-1)N+k'} \end{eqnarray*}$$

then J^(KN×KN)(u)σ⁽¹⁾(u) = J^(KN×K)(u)P⁽¹⁾(u). If we introduce another kernel

$$\begin{eqnarray*} &&\bar{k}_2(u,{\tilde{t}},t')=-\frac{1}{u^2}\mathbf{g}^T \exp(M_4(u){\tilde{t}})J^{(KN\times K)}(u)\exp(M_1(u)t') \mathbf{h}^{(K)} \end{eqnarray*}$$

and take the inverse Laplace transform ($u\rightarrow {\tilde {z}}=1$), we can finally express the output field during stage 4 in terms of the input field

$$\begin{equation}\fl E_{\rm out}({\tilde{t}}) =\int_0^{\mu\tau_{\rm d}} k_1({\tilde{t}},t') E_{\rm in}(\mu\tau_{\rm R}-t')\,{\rm d}t' +\int_0^{\mu\tau_{\rm p}} k_2({\tilde{t}},t') E_{\rm in}(\mu\tau_{\rm p}-t')\,{\rm d} t', \end{equation} \tag{ B.19 }$$

where the kernels k₁ and k₂ are inverse Laplace transforms of $\bar {k}_1$ and $\bar {k}_2$ evaluated at ${\tilde {z}}=1$.

B.5. Stage 5

In stage 5 the initial condition for σ is given by (B.16) with ${\tilde {t}}=\mu \tau _{\mathrm { d}}$. Defining

$$\begin{eqnarray*} &&L^{(KN\times KN)}(u) = \exp(M_4(u)\mu\tau_{\rm d})\exp(M_3(u)\mu\tau_{\rm s}),\\ &&B^{(KN\times KN)}(u) = \exp(M_4(u)\mu\tau_{\rm d})\exp(M_3(u)\mu\tau_{\rm s})\exp(M_2(u)\mu\tau_{\rm d}), \end{eqnarray*}$$

we can write the initial condition as

$$\begin{eqnarray*}&&\fl {\boldsymbol{\sigma}}^{(4)}(u) =\frac{{\rm i}}{u}\int_0^{\mu\tau_{\rm d}} L^{(KN\times KN)}(u)\exp(M_2(u)t') \mathbf{h}^{(KN)} E_{\rm in}(\mu\tau_{\rm R}-t')\,{\rm d}t'+B^{(KN\times KN)}(u){\boldsymbol{\sigma}}^{(1)}(u). \end{eqnarray*}$$

The electric field can be found using the Laplace transformed and discretized version of (B.1):

$$\begin{eqnarray*} &&\skew3\bar{E}(u,{\tilde{t}}) =\frac{{\rm i}}{u}\mathbf{g}_0^T\mathbf{P}(u,{\tilde{t}}), \end{eqnarray*}$$

where $\mathbf {P}(u,{\tilde {t}})=\exp (M_1(u){\tilde {t}})\mathbf {P}^{(4)}(u)$ is the solution to (B.6) with E_in = 0. The initial condition P⁽⁴⁾(u) is defined by (3). In vector notation this translates to blockwise summation of σ⁽⁴⁾(u) with respect to the index of the controlled broadening weighted by the appropriate value of the controlled broadening distribution. We can apply this blockwise summation to the matrices L^(KN×KN)(u) and B^(KN×KN)(u) directly. Hence we define L^(K×KN)(u) and B^(K×K)(u) with

$$\begin{eqnarray*} &&(L^{(K\times KN)}(u))_{j,k}=\sum_{j'=1}^N G(\tilde{\Delta}_{j'})\tilde{\Delta}_\delta (L^{(KN\times KN)}(u))_{(j-1)N+j',k},\\ &&(B^{(K\times K)}(u))_{j,k}=\sum_{j'=1}^N \sum_{k'=1}^N G(\tilde{\Delta}_{j'})\tilde{\Delta}_\delta (B^{(KN\times KN)}(u))_{(j-1)N+j',(k-1)N+k',} \end{eqnarray*}$$

where we also apply the same reduction for the matrix B^(KN×KN)(u) as for J^(KN×KN)(u) previously. We can now use kernels

$$\begin{eqnarray*} &&\bar{k}_3(u,{\tilde{t}},t')=-\frac{1}{u^2}\mathbf{g}_0^T \exp(M_1(u){\tilde{t}})L^{(K\times KN)}(u)\exp(M_2(u)t') \mathbf{h}^{(KN)},\\ &&\bar{k}_4(u,{\tilde{t}},t')=-\frac{1}{u^2}\mathbf{g}_0^T \exp(M_1(u){\tilde{t}})B^{(K\times K)}(u)\exp(M_1(u)t') \mathbf{h}^{(K)} \end{eqnarray*}$$

to express the output electric field in terms of the input field E_in. Taking the inverse Laplace transform ($u\rightarrow {\tilde {z}}=1$), the output field is

$$\begin{equation}\fl E_{\rm out}({\tilde{t}}) =\int_0^{\mu\tau_{\rm d}} k_3({\tilde{t}},t') E_{\rm in}(\mu\tau_{\rm R}-t')\,{\rm d}t' +\int_0^{\mu\tau_{\rm p}} k_4({\tilde{t}},t') E_{\rm in}(\mu\tau_{\rm p}-t')\,{\rm d}t', \end{equation} \tag{ B.20 }$$

where the kernels k₃ and k₄ are inverse Laplace transforms of $\bar {k}_3$ and $\bar {k}_4$ evaluated at ${\tilde {z}}=1$.

B.6. Concluding remarks

With the above results we have found the sought relation between the input and output field. The final expressions in (B.19) and (B.18) give the output field during stages 4 and 5, respectively. They can be combined into a single expression (11) if we define

$$\begin{eqnarray*} K_E({\tilde{t}},t')=\left\{\begin{array}{@{}ll@{}} k_1({\tilde{t}},t') & {\rm for}\quad {\tilde{t}}\leqslant\mu\tau_{\rm d},t'\leqslant\mu\tau_{\rm d},\\ k_2({\tilde{t}},t'-\mu\tau_{\rm d}) & {\rm for}\quad {\tilde{t}}\leqslant\mu\tau_{\rm d},t'>\mu\tau_{\rm d},\\ k_3({\tilde{t}}-\mu\tau_{\rm d},t') & {\rm for}\quad {\tilde{t}}>\mu\tau_{\rm d},t'\leqslant\mu\tau_{\rm d},\\ k_4({\tilde{t}}-\mu\tau_{\rm d},t'-\mu\tau_{\rm d}) & {\rm for}\quad {\tilde{t}}>\mu\tau_{\rm d},t'>\mu\tau_{\rm d}. \end{array}\right. \end{eqnarray*}$$

Here the shifts of the arguments of the kernels are performed because we go from considering E_out to be output fields only in stage 4 or only in stage 5 (where time was reset at the beginning of each stage) to having a single function E_out giving outgoing fields in both stages.

For inversion of the Laplace transform, we use numerical integration along the Talbot contour [32]. Since this method has exponential convergence, only a few discrete Laplace moments u are needed. This makes it possible to compute and store certain intermediate results for each u. For example the kernels k₁, k₂, k₃ and k₄ can be computed much faster if the matrices M₁(u), M₂(u) and M₄(u) are diagonalized. Hence we compute and store the eigenvalues and eigenvectors of those matrices, so that matrix exponentials can be evaluated by exponentiating the eigenvalues instead of the whole matrix.