B J Lemaire

Water quality in urban lakes: from continuous monitoring to forecasting Application to cyanobacteria dynamics in Lake Enghien (France)

PROLIPHYC Project

The Proliphyc project has developed a continuous in-situ monitoring system for cyanobacteria in freshwater ecosystems. It consists in a measurement buoy equipped with meteorological sensors and immersed probes to measure water quality parameters (see Fig. 1). This system is able to provide a long-term, high-frequency monitoring of lakes and reservoirs. Data set obtained can be used in order (Le Vu et al., 2010):

Study site: Lake of Enghien-les-Bains

In 2009, such monitoring system was implemented in Lake Enghien-les-Bains (Paris suburbs, France, see Fig. 2 and Table 1). Lake Enghien is an urban shallow lake that plays a significant role in the stormwater management of its watershed by storing up to 100,000m 3 of rainwater. The lake receives wastewater discharges from inappropriate connections in the stormwater network. This input results in a deterioration of the water quality and the lake is frequently affected by cyanobacteria blooms of Planktothrix agardhii, a potentially toxic cyanobacterium (Quiblier et al. 2008).

Water quality indicators

The time series can be used to infer indicators of cyanobacteria biomass, useful for lake management strategies. Three indicators, built from the raw time series, were proposed for Lake Enghien (see Fig. 3):

-Water temperature and oxygen saturation rate associated to daily maximal cyanobacteria concentration;

-Cyanobacteria daily variation rate

Cyanobacteria modelling

The coupled model DYRESM-CAEDYM (DYCD) was used for deterministic simulation of cyanobacteria dynamics in Lake Enghien. DYRESM is a one-dimensional numerical model for predicting the vertical distribution of temperature in lakes and reservoirs (Imerito 2007).

It was coupled to CAEDYM, the aquatic ecosystem model to simulate cyanobacteria dynamics [START_REF] Hamilton | Prediction of water quality in lakes and reservoirs. Part I -Model description[END_REF] 2) To collect surveillance data series to observe the general patterns of the aquatic ecosystem and to assess longterm changes;

3) To feed a statistical short-term forecasting model for early warning of cyanobacteria blooms; The results of both models showed good agreement with observed values. The NNM performances benefited from the high frequency of the measurements. As it is easier to implement, it provides for lake managers an affordable support for anticipating cyanobacteria blooms.

Conclusion

This continuous in-situ monitoring and forecasting system for cyanobacteria in freshwater ecosystems provides high-frequency time series, compulsory for addressing research issues. It allows researchers to implement deterministic models of phytoplankton dynamics and stakeholders to assess long-term changes of the aquatic ecosystem functioning. Moreover, the measured data can be continuously processed to infer daily, seasonal or annual water quality indicators of the lake. Finally, a statistical short-term forecasting model can supply a prevision of the cyanobacteria biomass and possibly an early warning of blooms. The indicators and the forecasts can be displayed within an information system easily reached by the citizens and/or the stakeholders.

Fig. 1 :

 1 Fig. 1 : PROLIPHYC buoy measured variables (*Chl-a total and 4 phytoplankton groups)

Fig. 2 :

 2 Fig. 2 : Lake Enghien: (a) Location in France, (b) Aerial vue of the lake and its neighbourhood (IAURIF 2008) Watershed Lake

Fig. 3 :

 3 Fig. 3 : Raw data: (a) Chlorophyll-a concentrations (b) Water temperature and O 2 concentration. Water quality indicators: (c) Cyanobacteria daily variation rate (d) Water temperature and O 2 saturation related to cyanobacteria daily maxima.

Fig. 4 :

 4 Fig. 4 : Dyresm-Caedym model structure (* Data from Proliphyc buoy)

( 1 )

 1 Fig. 6 : NNM and DYCD cyanobacteria forecasting in 2009: (a) 04-07 July (b) 11-14 July

Table 1 :

 1 Lake Enghien characteristics

Acknowledgements: This work is part of the Proliphyc project (Monitoring System of Phytoplankton Blooms -Application to Cyanobacteria) funded by the ANR-PRECODD program. We wish to acknowledge SIARE and Enghien municipality for their support.

Contact : silvat@leesu.enpc.fr bvl@leesu.enpc.fr


Modelling of cyanobacteria blooms in urban lakes

Application to Lake of Enghien-les-Bains

Context

The ecosystems in urban lakes are very vulnerable to human pressure because of their specific physical functioning, characterized by intermittent thermal stratification and low flow speeds. Furthermore, the high intake of nutrients in these environments encourages blooms of phytoplankton, including toxic cyanobacteria, disrupting their use and causing health problems (Catherine et al., 2008). This work focuses on the implementation of a predictive model of cyanobacteria blooms in urban lakes. It is part of a project whose objective is to develop a monitoring and warning system, in real-time, of phytoplankton blooms in freshwater ecosystems.

Silva T. (1,*) , Vinçon-Leite B. (1) , Lemaire B.J. (1) , Le Vu B. (1,2) , Escoffier N. (3) , Prévot F. (4) , Quiblier C. (3) ,Tassin B. (1) (*) silvat@leesu.enpc.fr

Lake of Enghien-les-Bains

Lake of Enghien-les-Bains is located in Val-d'Oise, France (48°58'N, 2°18'E) (Fig. 1). It is an urban shallow lake (maximum depth: 2.5 m; mean depth: 1.3 m; area: 41 ha; volume: 534 000 m 3 ), it plays a significant role in the stormwater management of its watershed (54 km 2 ; over 200 000 inhabitants) by storing up to 100000m 3 of rainwater. The lake receives wastewater discharges from inappropriate connections in the stormwater network. This input results in a deterioration of the water quality and the lake is frequently affected by cyanobacteria blooms (Quiblier et al. 2008). In July 2009, the cyanobacterial biomass has reached 300 µgChlaL -1 (Silva, 2010). The prevailing species, Planktothrix agardhii, is potentially toxic and is usually found in the first meters of the water column in shallow turbid eutrophic lakes.

First results

The DYRESM model was run for the period from 1 st April to 17 th August 2009, the hottest period of the year and therefore, more favorable to cyanobacteria blooms (Fig. 4A). The maximum and minimum layer thickness was defined, respectively 0.16m and 0.70 m. The modelling of the cyanobacteria dynamics was conducted during a shorter period, from 1 st to 20 th June 2009 when cyanobacteria are the prevailing phytoplankton group (Fig. 4B).

Discussion

The thermal model was able to adequately represent the daily cycles and the seasonal changes in water temperature (Fig 5A). Moreover, despite the shallowness of the lake, a thermal stratification at a depth of 2.0 m was noted during the summer (Fig 5B ).

Concerning the biological simulation, the negative value of the arithmetic mean indicates that the model underestimates the phytoplankton biomass (Fig 6A). It has been stated in several studies on cyanobacteria that nutrient concentration is very significant to phytoplankton growth (Quiblier et al. 2008). In an urban lake, the stormwater discharge may play a role in triggering the proliferation of phytoplankton because nutrient load can rapidly increase during rainfall events. Therefore, the

DYRESM-CAEDYM

DYRESM is a one-dimensional numerical model for predicting the vertical distribution of temperature and density in lakes and reservoirs. It is based on a Lagrangian layer scheme (Imerito 2007). DYRESM was coupled to CAEDYM, the aquatic ecosystem model to simulate the cyanobacteria dynamics (Hamilton and Schladow, 1997). The input data are the lake morphometry, inflows, outflows, meteorological forcing (wind speed, air temperature, solar radiation, rainfall, cloud cover and vapor pressure), and the initial conditions for all variables to model (Fig. 2).

Measurement buoy

The data required for modelling (except the lake morphometry, flows and the nutrients) were obtained from a freshwater-adapted measurement buoy installed on the lake in

France

Ile de France Knowledge about nutrient inputs is necessary to improve the model results, specially during rain events. To broaden the simulation period, other groups of phytoplankton and zooplankton community should be taken into account because of the competitive relations to be imposed on cyanobacteria. Model validation on a data set different from that used for calibration will make it usable as a tool for predicting short and medium term cyanobacteria blooms in urban lakes such as Lake of Enghien-les-Bains.

Conclusion and perspectives

These results show the advantage of coupling this type of model with high frequency measurements, to simulate the phytoplankton biomass in a shallow urban lake. The model showed a good performance in simulating the cyanobacteria dynamics in an urban shallow lake. Moreover, DYRESM-CAEDYM can be employed as a predictive model of phytoplankton blooms based on the chlorophyll and nutrient concentrations and on the weather forecast.

Acknowledgements: This work is part of the Proliphyc project (Monitoring System of Phytoplankton Blooms -Application to Cyanobacteria) funded by the ANR-PRECODD were obtained from a freshwater-adapted measurement buoy installed on the lake in November 2008. This buoy was located at the deepest region of the lake and it was equipped with meteorological sensors and immersed measuring probes of water quality parameters (Fig. 3). The measurements were performed in high frequency, every 30 min, at a depth varying from 0.50 to 1.0 m (for water parameters). Data were transmitted through GPRS protocol to a database as a daily email. If a shorter time step was necessary, the buoy could be remotely queried, which allows the continuous monitoring in a real-time of the physical-biological conditions of the lake.

The lake morphometry was obtained from studies carried out for the lake management by the Syndicat Intercommunal d'Assainissement de la Région d'Enghien-les-Bains (SIARE). The 2009 flows were not available, however, during dry weather, the inflow is approximately 8 100 m 3 day -1 (SIARE, 2004). Using a water balance, the outflow was set equal to 5 000 m 3 day -1 . Subsequently, in order to verify the reliability of the results, the model sensitivity to flows was checked.

Fig. 4 :

 4 Fig. 4 : Buoy measurements: water temperature (A) and phytoplankton biomass (B)

Fig. 5 :Fig. 1 :Fig

 51 Fig. 5 : Model results versus buoy measurements for water temperature (A) and water temperature profile simulated (B).

  Flow variation causes only insignificant changes in the model result. The nutrient concentrations (N and P) were obtained from a monthly monitoring of the lake physical-chemical parameters conducted from November 2007 to June 2008 (Marchandise 2008).

Fig. 6 :Fig. 3 :

 63 Fig. 6 : Model results versus buoy measurements for phytoplankton biomass (A) and rainfall during the simulation period (B).


PROLIPHYC PROJECT A measuring buoy

High-frequency data (15-30 min) Data teletransmission (daily)

Warning system

Proliphyc Project

5

Warning system

Based on the weather forecast and buoy measurements Short-term predictive modelling (1 week)

50 100 150 200 250 300 350 400 Chl-a concentration (ug L -1 ) Total chl-a Cyanobacteria chl-a

				Lake Enghien Lake Enghien Lake Enghien
	LAKE ENGHIEN-LES-BAINS MODELLING DYCD RESULTS
	41 ha, mean depth = 1.5 m Rainwater storing, sport, landscape Model results show good agreement with measurements DYRESM-CAEDYM model (Centre for Water Research, High frequency data thorough assesment of results University of Western Australia) Gal et al. (2009); Burger et al. (2008); Hornung (2002)
	Inappropriate wastewater discharges 1 phytoplankton group: cyanobacterium Planktothrix agardhii Medium-term simulations, hourly time step
	Planktothrix agardhii blooms Model calibration and validation
		Calibration period 1-
	450
	400	Calibration			Validation
				(IAURIF 2008)	(Catherine Quiblier 2008)
				(IAURIF 2008)
						(Escoffier 2008)
						6 7
	19/05 0	03/06	18/06	03/07	18/07	02/08

http://leesu.univ http://leesu.univ http://leesu.univ http://leesu.univ----paris paris paris paris----est.fr/proliphyc est.fr/proliphyc est.fr/proliphyc est.fr/proliphyc

16 June 2009 Validation period 2-16 July 2009 8 r = 0.92 rmse = 0.9°C r = 0.96 rmse = 15.6 µgL -1 r = 0.74 rmse = 19.8 µgL -1 r = 0.97 rmse = 0.6°C

								Lake Pampulha Lake Pampulha Lake Pampulha	Conclusions
	LAKE PAMPULHA WATERSHED AND LAKE MONITORING LAKE ECOLOGICAL MODEL CONCLUSIONS
	Belo Horizonte, Minas Gerais, Brazil 197 ha, 10 hm 3 , 5,0 m (mean depth) Belo Horizonte municipality data Lake Enghien High frequency data* Cyanobacteria species: Warning system Two main tributaries and lake outlet : Nostocaceae, Microcystis wesenbergii, Sphaerocavum sp. High frequency lake monitoring
		Former drinking water source Rain gauges Flow rate Cyanobacteria modelling: on-going Cyanobacteria dynamics modelling
	Urbanization + no sanitation Temperature, conductivity and Few data about the catchment
		infrastructure + no erosion control turbidity 35 infrastructure + no erosion control Automatic samplers: nutrients Automatic samplers: nutrients 30
	T (°C)	lake silting and storage capacity and dissolved oxygen 25 Lake Pampulha
		20	reduction Measurement buoy : 13/03 12/04 12/05 11/06 Management tool	11/07	10/08	09/09	09/10	08/11	08/12
	) Cyanobacteria (ug chl-a .L -1	Depth : 0.75 m Cyanobacteria blooms and macrophyte growth Touristic spot Fishing, sailing, landscape Temperature, conductivity, pH, O 2 and chl-a Meteorological forcing : INMET 20/08 09/09 29/09 19/10 08/11 28/11 18/12 0 100 200 300 400 22/02 13/03 02/04 22/04 12/05 01/06 21/06 0 50 100 150 200 250 300 High frequency watershed monitoring in progress 350 a b High frequency lake monitoring from September/2011 Cyanobacteria dynamics modelling preliminary results Watershed modelling and scenarios	10 12

*MAPLU project: Funded by FINEP (Brazilian Research Agency)

http://www.skyscrapercity.comi


Jovica Stanojevic

Valentina Parigi

Erwan Bimbard

Rosa Tualle-Brouri

Alexei Ourjoumtsev

Philippe Grangier

Controlling the quantum state of a single photon emitted from a single polariton

We investigate in detail the optimal conditions for a high fidelity transfer from a single polariton state to a single photon state and subsequent homodyne detection of the single photon. We assume that, using various possible techniques, the single polariton has initially been stored as a spin-wave grating in a cloud of cold atoms inside a low-finesse cavity. This state is then transferred to a single photon optical pulse using an auxiliary beam. We optimize the retrieval efficiency and determine the mode of the local oscillator that maximizes the homodyne efficiency of such a photon. We find that both efficiencies can have values close to one in a large region of experimental parameters.

I. INTRODUCTION

The generation and characterization of well-controlled single photon states have accomplished considerable progress during recent years. Beyond the usual nonclassical effects obtained using single photon counting and intensity correlation measurements [START_REF]Quantum Theory of Optical Coherence, Selected Papers and Lectures[END_REF], interference effects offer very efficient ways to get more information about the quality of single-photon wave packets. For instance, observing coalescence, that is the famous "Hong-Ou-Mandel" effect [START_REF] Hong | [END_REF], between two single photons emitted by two different atomic sources is a very good way to ensure that these two photons are indistinguishable [3].

Another avenue in order to fully characterize a single photon state is using homodyne detection, and applying quantum tomography techniques in order to reconstruct the Wigner function W (x, p) of the photon in phase space, where x and p are the quadrature operators of the quantized electric field [4]. Such a method provides a full characterization of the quantum state of the single photon, which is very intuitive because the quality of the single photon state, and especially its purity, directly translate into the negativity of the Wigner function.

In the optical domain quantum homodyne tomography has been implemented for various number states, including one- [4] and two- [5,6] photon Fock states. In these experiments, the usual method to produce a singlephoton state is to use parametric fluorescence emitting non-degenerate pairs of photons, and to detect one of them using a photon counter, such as an avalanche photodiode. Then the other beam is projected in the desired one-photon Fock state, which is analyzed by homodyne detection and then quantum tomography. It has to be noted that the quantum efficiency of the homodyne detection channel has to be very high: any loss at that stage degrades the purity of the single photon state by mixing it with vacuum, and quickly destroy the negativity of the Wigner function. This method works quite well, but it is intrinsically non-deterministic : the probability of the first (photon counting) event must be low, and then the photon cannot be emitted "on demand", when needed for applications.

Here we would like to investigate another scheme, FIG. 1: (Color online) An ensemble of N three-level atoms placed in a cavity is initially prepared in a single polariton state, with 1 atom in the metastable state |s and (N -1) atoms in the ground state |g . The atoms are coupled to the quantum field with collectively enhanced coupling w = g √ N , where g is the single photon Rabi frequency. By applying a classical control field (red) with Rabi frequency Ω(t) a singlephoton state is generated with well defined properties. The cavity detuning δ is used to optimize the photon retrieval efficiency.

which can, at least in principle, be made more deterministic. The idea is to emit the single photon from a so-called "polariton" state, where a single long-lived spin excitation is distributed over many atoms. It can be converted into a free-propagating single photon using a control laser beam (see Fig. 1). As we will show below, given the presence of the polariton the source is deterministic: the phase-matching condition between the readout laser beam, the atomic spin wave and the single photon leads to a collective enhancement effect which ensures that the photon is emitted in a well-defined spatial and temporal mode. For the whole scheme to be deterministic, the polariton should also be prepared in a deterministic way, which can be done for instance by using the so-called "Rydberg blockade" mechanism as explained in recent theoretical papers [7,8]. However, even when the preparation of the polariton is probabilistic, one may exploit the fact that the single photon can be emitted on demand, just when it is actually needed, within the relatively long coherence time of the spin wave. In this article we will not consider in more details how the polariton is prepared [9][10][11][12][13], but rather look at the following question: assuming that a single polariton is imprinted on an atomic cloud within an optical cavity, how efficiently can we turn it back into a single photon within a single mode? What will be the amplitude and frequency of the emitted single photon wave-packet, and is it possible to control it, and to mode-match it efficiently onto a local oscillator? In order to answer these questions, we will use as a starting point the theoretical analysis presented in [9], which has been shown to be in good agreement with experiments [14]. We will extend this analysis in order to characterize very carefully the homodyne efficiency which can be obtained for such a single-photon pulse, since this efficiency is crucial for quantum homodyne tomography. Finally we will discuss various experimental and theoretical considerations.

II. THEORY

We assume that the atoms are confined in an optical cavity and a single polariton state has been prepared so we only analyze the retrieval of a single photon. We are primarily interested in situations in which the collectively enhanced coupling between the atoms and a quantized cavity radiation mode is always one of the dominant energy scales, but not necessarily the only one. This is an interesting regime for experiments and it has not been studied in detail before. The underlying atom-cavity model has been already described [9] and thus we only briefly introduce it.

A. Cavity model

Various techniques and strategies for photon storage and retrieval have been analyzed using a relatively simple model [9]. It is assumed that a cloud of N three-level atoms (see Fig. 1) with two metastable states |g and |s , and an excited state |e are coupled to a classical and quantized single-mode fields. The classical (read) field couples the states |s and |e while the quantized field couples the states |g and |e . The coupling to the classical field is characterized by a slowly varying Rabi frequency Ω(t) and the collectively enhanced coupling w between the atoms and the quantized cavity radiation mode is w = g √ N , where g is the single photon Rabi frequency.

The assumption that, initially, a single excitation is shared by a very large number of atoms practically means that almost all atoms are in the ground state at all times so that the following approximate relations are justified: Σ N i=1 σ i gg ≈ N , and Σ N i=1 σ i ee ≈ Σ N i=1 σ i ss ≈ Σ N i=1 σ i es ≈ 0, where σ i uv = |u i i v|. Under these simplifications and in the rotating-wave approximation, we define the atomic polarization annihilation operator P (t) = N i=1 σ i ge / √ N , and spin-wave annihilation operator Ŝ(t) = N i=1 σ i gs / √ N . In the general case, these operators also include phase factors exp(-i k.r i ), where k is the laser wave-vector, and r i the position of atom number i. There phase factors can be absorbed by redefining the state |g i . In Appendix A we write the Langevin equations for the coupled propagation of the field annihilation operator Ê(t) and atomic operators P (t) and Ŝ(t) in the Heisenberg picture, and we show that it is actually enough to consider the c-number equations :

Ė = -(κ + iδ) E + iwP , (1) Ṗ
= -(γ + i∆) P + iwE + iΩS, (2)
Ṡ = iΩ * P, (3)
where κ and γ are respectively the cavity and polarization decay rates, while δ is the detuning between the laser-driven Raman light and the cavity, whereas ∆ the detuning between the driving laser and the atomic line (see Fig. 1).

It should be pointed out that the c-numbers E(t), P (t) and S(t) are not averages of the corresponding operators, but can be seen as quantum amplitudes, which allow one to evaluate normally ordered products of operators (see details in Appendix A). An intuitive justification for this approach is that, since there is only one excitation in the system, there will be only one emitted photon in the cavity mode; we are therefore calculating the amplitude of this single photon wave packet. The initial conditions are E(t 0 ) = P (t 0 ) = 0 and S(t 0 ) = 1, where t 0 indicates the beginning of the retrieval sequence.

The retrieval efficiency η is defined as the ratio of retrieved photons and stored excitations. Since the number of stored excitations is assumed to be one, η is equal to the number of retrieved photons

η = ∞ t0 dt|E out (t)| 2 , (4)
where

E out (t) = √ 2κ E(t). (5)
As said before, we do not need to use quantum operators in Eq. ( 4) since we never have more than one photon in the quantized radiation mode.

The propagation equations for (E R , P I , S R ) and (E I , P R , S I ) decouple for ∆ = δ = 0 and real Ω, where a R and a I are the real and imaginary part of a complex variable a. Consequently, E and S are real and P is imaginary in that case. For real Ω, there is the following mapping between the solutions for ±∆ E(-∆, -δ, t) = E * (∆, δ, t), (6) P (-∆, -δ, t) = -P * (∆, δ, t), (7) S(-∆, -δ, t) = S * (∆, δ, t).

(

) 8
For non-zero ∆, the interaction between the atoms and the quantized radiation mode shifts the cavity resonance.

Therefore, to compensate this shift and maximize the efficiency (4), we introduce some cavity detuning δ = 0. Relations ( 6)- (8) impose the mapping between the optimized cavity detunings δ opt as follows

δ opt (-∆) = -δ opt (∆). (9)
As explained in Section III, other mappings relevant for the optimization of homodyne detection (HD) also follow from Eqs. ( 6)- (8). All these mappings also have practical usage since they reduce numerical work by half.

B. Retrieval efficiency

We want to explore the parameter space by means of numerical calculations to find regions of high photon retrieval efficiency. In many experiments the emission of the single photon is detected by a counter, such as an avalanche photodiode, and the relevant parameter is η as defined by Eqs. 4 and 5. Here, we are also interested in using homodyne detection, which allows us to analyze fully the single photon wavepacket, and e.g. to reconstruct its Wigner function by using quantum optical tomography. Therefore we will also look for optimal conditions for homodyne detection.

Let us first review the expression for the maximal photon retrieval efficiency in a general case since our search for optimal retrieval conditions do not necessarily involve any special cases such as the bad-cavity limit (κ ≫ w). Actually, based on experimental considerations, the typical range of parameters we are interested in corresponds to w >> κ ∼ γ. We will show now that the efficiency maximum is then the same one derived in [9], though we are no more in the bad cavity limit.

A common way to simplify the problem is to try to find the first integrals of the equations of motion. Since our system is dissipative (κ, γ = 0), we may at least look for some relations analogous to the first integrals. The following relation gives just that

d dt |E| 2 + |P | 2 + |S| 2 = -2γ|P | 2 -2κ|E| 2 . ( 10
)
We are interested in the conditions for high retrieval efficiencies η. Using Eqs. ( 5) and (10), as well as the conditions

S(t 0 ) = 1, E(t 0 ) = E(∞) = P (t 0 ) = P (∞) = S(∞) = 0, one gets the following relation 1 = η + 2γ ∞ t0 |P (t)| 2 dt. (11)
The boundary condition S(∞) = 0 is only fulfilled for maximal efficiencies since the spin-wave decay rate is neglected. It is convenient to rewrite the integral in the last equation in the frequency space and then use the relation between the Fourier components P (ω) and E(ω) obtained from Eq. ( 1)

P (ω) = δ -ω -iκ w E(ω). (12)
Therefore

∞ t0 |P (t)| 2 dt = ∞ ∞ |P (ω)| 2 dω = κ 2 w 2 ∞ -∞ 1+ (ω -δ) 2 κ 2 |E(ω)| 2 dω, (13) which implies 2γ ∞ t0 |P (t)| 2 dt κγ w 2   2κ ∞ t0 |E(t)| 2 dt   = 1 C η, (14)
where C = w 2 /κγ is the cooperativity factor. Combining this condition and Eq. ( 11), we get the efficiency limit

η C 1 + C . ( 15
)
As shown in [9] for δ = 0, it is possible to achieve the highest efficiency C/(1 + C) in some particular regimes. Since the relation ( 15) is very general, introducing an optimized cavity detuning δ opt = 0 may allow us to optimize the efficiency, especially for low intensities of the driving field, but it will never exceed this limit. A necessary, but not sufficient, condition to have a maximal efficiency is that the term (ω -δ) 2 /κ 2 in Eq. ( 13) should be negligible. For δ = 0, this happens if the range of frequencies ω corresponding to non-vanishing P (ω) and E(ω) is much smaller than κ. In this case, achieving high efficiencies would require long read pulses or high Ω. However, according to Eq. ( 13), the efficiency can be improved significantly by introducing a cavity detuning δ that cancels ω corresponding to the maximum of |E(ω)|. A consequence is that δ opt is comparable with the spectral width of |P (ω)| and |E(ω)| so that we expect that both ω 2 /κ 2 and δ 2 opt /κ 2 are negligible for very long pulses. Because the term ω 2 /κ 2 in Eq. ( 13) originates from Ė(t) in Eq. ( 1), we may think that Ė(t) can be neglected for very high efficiencies. Dropping this term yields to the 'bad cavity' relation between E(t) and P (t)

-κE + iwP = 0, (16)
and therefore (κ, w > 0)

κ|E| = w|P |. (17)
In general, as shown in detail in Appendix B, for sufficiently long pulses and high efficiencies a better approximation can be obtained by noting that Eq. ( 1) is the first-order expansion of

-κE(t + δt) + iwP (t) = 0, ( 18
)
with respect to the short constant time delay δt = 1/κ. Equation ( 18) substitutes to Eq. ( 16) in the time domain, while in the frequency domain it becomes

-e iφω κE(ω) + iwP (ω) = 0 ( 19
)
with φ ω = -ωδt, and the relationship (17) remains valid.

C. Solution for large w

In this section, we assume that w is at least large but Ω and ∆ can be large as well. For w ≫ κ, δ, relation (12) shows that the ratio |P |/|E| is small as long as the pulsewidths of Ω, E, and P are of the same order of magnitude

|P | |E| ≪ 1, (20)
Because S ∼ 1 at short times (due to the initial condition S(t 0 ) = 1), Eqs. ( 20) and (2) lead to the following estimate of E for small ∆,

E ∼ - Ω(t) w . (21)
Let us introduce a new variable Q = E + pP and take p for which the following equation is held

Q = λQ + ipΩS. ( 22
)
There are two combinations of p, λ for which the last equation is fulfilled

λ ± = -1 2 γ +i∆+k ± i 4w 2 -(γ +i∆ -k) 2 , ( 23
)
p ± = i 2w γ +i∆-k ± i 4w 2 -(γ +i∆-k) 2 , (24)
where k = κ + iδ. We choose p, λ which facilitate the adiabatic elimination of Q in Eq. ( 22)

E + pP = - ip λ ΩS ≡ - w β 2 ΩS, (25)
where

β 2 = -iw λ p (26)
is introduced for convenience. Estimates (20) and (21) show that |P | is a much smaller quantity compared to both |E| and |S|. Therefore, dropping the small term pP in Eq. (25) has little effect to the relationship between much larger E and

S E = - w β 2 ΩS. ( 27
)
The last equation can be obtained in a different way and somewhat different approximations. For some parameters, we can use both λ ± and p ± to eliminate P in Eq. ( 25) and get relation ( 27) again with a new definition for β

β 2 = -iw λ + λ - λ + -λ - p + -p - p + p - . (28)
For our typical range of parameters, there is no practical advantage of using either Eq. (26) or Eq. ( 28) but it turns out that the definition (28) is more suitable for very

FIG. 2: (Color online) Time delay between P (t) and E (t). (a)

For high retrieval efficiencies and small Ω/Ω (i.e long pulses), the cavity field E is delayed with respect to P by a constant delay δt = 1/κ: E (t) ≈ iwP (t-δt)/κ. The numerical solution for P (solid red) is compared with P bd (t) ≡ -iκE (t)/w (solid green) and P shift (t) = P bd (t + δt) (dashed blue). (b) For shorter pulses (τ = 150 ns) this simple delay approximation breaks and P (t) and E (t) acquire different shapes.

long pulsewidths, as shown on Fig. 2. Since S ≈ 1 in the beginning of a pulse sequence, we see that, for initial times, the shape of |E| and |Ω| are always the same and |E| ≪ |Ω|. Taking the derivative of Eq. ( 27) with respect to t and using Eqs. ( 1) and (3) to eliminate Ṡ and P , we get

1 + |Ω| 2 β 2 Ė = Ω Ω -k |Ω| 2 β 2 E. ( 29
)
For real Ω, the solution of this equation consistent with the initial condition E(t 0 ) = 0 is relatively simple

E(t) = - w β Ω β 2 + Ω 2 exp -k t t0 Ω 2 β 2 + Ω 2 dt ′ . (30)
where the multiplicative constant (-w/β) has been deby taking limits of Eqs. (30) and (27) at t → t 0 and utilizing the initial condition S(t 0 ) = 1. Similarly, the solution for S is

S(t) = β β 2 + Ω 2 exp -k t t0 Ω 2 β 2 + Ω 2 dt ′ . ( 31
)
Finally, the solution for P in our approximation is

P (t) = iβ(kΩ + Ω) (β 2 + Ω 2 ) 3/2 exp -k t t0 Ω 2 β 2 + Ω 2 dt ′ . (32)
Let us notice that if |Ω| 2 ≪ |β| 2 , then E(t) ≈ -wΩ(t)/β 2 = -ipΩ(t)/λ, which means that Ω and E have the same pulse shape. However, according to Eq. (27) and (31), this also means that S ≈ 1 and hence the retrieval efficiency is very low. Solutions (30)-(32) are excellent approximations for large w and ∆ ≈ 0, they work rather well for moderate to large ∆, but for large detuning they become inaccurate, as illustrated in Fig. 3. The reason why the approximation worsens depends on which definition of β is used. For ∆ = δ = 0, both |λ ± | are the same but they can be very different for larger ∆. The same statement holds for |p ± |. Since the larger |p| corresponds to the larger |λ|, dropping the term pP in Eq. (25) becomes problematic for larger ∆ even though the adiabatic elimination itself works better. On the other hand, the problem with using Eq. (28) for larger ∆ is that the adiabatic elimination corresponding to the smaller |λ| becomes unreliable.

III. OPTIMIZATION OF HOMODYNE DETECTION

Homodyning consists in interfering two light fields: the local oscillator (LO), and the measured "quantum" beam (or signal beam), associated with mode annihilation operator â(t). The local oscillator is assumed to be monochromatic and to be in a coherent state |α , with α = |α|e -iθ where |α| 2 is the LO photon flux and θ is a controllable LO phase. More precisely, the signal from the homodyne detection (HD) is obtained by taking the difference of the two output photocurrents, resulting from the interference of the signal and LO. This difference pho-tocurrent i(t) is proportional to the following operator i(t) ∝ α * â(t) + α â † (t) /|α| = e iθ â(t)+e -iθ â † (t) (33) which is a quadrature operator with respect to the signal field, determined by the choice of θ.

In the present work, we will consider a situation with a continuous single-mode LO, but where the signal beam corresponds to a single photon, which is an essentially transient field, triggered by the "read-out" pulse sent onto the atoms in the cavity. This situation is somehow similar to the homodyne detection of transient homodyne signals following a photodetection event, as measured for instance in [16,17]. The idea is then to "extract" from the stationnary noise coming out from the homodyne detection a non-stationnary signal, through an appropriate time-domain filtering, which has to be matched as well as possible to the expected "single-photon wavepacket".

Two effects can essentially degrade this matching : the first one is the efficiency η introduced in the previous section, i.e. the probability that the photon is effectively emitted in the output mode â(t), and the second one is the quality of the interference between the photon and the LO, represented by an overlap factor χ defined below.

In order to examine the effect of the quality of the interference, we will neglect all technical mode-matching issues, and consider only the effect of the interference of the LO with the single-photon amplitude E(t) introduced in the previous section. Our goal is to optimally match the mode E det (t) measured by the homodyne detector to the mode of the single photon defined by E(t). The amplitude f (t) = |E det (t)| can be controlled either by modulating the intensity of the local oscillator or by rapidly sampling the signal and applying a numerical filter f (t) to the data, and in the following we will assume that the amplitudes of the two modes are identical:

f (t) = |E out (t)| √ η . (34)
The phase θ(t) of E det (t) includes a constant phase shift ϕ 0 , which can be introduced on the local oscillator using a piezoelectric actuator, and correspond to various choices for the measured quadrature. In the present situation of a single photon state, we don't expect this phase to have any effect on the measured signal. In general a rotating phase ω 0 t can be added by changing the LO frequency by ω 0 . More complex functions θ(t) can be realized in principle, for example by chirping the local oscillator frequency, but for simplicity we will assume that θ(t) = ϕ 0 + ω 0 t.

The mode-matching efficiency can be optimized by maximizing the variance of the homodyne detection output, integrated over the output pulse. Up to a multiplica-tive constant this variance is given by the expression

V = dtf (t)i(t) 2 (35) = dtdt ′ f (t)f (t ′ ) i(t)i(t ′ ) = dtdt ′ f (t)f (t ′ ) × ( â(t)â † (t ′ ) e iθ(t)-iθ(t ′ ) + â † (t)â(t ′ ) e iθ(t ′ )-iθ(t) + â(t)â(t ′ ) e iθ(t)+iθ(t ′ ) + â † (t)â † (t ′ ) e -iθ(t)-iθ(t ′ ) )
In the present case, the initial conditions are such that the output cavity field will be a diagonal mixture of zero and one photon, and therefore the average values â(t)â(t ′ ) and â † (t)â † (t ′ ) are zero for all times (see discussion in Appendix A). Since the calculated quantities provide normally ordered averages, one uses the commutation relations to get finally :

V = 1 + 2 dtdt ′ f (t)f (t ′ ) â † (t)â(t ′ ) e iθ(t ′ )-iθ(t) (36) = 1 + 2 I 0 where I 0 = dtf (t)E out (t)e -iω0t 2 . ( 37
)
As expected V = 1 (vacuum noise) if I 0 = 0, and V = 3 (single photon state) if I 0 = 1. It should be noticed that the phase ϕ 0 does not appear any more (as it could be expected), and one has thus only to optimize V with respect to ω 0 . The optimal ω 0 , denoted by ω opt , is obtained by maximizing the integral I 0 and is found numerically. Up to a constant prefactor, the optimized detector signal I max is equal to the maximal value of I 0 . Interpreting the integral in Eq. (37) as a scalar product of f (t)e iω0t and E out (t), and using the Cauchy-Schwarz inequality for real and normalized ( f = 1) filters, we get

I 0 ≤ E 2 = dt|E out (t)| 2 = η. ( 38
)
The upper bound I upper = η of I 0 , given by Eq. ( 38), can be used to define a measure χ of the optimization success, by the expression

χ = max(I 0 ) I upper = dtf (t)E out (t)e -iωoptt 2 η . ( 39
)
The mapping (6) implies the following mapping of the optimized frequency shift

ω opt (-∆) = -ω opt (∆), (40)
Nearly maximally allowed value η of I 0 is achieved by using the filtering function f (t), and by looking for the best possible linear (in time) approximation of the phase of E out (t) (see Fig. 4). As shown below, the regions of high values of η and χ have a very good overlap, which is important information for future experiments. Let us emphasize however that in practice, E out (t) in Eq. (39) has to be corrected by including experimental losses, which will result in lower overall detection efficiencies.

IV. RESULTS

In this section we present numerical results on optimal conditions for efficient retrieval and subsequent detection of a single photon. The efficiencies η and χ are evaluated using the numerical solution for E(t), assuming a readout pulse with a Gaussian shape (as a function of time). Based on a typical experimental design using Rb atoms in a cavity [18], we will take κ/2π = 9 MHz and γ/2π = 3 MHz as fixed values, and explore a parameter space consisting of (∆, Ω) points. Similarly, we will assume that the cooperativity factor is in the range C ∼ 100 [18], using for instance the two hyperfine states |g = 5s 1/2 , F = 2, M F = 2 and |s = 5s 87 Rb. The readout pulse duration is chosen short enough (150 ns) so that the motional decoherence can be neglected for a cold atom cloud [14,15].

1/2 , F = 1, M F = 1 of

A. Optimization of the cavity detuning

As pointed out in Section II B, if ∆ = 0 the frequency of the cavity mode can be shifted due to the presence of atoms. To compensate for this shift, we introduce a cavity detuning δ and numerically find δ = δ opt that maximizes either the efficiency η or the product χη. The most dramatic improvements occur at lower Ω, as shown by comparing Figs. 5 and6, and the dependence δ opt (∆, Ω) is shown in Fig. 7 for C = 100. The optimum cavity detuning for large Ω is rather small, but it further flattens the high efficiency "plateau" visible on Figs. 5 and6.

Optimizing χη gives the upper homodyne detection limit, according to Eq. (39). Ideally, the retrieval and detection efficiencies should be very close to each other

max[χ(δ, ω 0 )η(δ)] ∼ = max[η(δ)]. (41)
Since the right hand-side of the last relation only depends The efficiency χη is optimized with respect to frequency and cavity detuning (same parameters as in Fig. 5).

on δ, we can first find δ opt that optimizes η(δ) and then subsequently optimize χ(δ opt , ω 0 )η(δ opt ) with respect to ω 0 . This optimization procedure can be done very efficiently since only one parameter is varied at the time. This δ opt corresponds to the maximum retrieval efficiency and thus it is an useful parameter in itself. The values δ opt and ω opt obtained in this optimization step can then be used as an initial guess to find optimal δ and ω 0 for χη. However this second optimization step does not lead to any practical increase in χη, as long as the linear phase approximation is valid (see Fig. 4). In that case, the optimization of χη with respect to δ and ω 0 can be efficiently replaced by optimizing η(δ) with respect to δ and subsequently χ(δ opt , ω 0 ) with respect to ω 0 .

For example, in the case presented in Fig.

B. Calculation of the efficiencies η and χη

Numerical results for η and χη are shown in Fig. 5 and 6 for τ = 150 ns and C = 100. On Fig. 5 one has δ = 0 (i.e. no δ-optimization is performed), whereas Fig. 6 corresponds to the optimized cavity detuning δ opt (∆, Ω), shown on Fig. 7. We observe that both quantities are symmetric under the transformation ∆ → -∆, consistently with the mappings ( 6) and (7).

Without cavity detuning optimization (see Fig. 5), a plateau of nearly maximum retrieval efficiencies is reached for sufficiently large Rabi frequencies Ω of the read field. Increasing the pulse duration τ helps establishing the efficiency upper limit C/(C + 1) in a larger region in the parameter space, because the spectral distribution |E(ω)| becomes more narrow for longer pulses so that ∼ ω 2 /κ 2 in Eq. ( 13) is negligible. Increasing the cooperativity factor increases the efficiency limit C/(C+1) but also widens the plateau of nearly maximal efficiencies. Since we keep κ and γ fixed, larger C values imply larger couplings w. However, for ∆ ≈ δ ≈ 0, the growth rate of η as Ω increases is lower for larger w because, according to the solution (30), E is a function of the ratio Ω/β and β ∼ w. Consequently, we need to operate with more intense read fields for larger C in order to retrieve single photons with near certainty. For example, increasing C from 100 to 200 for ∆ = δ = 0 requires ∼ 40 % higher Rabi frequencies to reach the plateau region.

C. Optimized values of the efficiencies η and χη

In Fig. 6, we show η(∆, Ω) assuming that δ is optimized. As expected, the cavity optimization yields to large improvements in efficiency for smaller Ω. As Figs. 5a and 6a indicate, there is a region of relatively small ∆ FIG. 8: (Color online) Frequency shifts ν = ωopt/2π corresponding to optimized cavity detuning (the curve is antisymmetric with respect to the line ∆ = 0, ν = 0). These values of ν optimize the overall homodyne efficiency χη in a linear phase approximation (see Fig. 6). The parameters are the same as in Fig. 7.

and Ω in which the δ-optimization can significantly boost η to reach almost 1, as shown for instance in Fig. 4 a,b.

An overall effect of the δ-optimization to χη can be seen by comparing Figs. 5b and6b. As a general rule, the curve on Figs. 6b broadens and flattens, and high efficiency values become available for lower readout Rabi frequencies Ω. In general, increasing any of C and τ expands the region of maximal χη. Finally, Fig. 8 shows the behaviour of ν opt = ω opt (∆, Ω)/2π. The frequencies ν opt are rather small and of the order of a few MHz. For ∆ = 0 and real Ω, E(t) is real and negative and thus no phase correction is necessary. This could be a preferable experimental choice, provided that a high enough Ω is available. The function ω opt changes sign across the line ∆ = 0 due to the symmetry (40). As previously pointed out, the function ν opt exhibits distinguishable features for relatively small ∆ and Ω due to the cavity detuning optimization.

V. CONCLUSION

As a conclusion, we have calculated the amplitude and phase of a single-photon wave-packet, obtained by deterministically "reading out" a single polariton, which is assumed to be written beforehand onto a cloud of atoms within a low-finesse optical cavity. This "writing" can be done either non-deterministically by scattering a photon from the atoms [14], or deterministically by using Rydberg blockade techniques [7]. The results indicate that both the emission efficiency and the homodyne efficiency towards an adequate local oscillator can reach values very close to one. This opens the way towards the full quantum homodyne tomography of a single photon deterministically generated within a single mode.

As shown in ref. [9], the coupled propagation of the field annihilation operator Ê(t) and atomic operators P (t) and Ŝ(t) in the Heisenberg picture writes :

Ė = -(κ + iδ) Ê + iw P + √ 2κ Êin , (42)
Ṗ = -(γ + i∆) P + iw Ê + iΩ Ŝ + 2γ FP , (43) Ṡ = iΩ * P , (44)
where κ and γ are respectively the cavity and polarization decay rates, while δ is the detuning between the laser-driven Raman light and the cavity, and ∆ the detuning between the driving laser and the atomic line (see Fig. 1). Here we omit the spin-wave decay rate, but it can be reintroduced easily as shown in [9]. The (bosonic vacuum) Langevin noise terms are Êin and FP .

In the article, we treated Ê, P or Ŝ as complex numbers and ignored all noise terms. Here we want to explain what c-number quantity E(t), P (t) and S(t) represent and in which situations they can(not) be used. First, let us note that E(t), P (t) and S(t) are not expectation values of the corresponding operators because

Ê(t) = P (t) = Ŝ(t) = 0 (45)
is fulfilled at all times. The last expression is true because there is only one excitation in the system, which may be in either of the bosonic operators Ê, P or Ŝ. Our main goal in this appendix is to show that E(t), P (t) and S(t) determine the expectation values of normally ordered operators as follows

 † i (t) Âj (t) = A i (t) * A j (t), (46)
where Âi,j can be any of the operators Ê, P or Ŝ. We point out that the noise terms vanish in the propagation equations for ∂ t  † i (t) Âj (t) because we assume that the input noise terms Êin and FP are in the vacuum state. Then it is easy to check that all normally ordered noise correlations are zero [9]. This means that, as far as the averages  † i

where B i (t) and A j (t) satisfy Eqs. ( 47) and (48) with the initial conditions B i (t 0 ) = B i0 and A i (t 0 ) = A i0 . We can directly check that the products B i (t)A j (t) always satisfy the equation of motion for Bi (t) Âj (t) . However, to be solutions for Bi (t) Âj (t) , the products B i (t)A j (t) have to also satisfy the initial conditions. This is only possible if all Bi (t 0 ) Âj (t 0 ) are factorisable as required by Eq. ( 49). We emphasize that this factorization of Bi (t 0 ) Âj (t 0 ) is not a trivial condition for operators and it is not generally fulfilled.

In our case, the only initial nonzero expectation value is Ŝ † (t 0 ) Ŝ(t 0 ) so the initial averages are clearly factorisable. This means that the condition (49) is fulfilled and so is the relation (46). Therefore, we have shown that the expectation values of the normally ordered operators are exactly the same as for the products of the c-number quantities E(t), P (t) and S(t). For simplicity, S(t 0 ) = 1 is used because the only effect of another phase convention S(t 0 ) → S(t 0 )e iφ0 is that all E(t), P (t), and S(t) acquire the same additional phase factor e iφ0 which leaves the averages (46) unchanged.

In case where there is only one excitation in the system, the condition (49) is equivalent to the requirement that the initial state is a pure state. To show this, we consider here the evolution of the density matrix of the system (in the Schrödinger picture), rather than the evolution of the operators (in the Heisenberg picture) used up to now. We take as a basis the Fock states {|000 , |100 , |010 , |001 } for the bosonic operators Ê, P or Ŝ, where the states are labeled as |n E n P n S . The state |000 does not evolve under the Hamiltonian H of the considered three-level atoms (H |000 = 0), but only due to damping, and the Hamiltonian part of the evolution takes place within the subspace spanned by |100 , |010 , and |001 . Realistic atoms may require considering more states, but in our problem the calculation of this 3 × 3 sub-matrix is equivalent to the calculation of E(t), P (t) and S(t) as carried out in this paper. Denoting { Âi } = { Ê, P , Ŝ}, the density matrix elements are simply

r ij = Â † i Âj (t). ( 51
)
If the factorization condition (49) is satisfied, then r ij (t 0 ) = A * i (t 0 )A j (t 0 ), which implies ρ(t 0 ) is a projector, and therefore that the initial density matrix describes a pure state. The other implication stating that if the initial state is pure then Eq. ( 49) is satisfied follows directly from the definition (51). Even more importantly, due to

FIG. 3 :

 3 FIG.3:(Color online) Comparison between numerical (solid red) and analytical (dashed blue) solutions for the emitted signal field |E (t)|. The Rabi frequency Ω(t) of the driving field is also shown (dashed green), scaled down to match the pulse area of the corresponding cavity field |E (t)|. Let us note that the (very small) emitted pulse lies entirely within the envelope of the (very large) driving pulse, but its peak value does occur before the driving field peak value. The parameters are Ω/2π = 80 MHz, C = 200, δ = 0 MHz, and pulsewidth (FWHM) τ = 200 ns for several detunings. The agreement is practically exact for ∆ = 0, but it gradually worsens as ∆ increases.

FIG. 4 :

 4 FIG.4: (Color online) Linearization of the phase θE of the quantum field E (t). To get the highest homodyne efficiency, the phases of the quantum field and local oscillator have to be matched (see Section III). Assuming the phase θLO of the local oscillator to be linear in time, the optimal frequency ωopt of the LO and the optimal cavity detuning δopt that maximize the homodyne efficiency χη are numerically found for the generic parameters τ = 150 ns and C = 100 used in section IV. In panels a) and b), θE (red solid) and θLO (blue dashed) are shown for ∆/2π = 120 MHz and Ω/2π = 13 MHz. No cavity detuning optimization (δ = 0) is performed in case a) while an optimized δopt/2π = 15.5 MHz is used in case b). The homodyne efficiency χ starts to decrease for a given ∆ and sufficiently large Ω because θE becomes more complicated, as shown in panel c) for Ω/2π = 80 MHz and ∆/2π = 300 MHz, with δopt/2π = -2.8 MHz. Panel d) is the same as c), except for the cooperativity that becomes C=200, improving the phase matching quality.

FIG. 5 :

 5 FIG. 5: (Color online) Retrieval efficiency η (panel a) and overall homodyne efficiency χη (panel b) for non-optimized cavity detuning (i.e. δ = 0) for C = 100 and τ = 150 ns. The frequency shift ω0 is optimized to get the best χη for any given ∆ and Ω.

FIG. 6 :

 6 FIG. 6: (Color online) Retrieval efficiency η (panel a) and overall homodyne efficiency χη (panel b) for optimized cavity detuning.The efficiency χη is optimized with respect to frequency and cavity detuning (same parameters as in Fig.5).

  (4)b (resp. (4)c), the χη-optimization gives δ opt /2π = 15.5 MHz (resp. -2.8 MHz) while the η-optimization gives δ opt /2π =15.7 MHz (resp. -5.2 MHz). But even in the apparently less favourable case of Fig. (4)c, the products χη are almost the same and nearly maximal, because the decrease in χ is compensated by an increase in η, so that χη = 0.930 (resp. χη = 0.937) as shown on Fig. 4.

FIG. 7 :

 7 FIG. 7: (Color online) Optimized cavity detuning for τ = 150 ns and cooperativity C = 100, which corresponds to w/2π= 46.5 MHz. The cavity detuning δopt that maximizes the efficiency χη is found for any given (∆, Ω).

  t) Âj (t) are concerned, the propagation equations for ∂ t  † i (t) Âj (t) are the same ones as for ∂ t (A * i (t)A j (t)), which can be obtained from Eqs. (42)-(44) with all noise terms set to zero. Using this conclusion, the relation (46) follows from the following simple, but nonetheless important mathematical statement: If the evolution of operators Âi and Bi , 1 ≤ i ≤ n is linear and if the expectation values of all products Bi Âj at some t = t 0 can be factorized as Bi (t 0 ) Âj (t 0 ) = B i0 A j0 (49) then Bi (t) Âj (t) are factorisable at all times Bi (t) Âj (t) = B i (t)A j (t),

	Ȧi =	a is (t) Âs (t), ,	(47)
	s
	Ḃi =

s b is (t) Bs (t),

(48)

Acknowledgments This work is supported by the ERC Advanced Grant "DELPHI" (Grant Agreement Number 246669) of the European Research Council.

Eq. ( 50), the factorization is preserved at all times which implies that the density matrix takes the form ρ(t) = r(t)|000 000| + |ψ(t) ψ(t)| ,

where |ψ(t) = (E(t), P (t), S(t)) T is a (non-normalized) pure state. Therefore (up to a normalizing factor) the initial state |001 evolves into the superposition E(t) |100 + P (t) |010 + S(t) |001 . We see that the component of the single-photon state is exactly E(t).

Appendix B

We want to confirm Eq. ( 18) which asserts that, for sufficiently large pulsewidths τ , the quantum field E(t) simply follows P (t)

with a constant time delay δt = 1/κ. This relation is very similar to the relation between E(t) and P (t) in the bad cavity limit κ ≫ w (which is not fulfilled in our case) except for the delay δt. This similarity is related to the relation (17) in the frequency domain. In the time domain, we will show that it is closely related to the fact that the single photon escapes the cavity in a very early stage of the retrieval sequence.

Let us assume that the read pulses are parameterized as Ω(t) = Ω 0 f (t/τ ). By using the dimensionless quantities t/τ , kτ , . . ., the number of parameters in Eqs. ( 1)-( 3) is effectively reduced by one, which makes the scaled quantities more natural. We can interpret Eq. ( 18) as the definition of a new function δt(t). If δt(t)/τ is small at all times, then the conclusion δt = 1/κ directly follows that from Eqs. (1) and (18). Therefore, we only need to show that δt(t)/τ is indeed small at all times for sufficiently long read pulses. This can be explicitly shown using the solution (30). The implication that 1/κ becomes a short time scale in this limit also indicates that the single photon has been emitted very early in terms of relative times u = t/τ .

For long pulses, the real part of the argument of the exponential function in Eq. (30) acquires very large negative values which indicates that the photon emission has been terminated. The time interval during which the photon has been emitted can be estimated from the condition Re κτ u u0

Therefore, increasing τ at fixed Ω 0 results in a photon emission at shorter relative times u = t/τ . In the limit of very long pulses, all terms ∼ Ω/Ω are negligible. Consequently, according to Eq. ( 29),

holds at all times during the photon emission. Combining Eqs. ( 18), (1), and (55), we get that

is satisfied for all t during the photon emission. There are two logical possibilities consistent with the last equation: either δt(t)/τ is small at all times or δt(t) relates distant instants with the same |E| and phase θ E of the quantum field. The latter possibility is very restrictive and would be essentially an additional condition to the equations of motion ( 1)-( 3) since such distant correlations cannot naturally come from strictly local propagation equations. In fact, all panels in Fig. 4 clearly show that the points with the same |E| have different phases so that the only consistent possibility derived from Eq. ( 56) is that δt(t)/τ is small for all t and consequently δt = 1/κ.