Optimal photon generation from spontaneous Raman processes in cold atoms

We first review a derivation in [15] giving the dependency of the efficiency of the retrieval process on the relevant quantities in the physical setup. This yields an expression for the retrieval efficiency, that depends only on the shape of the spin wave from which the retrieval is performed, and on the optical depth of the relevant transition.

We emphasise that the work in [15] focuses on absorptive memory protocols where a field is first absorbed in an atomic medium, creating a spin wave that can be read out later to re-emit the field in a well defined spatio-temporal mode. To justify the relation to [15], in our proposal the spin wave creation is instead heralded by the detection of the write photon field, but the readout process is analogous, allowing us to make use of [15] to deduce the spin wave shapes that maximise the retrieval efficiency.

We consider a three-level atomic system in a Λ-configuration (see figure 1(a)) with spin excitations present in the form of $| g\rangle -| s\rangle$ coherences. In the situation where almost all the atoms remain in $| g\rangle$ and in a rotating frame, the backward wave propagation equation (see figure 1(a)) along with the Heisenberg–Langevin equations of motion yield

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{z}{{ \mathcal E }}_{r}(z,t) & = & -{\rm{i}}\sqrt{\displaystyle \frac{d{\gamma }_{{eg}}}{{cL}}}P(z,t),\\ {\partial }_{t}P(z,t) & = & -({\gamma }_{{eg}}+{\rm{i}}{\rm{\Delta }})P(z,t)+{\rm{i}}\sqrt{\displaystyle \frac{d{\gamma }_{{eg}}c}{L}}{{ \mathcal E }}_{r}(z,t)\\ & & +{\rm{i}}{{\rm{\Omega }}}_{R}(t)S(z,t)+{F}_{P}(z,t),\\ {\partial }_{t}S(z,t) & = & -{\gamma }_{0}S(z,t)+{\rm{i}}{{\rm{\Omega }}}_{R}^{* }(t)P(z,t)+{F}_{S}(z,t),\end{array}\end{eqnarray} \tag{ 1 }$$

where $P(z,t)=\sqrt{N}{\sigma }_{{ge}}(z,t){{\rm{e}}}^{-{\rm{i}}{\omega }_{1}\tfrac{L-z}{c}}$ and $S(z,t)=\sqrt{N}{\sigma }_{{gs}}(z,t){{\rm{e}}}^{-{\rm{i}}({\omega }_{1}-{\omega }_{2})\tfrac{L-z}{c}}$ are rescaled and slowly varying atomic operators (see appendix for details), with ω₁ (ω₂) referring to the energy transition of the $| e\rangle -| g\rangle$ ($| e\rangle -| s\rangle$) transition. γ_eg (γ₀) refers to the decay rate of the $| e\rangle -| g\rangle$ ($| g\rangle -| s\rangle$) transition. L denotes the length of the atomic sample and N the number of atoms within this sample. F_S and F_P indicate the noise operators associated to S and P, respectively. ${{\rm{\Omega }}}_{{\rm{R}}}$ (Δ) refers to the Rabi frequency (detuning) of the classical write control field on the $| e\rangle -| g\rangle$ transition, and ${{ \mathcal E }}_{{\rm{r}}}$ denotes the quantum field of the retrieval emission. The optical depth d characterises the absorption of resonant light in the sample, such that the outgoing light intensity is ${I}_{0}(z=L)={e}^{-2d}I(z=0)$, valid when the spectrum of the incoming light is well contained within the atomic bandwidth.

Here, we consider the situation where retrieval is completed well within the spin wave decoherence time, and thus ignore γ₀. In computing the spin and photon numbers, we also ignore the Langevin noise terms F_S and F_P as they appear in normal ordered form, and in the situation where almost all the atoms are in the ground state these do not contribute.

Defining first the reversed functions $\bar{P}(L-z,t)=P(z,t),\bar{S}(L-z,t)=S(z,t)$ and ${\bar{{ \mathcal E }}}_{{\rm{r}}}(L-z,t)\,={{ \mathcal E }}_{{\rm{r}}}(z,t)$, then taking the Laplace transforms of equations (1) from $L-z=z^{\prime} \to u$, we begin with the following set of transformed equations

$$\begin{eqnarray}&&{\bar{{ \mathcal E }}}_{{\rm{r}}}(u,t)={\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{{eg}}d}{{cL}}}\displaystyle \frac{1}{u}\bar{P}(u,t),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\partial }_{t}\bar{P}(u,t)=-\left[{\gamma }_{{eg}}\left(1+\displaystyle \frac{d}{{Lu}}\right)+{\rm{i}}{\rm{\Delta }}\right]\bar{P}(u,t)+{\rm{i}}{{\rm{\Omega }}}_{{\rm{R}}}(t)\bar{S}(u,t)\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\partial }_{t}\bar{S}(u,t)={\rm{i}}{{\rm{\Omega }}}_{{\rm{R}}}^{* }(t)\bar{P}(u,t).\end{eqnarray} \tag{ 4 }$$

From equations (3) and (4) we first obtain the following result

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}(\langle {\bar{P}}^{\dagger }({u}_{1},t)\bar{P}({u}_{2},t)+{\bar{S}}^{\dagger }({u}_{1},t)\bar{S}({u}_{2},t)\rangle )\\ \quad =\,{\gamma }_{{eg}}\left(-2-\displaystyle \frac{{d}}{{{Lu}}_{1}}-\displaystyle \frac{{d}}{{{Lu}}_{2}}\right)\langle {\bar{P}}^{\dagger }({u}_{1},t)\bar{P}({u}_{2},t)\rangle .\end{array}\end{eqnarray} \tag{ 5 }$$

With equation (2) we can then rewrite the number of emitted photons η in terms of P(u, t)

$$\begin{eqnarray*}\begin{array}{rcl}\eta & = & \displaystyle \frac{c}{L}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{t}\ \langle {{ \mathcal E }}_{r}^{\dagger }(z=0,t){{ \mathcal E }}_{r}(z=0,t)\rangle \\ & = & {\left.\displaystyle \frac{c}{L}{{ \mathcal L }}_{2}^{-1}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{t}\displaystyle \frac{{\gamma }_{{eg}}d}{{cL}}\displaystyle \frac{1}{{u}_{1}{u}_{2}}\langle {\bar{P}}^{\dagger }({u}_{1},t)\bar{P}({u}_{2},t)\rangle \right|}_{\displaystyle \genfrac{}{}{0em}{}{{z}_{1}^{{\prime} }\to L}{{z}_{2}^{{\prime} }\to L}},\end{array}\end{eqnarray*}$$

where ${{ \mathcal L }}_{2}^{-1}$ indicates the instruction to take the Laplace inverses of both u₁ and u₂ separately. With the use of equation (5) we can next rewrite $\langle {\bar{P}}^{\dagger }({u}_{1},t)\bar{P}({u}_{2},t)\rangle$ as a full derivative and perform the integral to get

$$\begin{eqnarray*}\begin{array}{rcl}\eta & = & {{ \mathcal L }}_{2}^{-1}\displaystyle \frac{d}{L}\displaystyle \frac{-1}{({u}_{1}+{u}_{2})d+2{{Lu}}_{1}{u}_{2}}\\ & & {{\times (\langle {\bar{P}}^{\dagger }({u}_{1},t)\bar{P}({u}_{2},t)\rangle +\langle {\bar{S}}^{\dagger }({u}_{1},t)\bar{S}({u}_{2},t)\rangle )| }_{t=0}^{\infty }| }_{\displaystyle \genfrac{}{}{0em}{}{{z}_{1}^{{\prime} }\to L}{{z}_{2}^{{\prime} }\to L}}\\ & = & {\left.\displaystyle \frac{1}{{L}^{2}}{{ \mathcal L }}_{2}^{-1}\displaystyle \frac{{dL}}{d({u}_{1}+{u}_{2})+2{{Lu}}_{1}{u}_{2}}\langle {\bar{S}}^{\dagger }({u}_{1},0)\bar{S}({u}_{2},0)\rangle \right|}_{\displaystyle \genfrac{}{}{0em}{}{{z}_{1}^{{\prime} }\to L}{{z}_{2}^{{\prime} }\to L}},\end{array}\end{eqnarray*}$$

where the last equality comes from the conditions we assume in a complete retrieval process, i.e. that we begin with only $| g\rangle -| s\rangle$ coherences and at the end of the process all atoms are in $| g\rangle$. By performing the inverse Laplace transforms one sees that for complete retrieval in the backward direction⁶ ,

$$\begin{eqnarray}\begin{array}{rcl}\eta & = & \displaystyle \frac{1}{L}{\displaystyle \int }_{0}^{L}{{\rm{d}}{z}}_{1}\displaystyle \frac{1}{L}{\displaystyle \int }_{0}^{L}{{\rm{d}}{z}}_{2}\ {k}_{r}(L-{z}_{1},L-{z}_{2})\\ & & \times \,\langle {S}^{\dagger }(L-{z}_{1},0)S(L-{z}_{2},0)\rangle ,\end{array}\end{eqnarray} \tag{ 6 }$$

where ${k}_{r}({z}_{1},{z}_{2})=\tfrac{d}{2}{{\rm{e}}}^{-d\tfrac{{z}_{1}+{z}_{2}}{2L}}{I}_{0}\left(\tfrac{d}{L}\sqrt{{z}_{1}{z}_{2}}\right)$ and I_n(x) indicates the modified nth Bessel function of the first kind. We proceed by considering the situation where there is originally a single spin wave in the sample (such that $\tfrac{1}{L}{\int }_{0}^{L}\,{S}^{\dagger }(z,0)S(z,0)\ {\rm{d}}{z}=1$), and thus interpret η as the efficiency of the retrieval process. The retrieval efficiency η is independent of the details of the read control field used, and is a result of the ratio between desired and undesired modes that are retrieved from the spin wave.

Having shown the dependence of the retrieval efficiency on the spin wave shape and optical depth, we now look to gain some intuition on how one might obtain the optimal retrieval efficiencies, by plotting the spin shapes that maximise the retrieval efficiency for given depths.

To do this, we recognise equation (6) as the continuous form of a product of discretised versions of k_r (in the form of a matrix) and $| S\rangle$ (in the form of a vector). Cast in this light, this integral can be computed by performing a matrix multiplication between the discretised versions of k_r and $| S\rangle$, and in this discrete approximation the optimal spin shape is the eigenvector of k_r with the largest eigenvalue. One can then interpolate the resulting vector to obtain optimised spin shapes, which are shown as solid coloured lines in figure 2.

Zoom In Zoom Out Reset image size

The best spin shapes for optimal retrieval show a clear spatial dependence with a bias (depending on the optical depth d) towards placing larger excitation probabilities towards the retrieval direction (backwards in this case). Intuitively, we see these shapes representing the best way to obtain constructive interference throughout the retrieval process. As the optical field is converted from the spin wave towards the retrieval direction, it benefits from encountering a higher excitation from the atoms it next encounters. We will denote the retrieval efficiencies from these optimal spin shapes as ${\eta }^{* }$.

In the previous section, we have outlined how the retrieval efficiency depends on the shape of the given spin excitation, and also how the optimal spin shapes can be computed. Here we propose a method of conveniently creating spin shapes that yield near-optimal retrieval efficiencies. In contrast to creating spin excitations using spontaneous Raman processes enabled by far-detuned write control fields, we explore the use of resonant control fields instead, which create spin excitations with significantly position-dependent excitation profiles. These profiles can be controlled by tuning the duration of the write control field, which is in turn related to the frequency spread. A shorter (longer) write field duration implies that it has a wider (sharper) spread in frequency, and this thus affects how quickly the write pulse is depleted within the sample.

We give a detailed derivation of the write process in appendix B.1. To summarise (see figure 1(b)), beginning with all atoms in the $| g\rangle$-level, we send a short rising exponential resonant write pulse with Rabi frequency ${{\rm{\Omega }}}_{W}(0,t)={{\rm{\Omega }}}_{W}^{\max }{{\rm{e}}}^{t/{\tau }_{W}}$ that does not significantly excite the atoms to the $| e\rangle$ level (${{\rm{\Omega }}}_{W}^{\max }{\tau }_{W}\ll 1$). If sent with a sufficiently short duration (τ_W ≪ 1/γ_eg, τ_W ≪ 1/γ_es) and shut off at t = 0, one can consider only the dynamics along the $| g\rangle -| e\rangle$ transition, and obtain atomic coherences of the form (see appendix B.2)

$$\begin{eqnarray}&&{\sigma }_{{ge}}(z,0)={{\rm{e}}}^{{{\rm{i}}{k}}_{{\rm{w}}}.z}{\theta }_{0}{{\rm{e}}}^{-\tfrac{\alpha z}{2}},\end{eqnarray} \tag{ 7 }$$

where ${\theta }_{0}={\rm{i}}\tfrac{{{\rm{\Omega }}}_{{\rm{W}}}^{\max }{\tau }_{{\rm{W}}}}{1+{\gamma }_{{eg}}{\tau }_{{\rm{W}}}},\alpha /2=d\tfrac{{\gamma }_{{eg}}{\tau }_{{\rm{W}}}}{1+{\gamma }_{{eg}}{\tau }_{{\rm{W}}}}\tfrac{1}{L}$ and k_w indicates the wave vector for the write photon, which is described using a quantum field ${{ \mathcal E }}_{{\rm{w}}}$. Immediately after the preparation, we look for the detection of write photons within a short detection window τ_d as a herald for single spin excitations. This avoids potential dephasing effects from the decoherence of the $| e\rangle$ level. In this short detection window of duration ${\tau }_{{\rm{d}}}\ll \min \left(\tfrac{1}{2{\gamma }_{{es}}},\tfrac{1}{2{\gamma }_{{eg}}}\right)$, and where ${\tau }_{{\rm{d}}}\ll {\left\{\bar{d}{\gamma }_{{es}}| {\theta }_{0}{| }^{2}\tfrac{1-{{\rm{e}}}^{-\alpha L}}{\alpha L}\right\}}^{-1}$, ensuring the number of emitted write photons n_w is much smaller than 1, we obtain (see appendices B.3 and B.4)

$$\begin{eqnarray}&&{n}_{{\rm{w}}}=(\bar{d}\ {\gamma }_{{es}}{\tau }_{{\rm{d}}})| {\theta }_{0}{| }^{2}\displaystyle \frac{1-{{\rm{e}}}^{-\alpha L}}{\alpha L},\end{eqnarray} \tag{ 8 }$$

where γ_es ($\bar{d}$) refers to the decay rate (optical depth) of the $| e\rangle -| s\rangle$ transition. The write photon number η_w is simply the product of $\bar{d}\ {\gamma }_{{es}}{\tau }_{{\rm{d}}}$ and the fraction of excited atoms (averaged across the sample).

In this same regime for τ_d, to leading order the corresponding heralded spin state is (see appendix B.5)

$$\begin{eqnarray}&&{S}^{\dagger }(z,{\tau }_{{\rm{d}}})=-{\rm{i}}\sqrt{\displaystyle \frac{\bar{d}{\gamma }_{{es}}c}{L}}{\theta }_{0}{{\rm{e}}}^{-\alpha z/2}{\int }_{0}^{{\tau }_{{\rm{d}}}}{{\rm{e}}}^{-{\gamma }_{0}({\tau }_{{\rm{d}}}-{t}_{{\rm{a}}})}{{ \mathcal E }}_{{\rm{w}}}(0,{t}_{{\rm{a}}})\,{{\rm{d}}{t}}_{{\rm{a}}},\end{eqnarray} \tag{ 9 }$$

which has an exponentially decaying spatial dependence from the z = 0 side of the sample. The extent of this spatial decay is characterised by α, which does depend on the given properties of the atomic sample, but can be controlled by varying the write control field duration τ_W. One can compare this class of heralded spin shapes (created from exponentially rising write pulses) to the optimal spin shapes in figure 2.

We now proceed with the retrieval process, and spell out the exact requirements for a certain implementation of retrieval—the fast π-pulse using a square waveform of duration τ_R. Once again, we focus on retrieval processes completed well within the spin wave decoherence time and performed under relevant experimental conditions. We thus ignore both the spin decoherence and Langevin noise terms in equation (1). Here we have implicitly assumed that the energy levels of the $| g\rangle$ and $| s\rangle$ levels are degenerate⁷ . See [15, 17] for details.

With a resonant square retrieve pulse in the backward direction (see figure 1(c)) one finds the following simple expression for the dynamics of the spin wave (details given in appendix A.1)

$$\begin{eqnarray}&&\ddot{\bar{S}}(u,t)+A\dot{\bar{S}}(u,t)+B\bar{S}(u,t)=0,\end{eqnarray} \tag{ 10 }$$

where $A={\gamma }_{{eg}}\left(1+\tfrac{d}{{Lu}}\right)$ and $B={{\rm{\Omega }}}_{{\rm{R}}}^{2}$ (for real Ω_R), and we have taken the Laplace transform $L-z=z^{\prime} \to u$.

In the regime⁸ where $2{{\rm{\Omega }}}_{{\rm{R}}}\gg {\gamma }_{{eg}}(1+d)$, we find $4B\gg {A}^{2}$, and obtain the following solution

$$\begin{eqnarray}&&\bar{S}(u,t)={{\rm{e}}}^{-{At}/2}\cos ({{\rm{\Omega }}}_{{\rm{R}}}t)\bar{S}(u,t={\tau }_{{\rm{d}}}),\end{eqnarray} \tag{ 11 }$$

which yields the following expression

$$\begin{eqnarray}\begin{array}{rcl}\bar{P}(u,t) & = & \displaystyle \frac{1}{{\rm{i}}{{\rm{\Omega }}}_{{\rm{R}}}}{\partial }_{{\rm{t}}}\bar{S}(u,t)\\ & = & \displaystyle \frac{{\rm{i}}}{{{\rm{\Omega }}}_{{\rm{R}}}}{{\rm{e}}}^{-\tfrac{A}{2}t}\left(\displaystyle \frac{A}{2}\cos ({{\rm{\Omega }}}_{{\rm{R}}}t)+{{\rm{\Omega }}}_{{\rm{R}}}\sin ({{\rm{\Omega }}}_{{\rm{R}}}t)\right)\bar{S}(u,t={\tau }_{{\rm{d}}}),\end{array}\end{eqnarray} \tag{ 12 }$$

where we then see that with a sufficiently fast π-pulse (such that 2Ω_Rτ_R = π) obeying γ_eg (1 + d) τ_R ≪ 2, one can convert S to P without loss, yielding

$$\begin{eqnarray}&&\bar{P}(u,{\tau }_{{\rm{R}}}+{\tau }_{{\rm{d}}})\approx {\rm{i}}\bar{S}(u,t={\tau }_{{\rm{d}}}).\end{eqnarray} \tag{ 13 }$$

The emitted read photon field can then be obtained by solving the set of equations in (1) after the fast read control field has ended (see appendix A.2), giving

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal E }}_{{\rm{r}}}(0,t) & = & {\rm{i}}\sqrt{\displaystyle \frac{{\gamma }_{{eg}}d}{{cL}}}{{\rm{e}}}^{-{\gamma }_{{eg}}t}{\displaystyle \int }_{0}^{L}{J}_{0}\left[2\sqrt{\displaystyle \frac{{\gamma }_{{eg}}d}{L}t(L-{z}_{1}^{{\prime\prime} })}\right]\\ & & P(L-{z}_{1}^{{\prime\prime} },{\tau }_{{\rm{R}}}+{\tau }_{{\rm{d}}})\,{{\rm{d}}{z}}_{1}^{{\prime\prime} }.\end{array}\end{eqnarray} \tag{ 14 }$$

Along with equation (13) and noting that ${\int }_{0}^{\infty }{{\rm{e}}}^{-\alpha x}{J}_{\nu }(2\beta \sqrt{x}){J}_{\nu }(2\gamma \sqrt{x}){\rm{d}}{x}=\tfrac{1}{\alpha }{I}_{\nu }\left(\tfrac{2\beta \gamma }{\alpha }\right)\exp \left(-\tfrac{{\beta }^{2}+{\gamma }^{2}}{\alpha }\right)$ [18], this emitted field then yields a retrieval efficiency given by equation (6).

We have seen that the proposed retrieval protocol yields a dependence on the spin shape, as described by equation (6). Hence we now compare the retrieval efficiencies attainable with our protocol and compare them to the optimal ones.

We can estimate the achievable efficiency of our protocol by choosing a write pulse duration such that the heralded spin shape best fits the optimal spin shape. A good approximation to this write pulse duration is well described in [19], and given by

$$\begin{eqnarray}&&{\tau }_{{\rm{W}}}^{\mathrm{approx}}=\displaystyle \frac{1}{{\gamma }_{{eg}}}\displaystyle \frac{1}{1+\tfrac{d}{2}}.\end{eqnarray} \tag{ 15 }$$

As we show below, this simple expression for the write pulse duration essentially produces the optimal efficiency available for a given optical depth. This is hence the write pulse duration we recommend.

We also compute the retrieval efficiencies η^fwd that would be obtained if the resonant write pulse of duration ${\tau }_{W}^{\mathrm{approx}}$ were to be followed by a co-propagating retrieve pulse instead. This would result in a situation where the spin wave would be far from optimal with respect to the retrieval direction. In figure 3, we compare the optimal efficiency ${\eta }^{* }$, the efficiencies η^res and η^fwd obtained with our proposal (from a spin wave created from a resonant exponential pulse with duration ${\tau }_{W}^{\mathrm{approx}}$) together with the efficiency of the standard approach using far off-resonant write pulses, for which the efficiency is bounded by the complete retrieval efficiency from a flat spin wave [15]

$$\begin{eqnarray}&&{\eta }^{\mathrm{off}-\mathrm{res}}=1-{{\rm{e}}}^{-d}({I}_{0}(d)+{I}_{1}(d)),\end{eqnarray} \tag{ 16 }$$

which we have verified numerically. This retrieval efficiency is valid for retrieval from both the forward and backward directions from a flat spin wave.

Zoom In Zoom Out Reset image size

Our proposal approaches optimal efficiencies, performing within ∼10⁻³ of ${\eta }^{* }$ and compares favourably with respect to the standard off-resonant case. The improvement in efficiency is dependent on the optical depth, and we present some values in table 1.

Table 1.  Comparison of retrieval efficiency from different spin shapes.

d	${\eta }^{\mathrm{fwd}}$	${\eta }^{\mathrm{off}-\mathrm{res}}$	${\eta }^{\mathrm{res}}$	${\eta }^{* }$
0.1	0.047 6	0.047 6	0.047 6	0.047 6
1	0.314 0	0.326 3	0.330 5	0.330 5
10	0.567 1	0.750 9	0.813 4	0.814 2
20	0.618 3	0.822 7	0.892 1	0.897 3
100	0.760 0	0.920 3	0.972 8	0.974 5

Finally, we have also investigated the on-resonance retrieval with an exponentially increasing write pulse using non optimised pulse widths $({\tau }_{W}={\gamma }_{{eg}}^{-1})$, but find a saturation of only ∼67 % of the retrieval efficiency for high optical depths.

For a single photon source to be useful, one needs to not only efficiently emit a single photon, but also to ensure that the said photon is emitted in a single mode. In our model we have assumed that the write and read photons are each in a single mode.

For an actual implementation, one way to check that a single mode for the write and read photons are collected and detected, is to perform an autocorrelation measurement with two detectors after a 50-50 beamsplitter (see appendix D). Under the assumption that the write and read photons are correlated through vacuum squeezing processes, this measurement allows us to determine the number of emission modes K, as it gives ${g}^{(2)}\sim 1+\tfrac{1}{K}$ [20] (valid in the absence of detector noise and for small emission probabilities).

We begin from the Hamiltonian $H={H}_{0}+V$ (see [15]), where we consider an atomic sample of length L, and a classical field sent from the z = L side of the sample. Choosing $| g\rangle$ to be the energy level reference for the atomic states, we have

$$\begin{eqnarray}&&{H}_{0}=\int {\rm{d}}\omega {\hslash }\omega {\hat{a}}_{\omega }^{\dagger }{\hat{a}}_{\omega }+\displaystyle \sum _{i=1}^{N}({\hslash }{\omega }_{s}{\sigma }_{{ss}}^{j}+{\hslash }{\omega }_{e}{\sigma }_{{ee}}^{j})\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}\begin{array}{rcl}V & = & -{\hslash }\displaystyle \sum _{i=1}^{N}\left({{\rm{\Omega }}}_{{\rm{R}}}\left(t-\displaystyle \frac{L-{z}_{i}}{c}\right){\sigma }_{{es}}^{i}{{\rm{e}}}^{-{\rm{i}}{\omega }_{2}t}{{\rm{e}}}^{+{\rm{i}}{\omega }_{2}\left(\tfrac{L-{z}_{i}}{c}\right)}\right.\\ & & \left.+g\sqrt{\displaystyle \frac{L}{2\pi c}}\displaystyle \int {\rm{d}}\omega \ {a}_{\omega }{{\rm{e}}}^{{\rm{i}}\omega \tfrac{L-{z}_{i}}{c}}{\sigma }_{{eg}}^{i}+{\rm{h}}{\rm{.}}{\rm{c}}.\right),\end{array}\end{eqnarray} \tag{ 18 }$$

where ${\sigma }_{\mu \nu }^{i}=| \mu {\rangle }_{i}\langle \nu |$ indicates atomic level operators for the ith atom, and a_w indicates the annihilation operator for the photonic mode at frequency ω. ω₂ (ω₁) indicates the frequency of the read control (photon) field, respectively. Note that we are considering resonant pulses, so we have ω₁ (ω₂) = ω_e (ω_s). Using

$$\begin{eqnarray*}\begin{array}{rcl}A & = & \displaystyle \sum _{i=1}^{N}[{\hslash }({\omega }_{1}-{\omega }_{2}){\sigma }_{{ss}}^{i}+{\hslash }{\omega }_{1}{\sigma }_{{ee}}^{i}]+{\hslash }{\omega }_{1}\displaystyle \int {\rm{d}}\omega \ {{ \mathcal E }}_{r}^{\dagger }(z,t){{ \mathcal E }}_{r}(z,t),\\ U & = & {{\rm{e}}}^{-{\rm{i}}{At}/{\hslash }},\end{array}\end{eqnarray*}$$

for the change of frame, then in the continuum limit, we obtain

$$\begin{eqnarray*}\begin{array}{rcl}{H}_{{\rm{new}}} & = & {U}^{\dagger }{HU}-A\\ & = & \displaystyle \int {\rm{d}}\omega \ {\hslash }\omega {a}_{\omega }^{\dagger }{a}_{\omega }-{\hslash }{\omega }_{1}\displaystyle \int {\rm{d}}z\ {{ \mathcal E }}_{r}^{\dagger }(z,t){{ \mathcal E }}_{r}(z,t)\\ & & +\displaystyle \frac{N}{L}\displaystyle \int {\rm{d}}z\{-{\hslash }{{\rm{\Omega }}}_{{\rm{R}}}(z,t){\sigma }_{{es}}(z,t){{\rm{e}}}^{+{\rm{i}}{\omega }_{2}\tfrac{L-z}{c}}+{\rm{H}}{\rm{.}}{\rm{c}}.\\ & & -g{{ \mathcal E }}_{r}(z,t){\sigma }_{{eg}}(z,t){{\rm{e}}}^{+{\rm{i}}{\omega }_{1}\tfrac{L-z}{c}}+{\rm{H}}{\rm{.}}{\rm{c}}.\},\end{array}\end{eqnarray*}$$

where we have defined a real Rabi frequency ${{\rm{\Omega }}}_{{\rm{R}}}(z,t)={{\rm{\Omega }}}_{{\rm{R}}}(t-\tfrac{L-z}{c})$ and ${{ \mathcal E }}_{{\rm{r}}}(z,t)\,=\sqrt{\tfrac{L}{2\pi c}}\int {\rm{d}}\omega \ {{\rm{e}}}^{{\rm{i}}{\omega }_{1}t}{a}_{\omega }{{\rm{e}}}^{{\rm{i}}(\omega -{\omega }_{1})\tfrac{L-z}{c}}$. Using the field propagation equation along with the Heisenberg–Langevin equations of motion, we have in a moving coordinate frame, ignoring spinwave decoherence and the noise terms, and also considering that σ_gg ∼ 1,

$$\begin{eqnarray*}\begin{array}{rcl}{\partial }_{z}{{ \mathcal E }}_{{\rm{r}}}(z,t) & = & -\displaystyle \frac{{\rm{i}}{g}\sqrt{N}}{c}P(z,t)\\ {\partial }_{{\rm{t}}}P(z,t) & = & -{\gamma }_{{eg}}P(z,t)+{\rm{i}}{g}\sqrt{N}{{ \mathcal E }}_{{\rm{r}}}(z,t)+{\rm{i}}{{\rm{\Omega }}}_{{\rm{R}}}(t)S(z,t)\\ {\partial }_{{\rm{t}}}S(z,t) & = & {\rm{i}}{{\rm{\Omega }}}_{{\rm{R}}}(t)P(z,t),\end{array}\end{eqnarray*}$$

where ${g}^{2}N=\tfrac{d{\gamma }_{{eg}}c}{L}$, $P(z,t)=\sqrt{N}{\sigma }_{{ge}}(z,t){{\rm{e}}}^{-{\rm{i}}{\omega }_{1}\tfrac{L-z}{c}}$ and $S(z,t)=\sqrt{N}{\sigma }_{{gs}}(z,t){{\rm{e}}}^{-{\rm{i}}({\omega }_{1}-{\omega }_{2})\tfrac{L-z}{c}}$. In the continuum limit, the spin and field operators obey the following commutation relations

$$\begin{eqnarray}&&[{\sigma }_{\alpha \beta }(z,t),{\sigma }_{\mu \nu }(z^{\prime} ,t)]=\displaystyle \frac{L}{N}\delta (z-z^{\prime} )({\delta }_{\beta \mu }{\sigma }_{\alpha \nu }(z,t)-{\delta }_{\nu \alpha }\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&[{{ \mathcal E }}_{{\rm{r}}}(z,t),{{ \mathcal E }}_{{\rm{r}}}^{\dagger }(z,t^{\prime} )]=\displaystyle \frac{L}{c}\delta (t-t^{\prime} ).\end{eqnarray} \tag{ 20 }$$

Rewriting equations (A.1) in the reverse direction e.g. ${\bar{{ \mathcal E }}}_{{\rm{r}}}(z^{\prime} ,t)={\bar{{ \mathcal E }}}_{{\rm{r}}}(L-z,t)={{ \mathcal E }}_{{\rm{r}}}(z,t)$, and taking the Laplace transform from $z^{\prime} \to u$ we obtain

$$\begin{eqnarray}\begin{array}{rcl} & & {\bar{{ \mathcal E }}}_{{\rm{r}}}(u,{t})=\displaystyle \frac{{\rm{i}}{g}\sqrt{N}}{{cu}}\bar{P}(u,{t})+\displaystyle \frac{1}{u}{\bar{{ \mathcal E }}}_{{\rm{r}}}(z^{\prime} =0,{t}),\\ & & {\partial }_{{\rm{t}}}\bar{P}(u,{t})=-{\gamma }_{{eg}}\bar{P}(u,{t})+{\rm{i}}{g}\sqrt{N}{\bar{{ \mathcal E }}}_{{\rm{r}}}(u,{\rm{t}})+{\rm{i}}{{\rm{\Omega }}}_{{\rm{R}}}({t})\bar{S}(u,{t}),\\ & & {\partial }_{{\rm{t}}}\bar{S}(u,{t})={\rm{i}}{{\rm{\Omega }}}_{{\rm{R}}}({\rm{t}})\bar{P}(u,{t}).\end{array}\end{eqnarray} \tag{ 21 }$$

We can combine these three equations into a single differential equation, where we have ignored the boundary term ${\bar{{ \mathcal E }}}_{{\rm{r}}}(z^{\prime} =0,t)$ since we send the read control field into the z = L side of the atoms. On resonance (Δ = 0), let $A={\gamma }_{{eg}}+\tfrac{{g}^{2}N}{{cu}}$ and $B={{\rm{\Omega }}}_{{\rm{R}}}^{2}$ to see

$$\begin{eqnarray}&&\ddot{\bar{S}}(u,t)+A\dot{\bar{S}}(u,t)+B\bar{S}(u,t)=0.\end{eqnarray} \tag{ 22 }$$

In the strong regime for the read control field, one requires $2| {{\rm{\Omega }}}_{{\rm{R}}}| \gg {\gamma }_{{eg}}(1+d)$, which implies

$$\begin{eqnarray*}\begin{array}{rcl}2{{\rm{\Omega }}}_{{\rm{R}}} & \gg & {\gamma }_{{eg}}\left(1+\displaystyle \frac{{dz}^{\prime} }{L}\right)\\ \Rightarrow \,2{{\rm{\Omega }}}_{{\rm{R}}} & \gg & {\gamma }_{{eg}}\left(1+\displaystyle \frac{d}{{Lu}}\right),\end{array}\end{eqnarray*}$$

which then yields the regime $4B\gg {A}^{2}$.

The solution for equation (22) in this regime is

$$\begin{eqnarray*}&&\bar{S}(u,t)={{\rm{e}}}^{-{At}/2}\cos ({{\rm{\Omega }}}_{{\rm{R}}}t){C}_{1}(u)+{{\rm{e}}}^{-{At}/2}\sin ({{\rm{\Omega }}}_{{\rm{R}}}t){C}_{2}(u),\end{eqnarray*}$$

where the initial condition implies

$$\begin{eqnarray*}&&\bar{S}(u,t)={{\rm{e}}}^{-{At}/2}\cos ({{\rm{\Omega }}}_{{\rm{R}}}t)\bar{S}(u,t=0).\end{eqnarray*}$$

One can then find the prepared polarisation in terms of the intial spin condition,

$$\begin{eqnarray*}\begin{array}{rcl}\bar{P}(u,t) & = & \displaystyle \frac{1}{{\rm{i}}{{\rm{\Omega }}}_{R}}{\partial }_{t}\bar{S}(u,t)\\ & = & \displaystyle \frac{{\rm{i}}}{{{\rm{\Omega }}}_{R}}{{\rm{e}}}^{-\tfrac{A}{2}t}\left[\displaystyle \frac{A}{2}\cos ({{\rm{\Omega }}}_{R}t)+{{\rm{\Omega }}}_{R}\sin ({{\rm{\Omega }}}_{R}t)\right]\bar{S}(u,t=0).\end{array}\end{eqnarray*}$$

In the limit where we have a sufficiently strong read control field $(2{{\rm{\Omega }}}_{{\rm{R}}}\gg \tfrac{\pi }{2}{\gamma }_{{eg}}(1+d))$, the π-pulse is completed quickly and we obtain a lossless preparation of $\bar{P}(u,t)$ from $\bar{S}(u,t=0)$ in the form

$$\begin{eqnarray}&&\bar{P}(u,{\tau }_{{\rm{R}}})={\rm{i}}\bar{S}(u,t=0).\end{eqnarray} \tag{ 23 }$$

Once the polarisation is prepared, we find the emission by solving for the dynamics in the absence of the laser,

$$\begin{eqnarray*}\begin{array}{rcl} & & {\partial }_{z}{\bar{{ \mathcal E }}}_{{\rm{r}}}(z,t)=-{\rm{i}}\displaystyle \frac{g\sqrt{N}}{c}\bar{P}(z,t),\\ & & ({\partial }_{t}+{\gamma }_{{eg}})\bar{P}(z,t)={\rm{i}}{g}\sqrt{N}{\bar{{ \mathcal E }}}_{{\rm{r}}}(z,t).\end{array}\end{eqnarray*}$$

Taking the Laplace transform from $L-z=z^{\prime} \to u$ and neglecting the boundary term since it does not contribute to the photon number, we have

$$\begin{eqnarray*}\begin{array}{rcl} & & {\bar{{ \mathcal E }}}_{{\rm{r}}}(u,t)={\rm{i}}\displaystyle \frac{g\sqrt{N}}{{cu}}\bar{P}(u,t),\\ & & ({\partial }_{t}+{\gamma }_{{eg}})\bar{P}(u,t)={\rm{i}}{g}\sqrt{N}{\bar{{ \mathcal E }}}_{{\rm{r}}}(u,t)=-\displaystyle \frac{{g}^{2}N}{{cu}}\bar{P}(u,t).\end{array}\end{eqnarray*}$$

This yields the evolution of P(u, t) after its preparation from S(u, t),

$$\begin{eqnarray}&&\bar{P}(u,t)={{\rm{e}}}^{-\left({\gamma }_{{eg}}+\tfrac{{g}^{2}N}{{cu}}\right)(t-{\tau }_{{\rm{R}}})}\bar{P}(u,{\tau }_{{\rm{R}}}),\end{eqnarray} \tag{ 24 }$$

and gives an emitted field of

$$\begin{eqnarray*}\begin{array}{rcl}{\bar{{ \mathcal E }}}_{{\rm{r}}}(z^{\prime} ,t) & = & {\rm{i}}\displaystyle \frac{g\sqrt{N}}{c}{{\rm{e}}}^{-{\gamma }_{{eg}}(t-{\tau }_{{\rm{R}}})}\\ & & \times \,{\displaystyle \int }_{0}^{z^{\prime} }{\rm{d}}{z}^{\prime\prime} {J}_{0}\left[2\sqrt{\displaystyle \frac{{g}^{2}N}{c}(t-{\tau }_{{\rm{R}}})(z^{\prime} -z^{\prime\prime} )}\right]\bar{P}(z^{\prime} ,{\tau }_{{\rm{R}}}),\end{array}\end{eqnarray*}$$

where ${J}_{n}[x]$ refers the nth Bessel function of the first kind. Now with $z^{\prime} =L-z$, we require the field at z = 0 for backward retrieval, and we finally obtain

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal E }}_{r}(0,t) & = & -\displaystyle \frac{g\sqrt{N}}{c}{{\rm{e}}}^{-{\gamma }_{{eg}}t}\\ & & \times \,{\displaystyle \int }_{0}^{L}{\rm{d}}{z}^{\prime\prime} {J}_{0}\left[2\sqrt{\displaystyle \frac{{g}^{2}N}{c}t(L-z^{\prime\prime} )}\right]S(L-z^{\prime\prime} ,0),\end{array}\end{eqnarray} \tag{ 25 }$$

where we have used equation (23) for a lossless preparation.

In the weak regime for the read control field, one requires $2| {{\rm{\Omega }}}_{{\rm{R}}}| \ll {\gamma }_{{eg}}$, which implies

$$\begin{eqnarray*}\begin{array}{rcl}2{{\rm{\Omega }}}_{R} & \ll & {\gamma }_{{eg}}\left(1+\displaystyle \frac{{dz}}{L}\right)\\ \Rightarrow 2{{\rm{\Omega }}}_{R} & \ll & {\gamma }_{{eg}}\left(1+\displaystyle \frac{d}{{Lu}}\right),\end{array}\end{eqnarray*}$$

which then yields the regime $4B\ll {A}^{2}$.

The solution for equation (22) in this regime is

$$\begin{eqnarray*}&&\bar{S}(u,t)={{\rm{e}}}^{-\tfrac{1}{2}(A+\sqrt{{A}^{2}-4B})t}{C}_{u}(1)+{{\rm{e}}}^{-\tfrac{1}{2}(A-\sqrt{{A}^{2}-4B})t}{C}_{u}(2).\end{eqnarray*}$$

When there is no laser (B = 0), there should be no spinwave decay since we have considered zero spin wave decoherence, so we set C_u(1) = 0 and obtain

$$\begin{eqnarray*}&&\bar{S}(u,t)={{\rm{e}}}^{-\tfrac{1}{2}(A-\sqrt{{A}^{2}-4B})t}\bar{S}(u,t=0).\end{eqnarray*}$$

Now, in this regime when the Rabi frequency is small, we have

$$\begin{eqnarray*}\begin{array}{rcl}{{\rm{e}}}^{-\tfrac{1}{2}(A-\sqrt{{A}^{2}-4B})t} & = & {{\rm{e}}}^{-\tfrac{1}{2}(A-A\sqrt{1-\tfrac{4B}{{A}^{2}}})t}\\ & \approx & {{\rm{e}}}^{-\tfrac{B}{A}t}\\ & = & {{\rm{e}}}^{-\tfrac{{{\rm{\Omega }}}^{2}}{{\gamma }_{{eg}}(1+\tfrac{d}{{Lu}})}t}\end{array}\end{eqnarray*}$$

This gives

$$\begin{eqnarray*}&&\bar{S}(u,t)={{\rm{e}}}^{-{Kt}\tfrac{1}{1+s/u}}\bar{S}(u,t=0),\end{eqnarray*}$$

where $K=\tfrac{{{\rm{\Omega }}}_{{\rm{R}}}^{2}}{{\gamma }_{{eg}}}$ and $s=\tfrac{d}{L}$. One can proceed to find $\bar{P}(u,{\rm{t}})=\tfrac{1}{{\rm{i}}{{\rm{\Omega }}}_{{\rm{R}}}}{\partial }_{{\rm{t}}}\bar{S}(u,{\rm{t}})$ and $\bar{{ \mathcal E }}(u,t)={\rm{i}}\tfrac{g\sqrt{N}}{{cu}}\bar{P}(u,t)$, giving

$$\begin{eqnarray*}&&\bar{{ \mathcal E }}(u,t)=-\displaystyle \frac{g\sqrt{N}}{c}\displaystyle \frac{K}{{{\rm{\Omega }}}_{{\rm{R}}}}\left[\displaystyle \frac{1}{u+s}{{\rm{e}}}^{-{Kt}+{Kt}\left(\tfrac{s}{s+u}\right)}\right]\bar{S}(u,t=0).\end{eqnarray*}$$

This yields

$$\begin{eqnarray}\begin{array}{rcl}\bar{{ \mathcal E }}(z^{\prime} ,t) & = & -\displaystyle \frac{g\sqrt{N}}{c}\displaystyle \frac{K}{{\rm{\Omega }}}{{\rm{e}}}^{-{Kt}}\\ & & \times \,{\displaystyle \int }_{0}^{z^{\prime} }{{\rm{e}}}^{-s(z^{\prime} -z^{\prime\prime} )}{I}_{0}(2\sqrt{{Kts}(z^{\prime} -z^{\prime\prime} )})\bar{S}(z^{\prime\prime} ,t=0){\rm{d}}{z}^{\prime\prime} .\end{array}\end{eqnarray} \tag{ 26 }$$

One can then compute the retrieval efficiency from a single spin wave, and this is found to yield the optimal retrieval efficiency.

$$\begin{eqnarray*}&&\begin{array}{l}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{t}\ \displaystyle \frac{c}{L}\langle {{ \mathcal E }}^{\dagger }(0,t){ \mathcal E }(0,t)\rangle \\ \quad =\,{\displaystyle \int }_{0}^{\infty }{\rm{d}}{t}\ \displaystyle \frac{d}{{L}^{2}}\displaystyle \frac{{{\rm{\Omega }}}_{R}^{2}}{{\gamma }_{{eg}}}{{\rm{e}}}^{-2{Kt}}{\displaystyle \int }_{0}^{L}{{\rm{d}}{z}}_{1}^{{\prime\prime} }\ {\displaystyle \int }_{0}^{L}{{\rm{d}}{z}}_{2}^{{\prime\prime} }\ {{\rm{e}}}^{-\tfrac{d}{L}(2L-{z}_{1}^{{\prime\prime} }-{z}_{2}^{{\prime\prime} })}\\ \times \,{I}_{0}\left[2\sqrt{{Kt}\displaystyle \frac{d}{L}(L-{z}_{1}^{{\prime\prime} })}\right]{I}_{0}\left[2\sqrt{{Kt}\displaystyle \frac{d}{L}(L-{z}_{2}^{{\prime\prime} })}\right]\\ \times \,\langle {S}^{\dagger }(L-{z}_{1}^{{\prime\prime} },t=0)S(L-{z}_{2}^{{\prime\prime} },t=0)\rangle \\ \quad =\,\displaystyle \frac{1}{L}{\displaystyle \int }_{0}^{L}{{\rm{d}}{z}}_{1}^{{\prime\prime} }\ \displaystyle \frac{1}{L}{\displaystyle \int }_{0}^{L}{{\rm{d}}{z}}_{2}^{\prime\prime} \ \displaystyle \frac{d}{2}{{\rm{e}}}^{-\tfrac{d}{2}\tfrac{(L-{z}_{1}^{{\prime\prime} })+(L-{z}_{2}^{{\prime\prime} })}{L}}\\ \times \,{I}_{0}\left[{\rm{d}}\sqrt{\displaystyle \frac{L-{z}_{1}^{{\prime\prime} }}{L}}\sqrt{\displaystyle \frac{L-{z}_{2}^{{\prime\prime} }}{L}}\right]\langle {S}^{\dagger }(L-{z}_{1}^{{\prime\prime} },t=0)S(L-{z}_{2}^{{\prime\prime} },t=0)\rangle \end{array}\end{eqnarray*}$$

where I_n[x] denotes the nth modified Bessel function of the first kind. We have made use of the fact that ${I}_{n}(z)={{\rm{i}}}^{-n}{J}_{n}({\rm{i}}{z})$ and also ${\int }_{0}^{\infty }{{\rm{e}}}^{-\alpha x}{J}_{\nu }(2\beta \sqrt{x}){J}_{\nu }(2\gamma \sqrt{x}){\rm{d}}{x}=\tfrac{1}{\alpha }{I}_{\nu }\left(\tfrac{2\beta \gamma }{\alpha }\right)\exp \left(-\tfrac{{\beta }^{2}+{\gamma }^{2}}{\alpha }\right)$.

The goal here is to first derive the expressions for the evolution of the atomic coherences in the write process. We begin from the Hamiltonian $\bar{H}={\bar{H}}_{0}+\bar{V}$

$$\begin{eqnarray}&&{\bar{H}}_{0}=\int {\rm{d}}\omega {\hslash }\omega {a}_{\omega }^{\dagger }{a}_{\omega }+\displaystyle \sum _{i=1}^{N}({\hslash }{\omega }_{s}{\sigma }_{{ss}}^{j}+{\hslash }{\omega }_{e}{\sigma }_{{ee}}^{j})\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}\begin{array}{rcl}\bar{V} & = & -{\hslash }\displaystyle \sum _{i=1}^{N}({{\rm{\Omega }}}_{{\rm{W}}}(t-{z}_{i}/c){\sigma }_{{eg}}^{i}{{\rm{e}}}^{-{\rm{i}}{\omega }_{1}(t-{z}_{i}/c)}\\ & & \left.+\bar{g}\sqrt{\displaystyle \frac{L}{2\pi c}}\displaystyle \int {\rm{d}}\omega \ {\hat{a}}_{\omega }{{\rm{e}}}^{{\rm{i}}\omega {z}_{i}/c}{\sigma }_{{es}}^{i}+{\rm{H}}{\rm{.}}{\rm{c}}.\right).\end{array}\end{eqnarray} \tag{ 28 }$$

Using

$$\begin{eqnarray*}\begin{array}{rcl}\bar{A} & = & \displaystyle \sum _{i=1}^{N}({\hslash }{\omega }_{s}{\sigma }_{{ss}}^{i}+{\hslash }{\omega }_{e}{\sigma }_{{ee}}^{i})+{\hslash }{\omega }_{2}\displaystyle \int {\rm{d}}z\ {{ \mathcal E }}_{w}^{\dagger }(z,t){{ \mathcal E }}_{w}(z,t),\\ \bar{U} & = & {{\rm{e}}}^{-{\rm{i}}\bar{A}t/{\hslash }}\end{array}\end{eqnarray*}$$

for the change of frame, then in the continuum limit, we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\bar{H}}_{{\rm{new}}} & = & {\bar{U}}^{\dagger }\bar{H}\bar{U}-\bar{A}\\ & = & \displaystyle \int {\rm{d}}\omega \ {\hslash }{a}_{\omega }^{\dagger }{a}_{\omega }-{\hslash }{\omega }_{2}\displaystyle \int {\rm{d}}z\ {{ \mathcal E }}_{w}^{\dagger }(z,t){{ \mathcal E }}_{w}(z,t)\\ & & +\displaystyle \frac{N}{L}\displaystyle \int {\rm{d}}z\{-{\hslash }{{\rm{\Omega }}}_{{\rm{W}}}(t-z/c){\sigma }_{{eg}}(z,t){{\rm{e}}}^{{\rm{i}}{\omega }_{1}z/c}+{\rm{H}}{\rm{.}}{\rm{c}}.\\ & & -g{{ \mathcal E }}_{w}(z,t){\sigma }_{{es}}(z,t){{\rm{e}}}^{{\rm{i}}{\omega }_{2}z/c}+{\rm{H}}{\rm{.}}{\rm{c}}.\},\end{array}\end{eqnarray} \tag{ 29 }$$

where we have defined ${{ \mathcal E }}_{{\rm{w}}}(z,t)=\sqrt{\tfrac{L}{2\pi c}}{{\rm{e}}}^{{\rm{i}}{\omega }_{2}t}\int d\omega \ {a}_{\omega }{{\rm{e}}}^{{\rm{i}}(\omega -{\omega }_{2})z/c}$.

Assuming a real Rabi frequency Ω_W, this yields the Heisenberg–Langevin equations as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{{\rm{t}}}{\sigma }_{{se}} & = & -{\gamma }_{{es}}{\sigma }_{{se}}+{\rm{i}}{{\rm{\Omega }}}_{W}{{\rm{e}}}^{{\rm{i}}{\omega }_{1}z/c}{\sigma }_{{sg}}\\ & & -{\rm{i}}\bar{g}{{ \mathcal E }}_{w}{{\rm{e}}}^{{\rm{i}}{\omega }_{2}z/c}({\sigma }_{{ee}}-{\sigma }_{{ss}})+{F}_{{se}}\\ {\partial }_{{\rm{t}}}{\sigma }_{{sg}} & = & -{\gamma }_{0}{\sigma }_{{sg}}+{\rm{i}}{{\rm{\Omega }}}_{W}{{\rm{e}}}^{-{\rm{i}}{\omega }_{1}z/c}{\sigma }_{{se}}\\ & & -{\rm{i}}\bar{g}{{ \mathcal E }}_{w}{{\rm{e}}}^{{\rm{i}}{\omega }_{2}z/c}{\sigma }_{{eg}}+{F}_{{sg}}\\ {\partial }_{{\rm{t}}}{\sigma }_{{eg}} & = & -{\gamma }_{{eg}}{\sigma }_{{eg}}-{\rm{i}}{{\rm{\Omega }}}_{W}{{\rm{e}}}^{-{\rm{i}}{\omega }_{1}z/c}{\sigma }_{{gg}}+{F}_{{eg}},\end{array}\end{eqnarray} \tag{ 30 }$$

where ω₁ (ω₂) indicates the frequency of the write control field (write photon field), respectively, and ${\bar{g}}^{2}N=\tfrac{\bar{d}{\gamma }_{{es}}c}{L}$.

During the write process we account for possible depletion of the write laser intensity, and hence do not assume Ω_W(r, t) to be constant throughout the sample. As a result of the laser we create coherences between the $| g\rangle -| e\rangle$ transition, which forms the initial state for the write photon field. Here we proceed to find the atomic coherences prepared as a result of our exponential shaped resonant write control field.

For a sufficiently short write control field, the dynamics of the field and the atoms can be described with the dynamics along the $| g\rangle -| e\rangle$ transition. Ignoring the noise terms on σ_ge and making the analogy between the classical and quantum fields on the $| g\rangle -| e\rangle$ transition,

$$\begin{eqnarray}\begin{array}{rcl}c{\partial }_{z}{{\rm{\Omega }}}_{{\rm{W}}}(z,t) & = & {{\rm{i}}{g}}^{2}N{\sigma }_{{ge}}(z,t){{\rm{e}}}^{-{\rm{i}}{\omega }_{1}z/c},\\ {\partial }_{t}{\sigma }_{{ge}} & = & -{\gamma }_{{eg}}{\sigma }_{{ge}}+{\rm{i}}{{\rm{\Omega }}}_{{\rm{W}}}(z,t){{\rm{e}}}^{+{\rm{i}}{\omega }_{1}z/c}{\sigma }_{{gg}}\\ & \approx & -{\gamma }_{{eg}}{\sigma }_{{ge}}+{\rm{i}}{{\rm{\Omega }}}_{{\rm{W}}}(z,t){{\rm{e}}}^{+{\rm{i}}{\omega }_{1}z/c},\end{array}\end{eqnarray} \tag{ 31 }$$

where we have assumed that almost all atoms remain in the $| g\rangle$ level.

Let us first assume a write control field with Rabi frequency Ω_W that begins at t = 0. Taking the Laplace transforms from $t\to w$, we find

$$\begin{eqnarray}&&{\partial }_{z}{{\rm{\Omega }}}_{{\rm{W}}}(z,w)={\rm{i}}\displaystyle \frac{{g}^{2}N}{c}{\sigma }_{{ge}}(z,w){{\rm{e}}}^{-{\rm{i}}{\omega }_{1}z/c},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\sigma }_{{ge}}(z,w)=\displaystyle \frac{1}{w+{\gamma }_{{eg}}}[{\rm{i}}{{\rm{\Omega }}}_{{\rm{W}}}(z,w){{\rm{e}}}^{{\rm{i}}{\omega }_{1}z/c}+{\sigma }_{{ge}}(z,t=0)].\end{eqnarray} \tag{ 33 }$$

Insert equation (33) into equation (32), and use the initial condition σ_ge(z, t = 0) = 0 to obtain

$$\begin{eqnarray*}&&{\partial }_{z}{{\rm{\Omega }}}_{W}(z,w)=-\displaystyle \frac{{g}^{2}N}{c}\left(\displaystyle \frac{1}{w+{\gamma }_{{eg}}}\right){{\rm{\Omega }}}_{W}(z,w),\end{eqnarray*}$$

yielding

$$\begin{eqnarray*}&&{{\rm{\Omega }}}_{{\rm{W}}}(z,w)={{\rm{e}}}^{-\tfrac{{g}^{2}N}{c}\left(\tfrac{1}{w+{\gamma }_{{eg}}}\right)z}{{\rm{\Omega }}}_{{\rm{W}}}(z=0,w).\end{eqnarray*}$$

Insert this into equation (33) to obtain

$$\begin{eqnarray*}&&{\sigma }_{{ge}}(z,w)=({{\rm{i}}{\rm{e}}}^{{\rm{i}}{\omega }_{1}z/c})\left[\displaystyle \frac{1}{w+{\gamma }_{{eg}}}{{\rm{e}}}^{-\tfrac{{\bar{g}}^{2}N}{c}\displaystyle \frac{1}{w+{\gamma }_{{eg}}}z}{{\rm{\Omega }}}_{{\rm{W}}}(z=0,w)\right].\end{eqnarray*}$$

After inverting the Laplace transform, we now shift the limits to consider a write control field with support on negative times, giving

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{ge}}(z,t) & = & ({{\rm{i}}{\rm{e}}}^{{\rm{i}}{\omega }_{1}z/c}){\displaystyle \int }_{-\infty }^{t}{{\rm{e}}}^{-{\gamma }_{{eg}}(t-{t}_{1}^{{\prime\prime} })}\\ & & \times \,{J}_{0}\left[2\sqrt{\displaystyle \frac{{\gamma }_{{eg}}d}{L}(t-{t}_{1}^{{\prime\prime} })z}\right]{{\rm{\Omega }}}_{{\rm{W}}}(z=0,{t}_{1}^{{\prime\prime} }){{\rm{d}}{t}}_{1}^{{\prime\prime} },\end{array}\end{eqnarray} \tag{ 34 }$$

where J_n(x) indicates the nth Bessel function of the first kind.

Thus, with an exponential write control field ${{\rm{\Omega }}}_{{\rm{W}}}(0,t)={{\rm{\Omega }}}_{{\rm{W}}}^{\max }{{\rm{e}}}^{t/{\tau }_{{\rm{W}}}}$ sent up to t = 0, we evaluate the atomic coherence at t = 0 with the help of ${\int }_{0}^{\infty }{{\rm{e}}}^{-{At}}{J}_{0}[2\sqrt{{Bt}}]{\rm{d}}{t}=\tfrac{1}{A}{{\rm{e}}}^{-B/A}$ and finally obtain

$$\begin{eqnarray}&&{\sigma }_{{ge}}(z,0)={{\rm{e}}}^{{\rm{i}}{\omega }_{1}z/c}{\theta }_{0}{{\rm{e}}}^{-\tfrac{\alpha z}{2}},\end{eqnarray} \tag{ 35 }$$

where ${\theta }_{0}={\rm{i}}\tfrac{{{\rm{\Omega }}}_{\max }{\tau }_{{\rm{W}}}}{1+{\gamma }_{{eg}}{\tau }_{{\rm{W}}}}$ and $\alpha /2=d\tfrac{{\gamma }_{{eg}}{\tau }_{{\rm{W}}}}{1+{\gamma }_{{eg}}{\tau }_{{\rm{W}}}}\tfrac{1}{L}$.

After the preparation of atomic coherences, we begin to see spontaneous emission from the $| e\rangle$ level. Along with the field propagation equation, the relevant Heisenberg–Langevin equations are

$$\begin{eqnarray*}\begin{array}{rcl} & & c{\partial }_{z}{{ \mathcal E }}_{{\rm{w}}}={\rm{i}}\bar{g}{{N}{\rm{e}}}^{-{\rm{i}}{\omega }_{2}z/c}{\sigma }_{{se}}(z,t),\\ & & {\partial }_{t}{\hat{\sigma }}_{{es}}=-{\gamma }_{{es}}{\sigma }_{{se}}-{\rm{i}}\bar{g}{{ \mathcal E }}_{w}{{\rm{e}}}^{{\rm{i}}{\omega }_{2}z/c}({\sigma }_{{ee}}-{\sigma }_{{ss}})+{F}_{{se}}.\end{array}\end{eqnarray*}$$

Defining ${Q}^{\dagger }=\sqrt{N}{{\rm{e}}}^{-{\rm{i}}{\omega }_{2}z/c}{\sigma }_{{se}}$, we will consider the write emission for short detection times. Using (35) we thus replace σ_ee − σ_ss with its mean value at position z and t = 0 to obtain

$$\begin{eqnarray}\begin{array}{rcl}c{\partial }_{z}{{ \mathcal E }}_{{\rm{w}}}(z,t) & = & {\rm{i}}\bar{g}\sqrt{N}{Q}^{\dagger }(z,t),\\ {\partial }_{t}{Q}^{\dagger }(z,t) & = & -{\gamma }_{{es}}{Q}^{\dagger }(z,t)-{\rm{i}}\bar{g}\sqrt{N}| {\theta }_{0}{| }^{2}{{\rm{e}}}^{-\alpha z}{{ \mathcal E }}_{{\rm{w}}}(z,t)\\ & & +{F}_{Q}^{\dagger }(z,t).\end{array}\end{eqnarray} \tag{ 36 }$$

Performing first the Laplace transform in space ($z\to s$)

$$\begin{eqnarray*}&&\begin{array}{l}s{{ \mathcal E }}_{{\rm{w}}}(s,t)-{{ \mathcal E }}_{{\rm{w}}}(z=0,t)=A\ {Q}^{\dagger }(s,t),\\ {\partial }_{{\rm{t}}}{Q}^{\dagger }(s,t)=-{\gamma }_{{es}}{Q}^{\dagger }(s,t)+B{{ \mathcal E }}_{{\rm{w}}}(s+\alpha ,t)\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +{F}_{Q}^{\dagger }(s,t),\end{array}\end{eqnarray*}$$

and then in time ($t\to \omega$), we get

$$\begin{eqnarray*}\begin{array}{rcl} & & s{{ \mathcal E }}_{{\rm{w}}}(s,\omega )-{{ \mathcal E }}_{{\rm{w}}}(z=0,\omega )=A\ {Q}^{\dagger }(s,\omega ),\\ & & {Q}^{\dagger }(s,\omega )=\displaystyle \frac{1}{{\gamma }_{{es}}+\omega }B{{ \mathcal E }}_{{\rm{w}}}(s+\alpha ,\omega )+{F}_{Q}^{\dagger }(s,\omega )+{Q}^{\dagger }(s,t=0),\end{array}\end{eqnarray*}$$

where $A={\rm{i}}\tfrac{\bar{g}\sqrt{N}}{c}$ and $B=-{\rm{i}}\bar{g}\sqrt{N}{\theta }_{0}^{2}$.

Substituting the second line into the first, we eliminate Q(s, ω) and are left with a boundary term in Q:

$$\begin{eqnarray*}\begin{array}{rcl}s{ \mathcal E }(s,\omega )-{ \mathcal E }(z=0,\omega ) & = & \left(\displaystyle \frac{A}{{\gamma }_{{es}}+\omega }\right)[B{ \mathcal E }(s+\alpha ,\omega )\\ & & +{F}_{Q}^{\dagger }(s,\omega )+{Q}^{\dagger }(s,t=0)].\end{array}\end{eqnarray*}$$

The following formula also holds with a shift from s to s + α:

$$\begin{eqnarray*}&&\begin{array}{l}(s+\alpha ){{ \mathcal E }}_{{\rm{w}}}(s+\alpha ,\omega )-{{ \mathcal E }}_{{\rm{w}}}(z=0,\omega )\\ \quad =\,\left(\displaystyle \frac{A}{{\gamma }_{{es}}+\omega }\right)[B{{ \mathcal E }}_{{\rm{w}}}(s+2\alpha ,\omega )+{F}_{Q}^{\dagger }(s+\alpha ,\omega )+{Q}^{\dagger }(s+\alpha ,t=0)].\end{array}\end{eqnarray*}$$

By substituting ${{ \mathcal E }}_{{\rm{w}}}({s}+\alpha ,\omega )$ into the previous equation we can find ${ \mathcal E }(s,\omega )$ in terms of ${{ \mathcal E }}_{{\rm{w}}}({s}+2\alpha ,\omega )$, and by taking the substitution into the nth step we have

$$\begin{eqnarray*}\begin{array}{rcl}{{ \mathcal E }}_{{\rm{w}}}(s,\omega ) & = & K{(\omega )}^{n}D(n){{ \mathcal E }}_{{\rm{w}}}(s+n\alpha ,\omega )\\ & & +\displaystyle \frac{1}{B}\displaystyle \sum _{j=1}^{n}K{(\omega )}^{j}D(j){F}_{Q}^{\dagger }(s+(j-1)\alpha ,\omega )\\ & & +\displaystyle \frac{1}{B}\displaystyle \sum _{j=1}^{n}K{(\omega )}^{j}D(j){Q}^{\dagger }(s+(j-1)\alpha ,t^{\prime} =0)\\ & & +{[K(\omega )]}^{-1}\displaystyle \sum _{j=1}^{n}[K{(\omega )}^{j}D(j)]{{ \mathcal E }}_{{\rm{w}}}(z^{\prime} =0,\omega ),\end{array}\end{eqnarray*}$$

where $K(\omega )=\tfrac{{AB}}{{\gamma }_{{es}}+\omega }$ and $D(n)={\prod }_{k=0}^{n-1}\tfrac{1}{s\,+\,k\alpha }$.

Taking the limit of $n\to \infty$, the first term disappears, and we proceed to perform the inverse transform $s\to z$. With a shift in the index $j,{{ \mathcal L }}^{-1}[D(j+1)]=\tfrac{1}{j!}{\left(\tfrac{1-{{\rm{e}}}^{-\alpha z}}{\alpha }\right)}^{j}$ and the shifting property of the Laplace Transform,

$$\begin{eqnarray}\begin{array}{rcl} & & {{ \mathcal E }}_{w}(z,\omega )\\ & = & \displaystyle \frac{1}{B}\displaystyle \sum _{j=0}^{\infty }K{(\omega )}^{j+1}{\displaystyle \int }_{0}^{z}\displaystyle \frac{1}{j!}{\left(\displaystyle \frac{1-{{\rm{e}}}^{-\alpha (z-z^{\prime\prime} )}}{\alpha }\right)}^{j}{{\rm{e}}}^{-j\alpha z^{\prime\prime} }{F}_{Q}^{\dagger }(z^{\prime\prime} ,\omega )\ {\rm{d}}{z}^{\prime\prime} \\ & & +\displaystyle \frac{1}{B}\displaystyle \sum _{j=0}^{\infty }K{(\omega )}^{j+1}{\displaystyle \int }_{0}^{z}\displaystyle \frac{1}{j!}{\left(\displaystyle \frac{1-{{\rm{e}}}^{-\alpha (z-z^{\prime\prime} )}}{\alpha }\right)}^{j}{{\rm{e}}}^{-j\alpha z^{\prime\prime} }{Q}^{\dagger }(z^{\prime\prime} ,t=0)\ {\rm{d}}{z}^{\prime\prime} \\ & & +\displaystyle \frac{1}{K(\omega )}\displaystyle \sum _{j=0}^{\infty }K{(\omega )}^{j+1}\displaystyle \frac{1}{j!}{\left(\displaystyle \frac{1-{{\rm{e}}}^{-\alpha (z)}}{\alpha }\right)}^{j}{{ \mathcal E }}_{{\rm{w}}}(z^{\prime} =0,\omega )\\ & = & \displaystyle \frac{A}{{\gamma }_{{es}}+\omega }{\displaystyle \int }_{0}^{z}{{\rm{e}}}^{\left[\tfrac{1}{{\gamma }_{{es}}+\omega }M(z,z^{\prime\prime} ){{\rm{e}}}^{-\alpha z^{\prime\prime} }\right]}{F}_{Q}^{\dagger }(z^{\prime\prime} ,\omega ){\rm{d}}{z}^{\prime\prime} \\ & & +\displaystyle \frac{A}{{\gamma }_{{es}}+\omega }{\displaystyle \int }_{0}^{z}{{\rm{e}}}^{\left[\tfrac{1}{{\gamma }_{{es}}+\omega }M(z,z^{\prime\prime} ){{\rm{e}}}^{-\alpha z^{\prime\prime} }\right]}{Q}^{\dagger }(z^{\prime\prime} ,t^{\prime} =0){\rm{d}}{z}^{\prime\prime} \\ & & +{{\rm{e}}}^{\tfrac{1}{{\gamma }_{{es}}+\omega }M(z,0)}{{ \mathcal E }}_{w}(z=0,\omega ),\end{array}\end{eqnarray} \tag{ 37 }$$

where $M(z^{\prime} ,z^{\prime\prime} )=\tfrac{{AB}}{\alpha }[1-{{\rm{e}}}^{-\alpha (z^{\prime} -z^{\prime\prime} )}]$.

Finally, noting that

$$\begin{eqnarray*}\begin{array}{rcl} & & {{ \mathcal L }}^{-1}\left[\displaystyle \frac{1}{{\gamma }_{{es}}+\omega }{{\rm{e}}}^{\tfrac{A}{{\gamma }_{{es}}+\omega }}\right]={{\rm{e}}}^{-{\gamma }_{{es}}t}[{\sqrt{{At}}}^{(-1)}{I}_{1}(2\sqrt{{At}})+{I}_{2}(2\sqrt{{At}})],\\ & & {{ \mathcal L }}^{-1}[{{\rm{e}}}^{\tfrac{A}{{\gamma }_{{es}}+\omega }}]={{\rm{e}}}^{-{\gamma }_{{es}}t}\left[\sqrt{\displaystyle \frac{A}{t}}{I}_{1}(2\sqrt{{At}})+\delta (t)\right],\end{array}\end{eqnarray*}$$

we get

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal E }}_{w}(z,t) & = & A{\displaystyle \int }_{0}^{z}{\displaystyle \int }_{0}^{t}{{\rm{e}}}^{-{\gamma }_{{es}}(t-{t}_{1}^{{\prime\prime} })}{H}_{1}[\alpha ,z,{z}_{1}^{{\prime\prime} },t,{t}_{1}^{{\prime\prime} }]{F}_{Q}^{\dagger }({z}_{1}^{{\prime\prime} },{t}_{1}^{{\prime\prime} }){{\rm{d}}{t}}_{1}^{{\prime\prime} }{{\rm{d}}{z}}_{1}^{{\prime\prime} }\\ & & +A{\displaystyle \int }_{0}^{z}{{\rm{e}}}^{-{\gamma }_{{es}}(t)}{H}_{1}[\alpha ,z,{z}_{1}^{{\prime\prime} },t,0]{Q}^{\dagger }({z}_{1}^{{\prime\prime} },0){{\rm{d}}{z}}_{1}^{{\prime\prime} }\\ & & +{\displaystyle \int }_{0}^{t}{{\rm{e}}}^{-{\gamma }_{{es}}(t-t^{\prime\prime} )}{H}_{2}(\alpha ,z,0,t,t^{\prime\prime} ){{ \mathcal E }}_{w}(0,t^{\prime\prime} ){\rm{d}}{t}^{\prime\prime} \\ & & +{{ \mathcal E }}_{w}(0,t),\end{array}\end{eqnarray} \tag{ 38 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{H}_{1}(\alpha ,{z}_{1},{z}_{2},{t}_{1},{t}_{2}) & = & {I}_{0}[2\sqrt{M({z}_{1},{z}_{2}){{\rm{e}}}^{-\alpha {z}_{2}}({t}_{1}-{t}_{2})}],\\ {H}_{2}(\alpha ,{z}_{1},{z}_{2},{t}_{1},{t}_{2}) & = & \sqrt{\displaystyle \frac{M({z}_{1},{z}_{2})}{{t}_{1}-{t}_{2}}}{I}_{1}[2\sqrt{M({z}_{1},{z}_{2}){{\rm{e}}}^{-\alpha {z}_{2}}({t}_{1}-{t}_{2})}].\end{array}\end{eqnarray*}$$

Computing the photon flux requires the commutation relations for Q and a 2-point noise correlation function involving F_Q. In a short time window τ_d where σ_ee − σ_ss is not changing, and with the Einstein relations (see Ch 15.5 of [23]), the Langevin equations for system operators can be written

$$\begin{eqnarray}&&{\dot{A}}_{\mu }={D}_{\mu }(t)+{F}_{\mu }(t).\end{eqnarray} \tag{ 39 }$$

The corresponding memoryless noise correlations for operators μ and ν are such that

$$\begin{eqnarray}&&\langle {F}_{\mu }(t^{\prime} ){F}_{\nu }(t^{\prime\prime} )\rangle =2\langle {D}_{\mu \nu }\rangle \delta (t^{\prime} -t^{\prime\prime} ),\end{eqnarray} \tag{ 40 }$$

where

$$\begin{eqnarray}&&2\langle {D}_{\mu \nu }\rangle =-\langle {A}_{\mu }{D}_{\nu }\rangle -\langle {D}_{\mu }{A}_{\nu }\rangle +\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\langle {A}_{\mu }{A}_{\nu }\rangle .\end{eqnarray} \tag{ 41 }$$

Thus, identifying terms in equation (36) with terms in equation (39), we make use of

$$\begin{eqnarray}\begin{array}{rcl}[Q(z,t),{Q}^{\dagger }(z^{\prime} ,t)] & = & N[{\sigma }_{{\rm{e}}{s}}(z,t),{\sigma }_{{se}}(z^{\prime} ,t)]\\ & = & L\delta (z-z^{\prime} )| {\theta }_{0}{| }^{2}{{\rm{e}}}^{-\alpha z^{\prime} },\end{array}\end{eqnarray} \tag{ 42 }$$

then we make use of the fact that $\langle {Q}^{\dagger }(z,t)Q(z^{\prime} ,t)\rangle$ right after our preparation of atomic coherences is zero, giving $\langle Q(z,t){Q}^{\dagger }(z^{\prime} ,t)\rangle =L\delta (z-z^{\prime} )| {\theta }_{0}{| }^{2}{{\rm{e}}}^{-\alpha z^{\prime} }$.

Then one obtains

$$\begin{eqnarray}&&2\langle {D}_{Q,{Q}^{\dagger }}\rangle =2{\gamma }_{{es}}L| {\theta }_{0}{| }^{2}{{\rm{e}}}^{-\alpha z}\delta (z-z^{\prime} ),\end{eqnarray} \tag{ 43 }$$

yielding

$$\begin{eqnarray}&&\langle {F}_{Q}(z,t){F}_{Q}^{\dagger }(z^{\prime} ,t^{\prime} )\rangle =2{\gamma }_{{es}}L| {\theta }_{0}{| }^{2}{{\rm{e}}}^{-\alpha z}\delta (z-z^{\prime} )\delta (t-t^{\prime} ),\end{eqnarray} \tag{ 44 }$$

valid when σ_ee − σ_ss is not changing.

This yields a photon flux of

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{c}{L}\langle {{ \mathcal E }}_{{\rm{w}}}^{\dagger }(L,t){{ \mathcal E }}_{{\rm{w}}}(L,t)\rangle & = & \displaystyle \frac{c}{L}\displaystyle \frac{{\bar{g}}^{2}N}{{c}^{2}}{\displaystyle \int }_{0}^{L}{{\rm{e}}}^{-2{\gamma }_{{es}}t}{H}_{1}{[\alpha ,L,{z}_{1}^{{\prime\prime} },t,0]}^{2}V| {\theta }_{0}{| }^{2}{{\rm{e}}}^{-\alpha {z}_{1}^{{\prime\prime} }}{{\rm{d}}{z}}_{1}^{{\prime\prime} }\\ & & +\displaystyle \frac{c}{L}\displaystyle \frac{{\bar{g}}^{2}N}{{c}^{2}}{\displaystyle \int }_{0}^{t}{\displaystyle \int }_{0}^{L}{{\rm{e}}}^{-2{\gamma }_{{es}}(t-{t}_{1}^{{\prime\prime} })}{H}_{1}{[\alpha ,L,{z}_{1}^{{\prime\prime} },t,{t}_{1}^{{\prime\prime} }]}^{2}\\ & & \times \,2{\gamma }_{{es}}L| {\theta }_{0}{| }^{2}{{\rm{e}}}^{-\alpha {z}_{1}^{{\prime\prime} }}{{\rm{d}}{z}}_{1}^{{\prime\prime} }{{\rm{d}}{t}}_{1}^{{\prime\prime} }.\end{array}\end{eqnarray*}$$

For sufficiently short detection times ${\tau }_{{\rm{d}}}\ll \tfrac{1}{2{\gamma }_{{es}}}$, the noise contribution (second term) can be ignored, and furthermore when the photon number is much smaller than 1 $({\tau }_{{\rm{d}}}\ll {\left\{\tfrac{{\bar{g}}^{2}N}{c}| {\theta }_{0}{| }^{2}\left(\tfrac{1-{{\rm{e}}}^{-\alpha L}}{\alpha }\right)\right\}}^{-1})$ we can consider just the leading term in the series expansion, and observe a constant flux

$$\begin{eqnarray}\begin{array}{rcl} & & \displaystyle \frac{c}{L}\langle {{ \mathcal E }}_{{\rm{w}}}^{\dagger }(L,{\tau }_{{\rm{d}}}){{ \mathcal E }}_{{\rm{w}}}(L,{\tau }_{{\rm{d}}})\rangle \\ & = & \displaystyle \frac{{\bar{g}}^{2}N}{c}| {\theta }_{0}{| }^{2}{\displaystyle \int }_{0}^{L}{({I}_{0}[2\sqrt{M[L,{z}_{1}^{{\prime\prime} }]{{\rm{e}}}^{-\alpha {z}_{1}^{{\prime\prime} }}{\tau }_{{\rm{d}}}}])}^{2}{{\rm{e}}}^{-\alpha {z}_{1}^{{\prime\prime} }}{{\rm{d}}{z}}_{1}^{{\prime\prime} }\\ & \approx & \displaystyle \frac{{\bar{g}}^{2}N}{c}| {\theta }_{0}{| }^{2}{\displaystyle \int }_{0}^{L}(1+2M[L,{z}_{1}^{{\prime\prime} }]{{\rm{e}}}^{-\alpha {z}_{1}^{{\prime\prime} }}{\tau }_{{\rm{d}}}+O({\tau }_{{\rm{d}}}^{2})){{\rm{e}}}^{-\alpha {z}_{1}^{{\prime\prime} }}{{\rm{d}}{z}}_{1}^{{\prime\prime} }\\ & = & \displaystyle \frac{{\bar{g}}^{2}N}{c}| {\theta }_{0}{| }^{2}\displaystyle \frac{1-{{\rm{e}}}^{-\alpha L}}{\alpha }.\end{array}\end{eqnarray} \tag{ 45 }$$

We therefore obtain a photon number of $\tfrac{{\bar{g}}^{2}N}{c}| {\theta }_{0}{| }^{2}\tfrac{1-{{\rm{e}}}^{-\alpha L}}{\alpha }{\tau }_{{\rm{d}}}$ within this short detection window. This is precisely the excited atom fraction multiplied by $\bar{d}{\gamma }_{{es}}{\tau }_{{\rm{d}}}$, since the fraction of atoms that were excited after the write pulse is $\tfrac{1}{L}{\int }_{0}^{L}\langle {\sigma }_{{ee}}(z,0)\rangle {\rm{d}}{z}=\tfrac{1}{L}{\int }_{0}^{L}| {\theta }_{0}{| }^{2}{{\rm{e}}}^{-\alpha z}{\rm{d}}{z}=| {\theta }_{0}{| }^{2}\tfrac{1-{{\rm{e}}}^{-\alpha L}}{\alpha L}$.

We start with the description of the spin operator from equation (30), by defining $S=\sqrt{N}{\sigma }_{{gs}}{{\rm{e}}}^{-{\rm{i}}({\omega }_{1}-{\omega }_{2})z/c}$ and replacing σ_eg(z, t) by its mean value ${\theta }_{0}^{* }{{\rm{e}}}^{-\alpha z/2}{{\rm{e}}}^{-{\rm{i}}{\omega }_{1}z/c}$

$$\begin{eqnarray*}\begin{array}{rcl}({\partial }_{{\rm{t}}}+{\gamma }_{0}){\hat{S}}^{\dagger }(z,t)-{F}_{S}^{\dagger }(z,t) & = & -{\rm{i}}\bar{g}\sqrt{N}{{\rm{e}}}^{{{\rm{i}}{k}}_{{\rm{w}}}{\rm{r}}}{{ \mathcal E }}_{{\rm{w}}}(z,t){\sigma }_{{eg}}\\ & = & -{\rm{i}}\bar{g}\sqrt{N}{{ \mathcal E }}_{{\rm{w}}}(z,t)[{\theta }_{0}^{* }{{\rm{e}}}^{-\alpha z/2}].\end{array}\end{eqnarray*}$$

Take the Laplace transform from $t\to \omega$ to see

$$\begin{eqnarray*}&&(\omega +{\gamma }_{0}){S}^{\dagger }(z,\omega )-{S}^{\dagger }(z,t=0)-{F}_{S}^{\dagger }(z,\omega )=C(z){{ \mathcal E }}_{w}(z,\omega ),\end{eqnarray*}$$

where $C(z)=-{\rm{i}}\bar{g}\sqrt{N}{\theta }_{0}^{* }{{\rm{e}}}^{-\alpha z/2}$. Then we have

$$\begin{eqnarray}&&{S}^{\dagger }(z,\omega )=\displaystyle \frac{C(z)}{\omega +{\gamma }_{0}}{{ \mathcal E }}_{{\rm{w}}}(z,\omega )+\displaystyle \frac{1}{\omega +{\gamma }_{0}}{S}^{\dagger }(z,t=0)+\displaystyle \frac{1}{\omega +{\gamma }_{0}}{F}_{S}^{\dagger }(z,\omega ),\end{eqnarray} \tag{ 46 }$$

and noting that ${{ \mathcal L }}^{-1}\left[\tfrac{1}{\omega +{\gamma }_{0}}\right]={{\rm{e}}}^{-{\gamma }_{0}t}$ yields

$$\begin{eqnarray}\begin{array}{rcl}{S}^{\dagger }(z,t) & = & C(z){\displaystyle \int }_{0}^{t}{{\rm{e}}}^{-{\gamma }_{0}(t-t^{\prime} )}{{ \mathcal E }}_{w}(z,t^{\prime} ){\rm{d}}{t}^{\prime} \\ & & +{{\rm{e}}}^{-{\gamma }_{0}t}{S}^{\dagger }(z,t=0)+{\displaystyle \int }_{0}^{t}{{\rm{e}}}^{-{\gamma }_{0}(t-t^{\prime} )}{F}_{S}^{\dagger }(z,t^{\prime} ){\rm{d}}{t}^{\prime} ,\end{array}\end{eqnarray} \tag{ 47 }$$

where the field expression ${ \mathcal E }$ from the previous subsection is required. Ignoring terms that do not show up in the normal ordered $\langle {S}^{\dagger }S\rangle$, we have

$$\begin{eqnarray}\begin{array}{rcl}{S}^{\dagger }(z,t) & = & C(z){\displaystyle \int }_{0}^{t}{\rm{d}}{t}^{\prime} {{\rm{e}}}^{-{\gamma }_{0}(t-t^{\prime} )}\\ & & \times \,\left({\displaystyle \int }_{0}^{t^{\prime} }{\rm{d}}{t}^{\prime\prime} {{\rm{e}}}^{-{\gamma }_{{es}}(t^{\prime} -t^{\prime\prime} )}{H}_{2}(\alpha ,z,0,t^{\prime} ,t^{\prime\prime} ){{ \mathcal E }}_{{\rm{w}}}(0,t^{\prime\prime} )+{{ \mathcal E }}_{{\rm{w}}}(0,t^{\prime} )\right).\end{array}\end{eqnarray} \tag{ 48 }$$

Computing $\langle {S}^{\dagger }S\rangle$ requires the commutator $[{{ \mathcal E }}_{{\rm{w}}}(z,t),{{ \mathcal E }}_{{\rm{w}}}^{\dagger }(z^{\prime} ,t^{\prime} )]=L\delta [z-z^{\prime} -c(t-t^{\prime} )]$ and yields 4 terms. In the short time window where one can ignore the atomic dephasing $\left({\tau }_{{\rm{d}}}\ll \tfrac{1}{2{\gamma }_{0}},\tfrac{1}{2{\gamma }_{{es}}}\right)$, and also where the photon number is much smaller than 1 $\left({\tau }_{{\rm{d}}}\ll {\left\{\tfrac{{\bar{g}}^{2}N}{c}| {\theta }_{0}{| }^{2}\left(\tfrac{1-{{\rm{e}}}^{-\alpha L}}{\alpha }\right)\right\}}^{-1}\right)$, only one term dominates (the term independent of H₂). The number of spins is then equivalent to the photon number

$$\begin{eqnarray*}&&\displaystyle \frac{1}{L}{\int }_{0}^{L}\langle {S}^{\dagger }(z,{\tau }_{{\rm{d}}})S(z,{\tau }_{{\rm{d}}})\rangle {\rm{d}}{z}\approx \displaystyle \frac{{\bar{g}}^{2}N}{c}| {\theta }_{0}{| }^{2}\displaystyle \frac{1-{{\rm{e}}}^{-\alpha L}}{\alpha }{\tau }_{{\rm{d}}}.\end{eqnarray*}$$

To ensure that a single photon source is single mode in all degrees of freedom, one would have to verify that the outgoing emission is not produced in a combination of modes. In the system that we consider, one cause for multi-mode emission are multiple two mode squeezing processes occurring during the initial write process. Thereafter, the subsequent retrieval process yields a read photon in more than one mode. Here, we include a short section to explain how the unconditional autocorrelation measurement (g⁽²⁾(0)) scales with the number of driven mode pairs [20], and this allows one to verify that no higher-number two mode squeezing processes have occured.

Let us first consider the state created by K vacuum squeezing processes

$$\begin{eqnarray}&&{\rho }_{{\rm{multi}}}={(1-p)}^{K}\left[{{\rm{e}}}^{\sqrt{p}{\sum }_{m=1}^{K}{a}_{m}^{\dagger }{b}_{m}^{\dagger }}| {\rm{\Omega }}\rangle \langle {\rm{\Omega }}| {{\rm{e}}}^{\sqrt{p}{\sum }_{n=1}^{K}{a}_{n}{b}_{n}}\right],\end{eqnarray} \tag{ 49 }$$

where $| {\rm{\Omega }}\rangle$ indicates the vacuum in all modes. Now consider a detector that sees all the K modes, giving a number operator of the form

$$\begin{eqnarray}&&{\hat{N}}_{K}=\displaystyle \sum _{c=1}^{K}{a}_{c}^{\dagger }{a}_{c}.\end{eqnarray} \tag{ 50 }$$

This yields

$$\begin{eqnarray}&&\mathrm{Tr}\left[\displaystyle \sum _{c,d=1}^{K}{a}_{c}^{\dagger }{a}_{d}^{\dagger }{a}_{c}{a}_{d}{\rho }_{{\rm{multi}}}\right]={\left(\displaystyle \frac{p}{1-p}\right)}^{2}({K}^{2}+K),\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}&&\mathrm{Tr}\left[\displaystyle \sum _{c=1}^{K}{a}_{c}^{\dagger }{a}_{c}{\rho }_{{\rm{multi}}}\right]=\displaystyle \frac{p}{1-p}K,\end{eqnarray} \tag{ 52 }$$

giving the unconditional autocorrelation of the a modes

$$\begin{eqnarray}\begin{array}{rcl}{g}_{{\rm{multi}}}^{(2)} & = & \displaystyle \frac{\left\langle {\displaystyle \sum }_{c,d=1}^{K}{a}_{c}^{\dagger }{a}_{d}^{\dagger }{a}_{c}{a}_{d}\right\rangle }{{\left\langle {\displaystyle \sum }_{c=1}^{K}{a}_{c}^{\dagger }{a}_{c}\right\rangle }^{2}}\\ & = & 1+\displaystyle \frac{1}{K},\end{array}\end{eqnarray} \tag{ 53 }$$

which leads us to check if our photon field is indeed consistent with that coming from a single-pair two-mode squeezing process.

We thus compute the unconditional autocorrelation function of the read photon field at time t

$$\begin{eqnarray*}&&{g}_{{\rm{read}}}^{2}=\displaystyle \frac{\langle {{ \mathcal E }}_{r}^{\dagger }(0,t){{ \mathcal E }}_{r}^{\dagger }(0,t){{ \mathcal E }}_{r}(0,t){{ \mathcal E }}_{r}(0,t)\rangle }{\langle {{ \mathcal E }}_{r}^{\dagger }(0,t){{ \mathcal E }}_{r}(0,t){\rangle }^{2}},\end{eqnarray*}$$

In the regime we consider, where we have a short detection time and a fast readout, developing the numerator of the g⁽²⁾ function leads to the term

$$\begin{eqnarray*}\begin{array}{rcl} & & \langle {{ \mathcal E }}_{{\rm{w}}}(0,{t}_{a}){{ \mathcal E }}_{{\rm{w}}}(0,{t}_{b}){{ \mathcal E }}_{{\rm{w}}}^{\dagger }(0,{t}_{c}){{ \mathcal E }}_{{\rm{w}}}^{\dagger }(0,{t}_{{\rm{d}}})\rangle \\ & = & \langle {{ \mathcal E }}_{{\rm{w}}}(0,{t}_{a})\left[{{ \mathcal E }}_{{\rm{w}}}^{\dagger }(0,{t}_{c}){{ \mathcal E }}_{{\rm{w}}}(0,{t}_{b})+\displaystyle \frac{L}{c}\delta ({t}_{b}-{t}_{c})\right]{{ \mathcal E }}_{{\rm{w}}}(0,{t}_{{\rm{d}}})\rangle \\ & = & {\left(\displaystyle \frac{L}{c}\right)}^{2}\delta ({t}_{a}-{t}_{c})\delta ({t}_{b}-{t}_{{\rm{d}}})+{\left(\displaystyle \frac{L}{c}\right)}^{2}\delta ({t}_{a}-{t}_{{\rm{d}}})\delta ({t}_{b}-{t}_{c}),\end{array}\end{eqnarray*}$$

which yields

$$\begin{eqnarray*}&&{g}_{{\rm{read}}}^{2}=\displaystyle \frac{2\langle {{ \mathcal E }}_{r}^{\dagger }(0,t){{ \mathcal E }}_{r}(0,t){\rangle }^{2}}{\langle {{ \mathcal E }}_{r}^{\dagger }(0,t){{ \mathcal E }}_{r}(0,t){\rangle }^{2}}=2\end{eqnarray*}$$

as it should, since we assume a mono-mode emission (K = 1). Here we have used equation (25) and the leading term of equation (48).