Rosa Tualle-Brouri

Alexei Ourjoumtsev

Aurelien Dantan

Philippe Grangier

Martijn Wubs

Anders S Sørensen

A multimode model for projective photon-counting measurements

Keywords: numbers: 03.67.-a, 42.50.Dv, 03.65.Wj

We present a general model to account for the multimode nature of the quantum electromagnetic field in projective photon-counting measurements. We focus on photon-subtraction experiments, where non-gaussian states are produced conditionally. These are useful states for continuous-variable quantum information processing. We present a general method called mode reduction that reduces the multimode model to an effective two-mode problem. We apply this method to a multimode model describing broadband parametric downconversion, thereby improving the analysis of existing experimental results. The main improvement is that spatial and frequency filters before the photon detector are taken into account explicitly. We find excellent agreement with previously published experimental results, using fewer free parameters than before, and discuss the implications of our analysis for the optimized production of states with negative Wigner functions.

I. INTRODUCTION

The ability to prepare and measure specific quantum states of the light is the keystone of many quantum information processing (QIP) protocols. These states can be described either with discrete variables in terms of photons, or with continuous variables in terms of waves. In the latter case, the physical quantities of interest are the amplitude and the phase of the light wave, or their Cartesian counterparts called quadratures x and p. A very convenient representation of the quantum state is then provided by the Wigner function W (x, p), which corresponds to a quasi-probability distribution of the quadratures, 'quasi-' because W may assume negative values.

An important task for QIP is the ability to undo effects of decoherence by 'distillation': to obtain a single quantum state that is more pure from two or more copies that have undergone decoherence. Since states of light with gaussian Wigner functions cannot be distilled with gaussian operations [1,2], one is left with two strategies: either to distill gaussian states with non-gaussian operations, or to distill non-gaussian states with gaussian operations [3]. This paper is a contribution to the latter strategy, and focuses on the preparation of the nongaussian states rather than on their distillation.

The negativity of the Wigner function is a standard figure of merit, quantifying at the same time how nongaussian and how non-classical a quantum state is [4,5]. One way of obtaining non-gausian states is by conditional photon-counting measurements, as first proposed by Dakna et al. [6]. It was soon realized that such condi- * E-mail: rosa.tualle-brouri@institutoptique.fr tional measurements can improve quantum teleportation of continuous variables [7]. In recent years, several experiments [8,9,10,11,12,13,14,15] combining continuousand discrete-variable tools allowed for preparing and observing quantum states of free-propagating light with negative Wigner functions [16].

Many of these experiments are based on the use of a squeezed vacuum produced by parametric fluorescence, which involves many optical modes [17]. This multimode nature is exhibited in both continuous-wave (CW) operation, using optical parametric oscillators (OPO) below threshold [18,19], and in pulsed experiments with a single-pass high amplification. In order to make accurate predictions, it is therefore crucial to develop multimode theoretical models. This was done in [18,19,20,21] for setups using an OPO, and in [22,23] for pulsed 1-photon Fock state tomography in a case of very low squeezing, when Fock states expansion are limited to 1 photon only. However, these models do not fully account for all phenomena linked to the non-constant space and time profiles of the modes under study, the spatial pulse profile in the transverse direction for example, and they especially do not account for gain-induced distortions in the parametric amplification process [24]. As we will see later on, these phenomena are one main signature of this multimode nature, and are critical in the case of single-pass pulsed experiments.

In this paper we propose an alternative general framework to describe the generation of squeezed light, with a twofold goal: first, to introduce a method that reduces a multimode model to an effective two-mode description. Second, to show that a specific spatio-temporal multimode model for photon-subtraction experiments and our mode-reduction analysis thereof give an improved understanding of state-of-the-art photon-subtraction exper-

FIG. 1:

Projective photon-counting measurement: in our model we assume that we start with a multimode vacuum input state which is subject to a general multimode Bogoliubov transform U . After the transformation a few modes are filtered out and result in a detection event in the avalanche photodiode (APD). Such a detection event in the APD prepares a state in the single mode which is mode analyzed by homodyne detection.

iments.

In Sec. II we show how to reduce a multi-mode model to an effective two-mode model. This mode-reduction procedure is then applied in Sec. III to give an improved analysis of the photon-subtraction experiments of Ref. [10]. We discuss the method and its application and conclude in Sec. IV. Some technicalities are deferred to two Appendices.

II. REDUCTION OF MULTIMODE MODEL A. General multimode model

As a starting point we have a complete set of spatial, temporal, or spatio-temporal optical modes, in terms of which the light propagation can be described. The modes have field operators âj and â † j with bosonic commutation relations [â j , â † k ] = δ jk . The âj may also stand for continuous operators like â(t), where t is time, with commutation relation [â(t), â † (t ′ )] = δ(t -t ′ ), in which case sums over modes are to be replaced with integrals. To simplify the notation even further, we introduce â, the column vector of all the âj . Similarly â † is the row vector which has operators â † j as components. We shall also need the row vector âT and the column vector â * = â †T . Similar notation will be used for other vectors.

In this section we consider a general unitary transformation U in which output quadratures linearly depend on input quadratures, thus preserving the gaussian nature of the quantum fields. In fact, the only non-gaussian operation will be the projective measurement, corresponding to a detection event in a subset of the output modes by an avalanche photodiode (APD), see Figure 1. We will assume the evolution operator U in Fig. 1 to be a Bogoliubov transformation, where the output field operators depend linearly on the input ones:

âout = U † â U = u â + v â * , (1)
where u and v are two matrices which satisfy

uu † -vv † = u † u -v † v = 1, (2a) uv T -vu T = 0, (2b)
whereby the commutation relations of the field operators are preserved. The output state of the light can be characterized by doing homodyne measurements on a normalized mode described by ψ h , which has the mode operator

âh = ψ † h â. ( 3
)
This mode can be defined as the mode that perfectly matches the local oscillator of the homodyne detector.

After the Bogoliubov transform (1) it is described by

âh,out = U † âh U = ψ † h (u â + v â * ). ( 4
)
Without conditioning upon photon detection, and since the gaussian nature of the initial vacuum state is preserved by the Bogoliubov transform, the homodyne measurements will show gaussian Wigner functions corresponding to squeezed vacuum states of light.

The point is now how to describe the output state conditional upon detection of a photon by the APD, which is a projective measurement. Analogous to the homodyne detection, we assume the photon detection mode j to be described by a normalized state φ d (j), with corresponding field operator

âd (j) = φ † d (j) â. ( 5
)
After the time evolution described by Eq. ( 1), the output at the photon detector becomes:

âd,out (j) = U † âd (j)U = φ † d (j)(u â + v â * ). (6)
If we only know that a photon has been detected but not in which detection mode, then we should average the conditional output state over all these detection modes, as is shown in more detail below.

At this point we have a multimode output state that in principle we can update given the detection of a photon in the detector j. In practice it is convenient to first simplify the multimode expressions ( 4) and (6).

B. Mode reduction

The first and central step in the mode-reduction procedure is to rewrite the homodyne mode Eq. ( 4) in the form in terms of new mode operators a ( †) 0,1 with standard commutation relations. The coefficients η, α, and β in Eq. ( 7) are found as follows. Besides having standard commutations, â0 and â1 in Eq. ( 7) should annihilate the vacuum state, so that â0 must contain all annihilation operators in Eq. ( 4):

âh,out = η â0 + α â † 0 + β â † 1 , (7)
η â0 = ψ † h u â, (8)
with η fixed up to a phase factor by [â 0 , â † 0 ] = 1. We choose η to be real-valued and positive so that

η = ψ † h uu † ψ h . (9)
Furthermore, from

[â 0 , â † 1 ] = 0 it follows that α = [â 0 , âh,out ] = 1 η ψ † h u v T ψ * h . (10)
A complex value for α can be removed with a redefinition of the phase of the homodyne mode, which just means that we can assume that the state is squeezed in the x or p direction. Finally, β can be found from the fact that âh,out and â † h,out have commutator one, as have â1 and â † 1 :

β = η 2 -α 2 -1. (11)
Here again we used the freedom to choose β real-valued and nonnegative. This completes the mathematics of the mode reduction of the multimode homodyning signal. The physical argument that two and only two modes should remain goes as follows. The squeezed vacuum after the Bogoliubov transform can only be a centered gaussian state, hence it is fully described by only the variances V x and V p . The squeezed vacuum output can therefore be modeled [START_REF] Gardiner | Quantum Noise[END_REF][START_REF] Wenger | [END_REF] by a perfect single-mode degenerate optical parametric amplifier (DOPA) with squeezing parameter r, followed by a perfect non-degenerate optical parametric amplifier (NDOPA) with squeezing parameter gr, as presented in Figure 2. This is a two-mode model, with a Hilbert space H 2 . The values of the parameters r and g can be deduced from the two independent coefficients in Eqs. (9)(10)(11), using: η = cosh(r) cosh(gr) (12a) α = sinh(r) cosh(gr) (12b) β = sinh(gr).

(12c)

A related squeezing parameter that we will also use in the following is s ≡ exp(-2r).

C. Conditioning upon photon detection

We now condition upon the measurement of a click in the photon detector (APD). We assume to be in the limit that the average number of photons per pulse entering the photon detector is much less than one. Then a single click in the detector corresponds to the detection of a single photon.

One can make a reduced-mode description of the photon detection operator âd,out (j) of Eq. ( 5), analogous to Eq. ( 7). The operator âd,out (j) can be expanded into a part acting on H 2 , plus a component âd⊥ (j) acting on the complementary space H ⊥ orthogonal to

H 2 : âd,out (j) = γ j â † 0 + δ j â † 1 + ǫ j â0 + κ j â1 + âd⊥ (j). ( 13
)
Note that âd⊥ (j) on the right-hand side contains all the terms acting on H ⊥ , including creation operators. The coefficients in (13) can again be found by taking commutators, for example:

γ j = [â 0 , âd,out (j)] = 1 η ψ † h u v T φ * d (j) (14a)
δ j = [â 1 , âd,out (j)] = 1 β * [ φ † d (j) v v † ψ h -α * γ j ]. (14b)
In general, a click recorded in the APD corresponds to the measurement of at least one photon. In the limit of low detection probability, the action of the detection is the subtraction of a single photon. Note that this assumption is practically always obeyed if the parameter j is in a continuum (like in the case of spectral filtering), since the probability to have two photons exactly in the same mode is then negligible. Henceforth we assume to be in this limit. In the Heisenberg picture, a photon detection then corresponds to the application of the operator âd,out (j) to the initial state (i.e. to the vacuum state), followed by the normalization of the result:

|Ψ(j) = P -1/2 j âd,out (j)|0 ( 15
)
where P j is the detection probability for the mode j:

P j = 0|â † d,out (j)â d,out (j)|0 = 0|â † d⊥ (j)â d⊥ (j)|0 + |γ j | 2 + |δ j | 2 (16) = -→ φ † d (j)vv † -→ φ d (j).
Below we will use that the vacuum expectation value of nd⊥ (j) = â † d⊥ (j)â d⊥ (j) can be expressed as

P j -|γ j | 2 - |δ j | 2 .
Before continuing, it can be instructive to recall the concision allowed by the Heisenberg picture. In a singlemode problem, a photon-subtracted squeezed state is equivalent to a squeezed single-photon state: this case corresponds to the simple Bogoliubov transform âout = cosh(r)â in + sinh(r)â † in , which directly gives a pure 1photon state (after normalization) when applied to the vacuum. As the states do not evolve in the Heisenberg picture, they all can be considered as 'input' states; but when measured using output quadratures, this 1-photon state will appear to be squeezed.

We are going to use the same approach in the multimode case. One can first note that the conditioned state |Ψ(j) in Eq. ( 15) is already a 1-photon state. This state, however, does not belong to H 2 only, so that measurements output are not so obvious to compute. In fact, we are solely interested in expectation values g(a h,out , a † h,out ) of operators describing the output that is measured in the homodyne detector. Such expectation values can be written as

g(a h,out , a † h,out ) = ψ j | g(â h,out , â † h,out ) |ψ j (17a)
= Tr{g(â h,out , â † h,out )|ψ j ψ j |}.(17b) In Eq. (17b), the trace can be separated into a trace over H 2 and a trace over H ⊥ , and the latter does not act on the function g(â h,out , â † h,out ), whose expectation value is then:

= Tr 2 {g(â h,out , â † h,out )Tr ⊥ |ψ j ψ j |} (18a) = Tr 2 {g(â h,out , â † h,out )ρ j } (18b)
All quantities of interest can therefore be deduced from the input state reduced density matrix ρ j , acting in H 2 , and the crucial advantage of mode reduction is to allow a simple expression for this matrix: Writing |0 = |00 ⊗ |0 ⊥ , where |00 and |0 ⊥ are the ground states of H 2 and H ⊥ , respectively, and using Eqs. (13,15), we directly obtain:

ρ j = Tr ⊥ |ψ j ψ j | (19)
= (1 -ξ j ) |00 00| + ξ j â † θ (j)|00 00|â θ (j), in terms of the modal purity

ξ j = |γ j | 2 + |δ j | 2 P j , (20)
and where â † θ (j) = cos θ j â † 0 + sin θ j â † 1 , with tan θ j = δ j /γ j (21) is an operator that creates a single photon in a superposition of mode 0 and mode 1. The state (19) produced from a detection event in mode j is a mixed state, mixing vacuum and a single-photon state with weight ξ j . Without conditioning or in the limit ξ j → 0, we have ρ j = |00 00|.

D. Wigner functions

Squeezed vacuum.-Before determining the output Wigner function corresponding to the conditional state (19), it is instructive to first determine the Wigner function of the output state in the simplest experimental situation, where we ignore the photon detector. The input state is then ρ = |00 00|. We will make use of the standard Wigner functions of the vacuum W 0 (x, p) = exp(-r 2 )/π and of single-photon states

W 1 (x, p) = (2r 2 -1) exp(-r 2 )/π, both with r 2 = x 2 + p 2 . Clearly, W [ρ](x 0 , p 0 ; x 1 , p 1 ) equals W 0 (x 0 , p 0 )W 0 (x 1 , p 1 ).
In order to obtain the output Wigner function for the homodyne mode, we wish to express W [ρ] as a function of x h,out , p h,out (defined from the output homodyne mode âh,out given by Eq. ( 7)). This requires the introduction of another mode â1,out orthogonal to a h,out , so that the transformation (â 0 , â1 ) → (â h,out , â1,out ) is symplectic (i.e. commutation relations are preserved). Using the model of Fig. 2, one can choose â1,out of the form:

â1,out = β α â0 + η â † 0 1 + β 2 + 1 + β 2 â1 . ( 22
)
This form is by no means unique, but this does not pose a problem since mode 1 out will eventually be integrated out. One can now invert the relations (7,22), thereby expressing x 0,1 , p 0,1 as a function of x h,out , p h,out , x 1,out , and p 1,out . After tracing over mode 1 out , which amounts to integrating over x 1,out and p 1,out , we find the output signal entering the homodyne detector to be a squeezed vacuum state with a gaussian Wigner function

W 0,sqz (x, p) = 1 π V x V p exp - x 2 V x - p 2 V p , (23)
where x and p stand for x h,out and p h,out and with variances

V x = (η + α) 2 + β 2 (24a) V p = (η -α) 2 + β 2 . ( 24b
)
We can now invert Eq. ( 24) and rewrite the three modereduction parameters α, η, and β in Eq. ( 7) in terms of the variances, giving:

η = V p + V x + 2 2 V x + V p + 2 , (25a)
α = V x -V p 2 V x + V p + 2 , ( 25b
)
β = 2 V x V p -1 2 V x + V p + 2 . ( 25c
)
In the following, we will keep writing α, η, and β to shorten notation. It should be kept in mind, however, that Eq. ( 25) directly expresses these parameters in terms of the measurable variances V x,p of the squeezed vacuum.

In particular, β vanishes for minimal-uncertainty states.

Notice also that the parametrization for the modereduction parameters [START_REF] Gardiner | Quantum Noise[END_REF] is equivalent to the one in Eq. ( 12) in terms of squeezing parameters r and g. Thus r and g can be expressed in terms of the variances V x,p , and vice versa.

Photon-subtracted squeezed vacuum.-As for the squeezed vacuum, we now calculate the Wigner function for the photon-subtracted squeezed vacuum, starting with the initial state (19). The mode ( 21) has a onephoton excitation in state (19). The orthogonal mode with creation operator â † θ+π/2 (j) is not excited. Hence the Wigner function corresponding to the state (19) is

W (x 0 , p 0 ; x 1 , p 1 ) = (1 -ξ j ) W 0 (x 0 , p 0 )W 0 (x 1 , p 1 ) (26) + ξ j W 1 (x θj , p θj )W 0 (x θj +π/2 , p θj+π/2 ).
Note that in this expression, quadratures x θj ,θj+π/2 , p θj,θj +π/2 can be easily expressed as functions of quadratures x 0,1 , p 0,1 using ( 21). As before, the Wigner function for the output signal is found by using the symplectic transformation defined by Eqs. ( 7) and (22). By tracing again over the mode 1 out , we obtain (see appendix A)

W j (x, p) = C j + 2A j x 2 V 2 x + 2B j p 2 V 2 p + D j xp V x V p W 0,sqz . (27)
The constants in this Wigner function are given by

A j = P -1 j |γ j (η + α) + δ j β| 2 , (28a)
B j = P -1 j |γ j (η -α) -δ j β| 2 , (28b)
C j = 1 -A j /V x -B j /V p , (28c)
D j = -8P -1 j Im(γ * j δ j )ηβ. (28d)
Note that the Wigner function [START_REF] Leonhardt | Measuring the Quantum state of Light[END_REF] of the photonsubtracted squeezed state differs from the Wigner function of the squeezed vacuum W 0,sqz (x, p) of Eq. ( 23) only because of the polynomial in x and p between the large brackets. The same quantities V x,p as in Eq. ( 24) show up, with or without conditioning. Since in general W (x, p) ≥ -π -1 [START_REF] Leonhardt | Measuring the Quantum state of Light[END_REF], we find the condition C j ≥ -1. Averaged Wigner functions.-Practical detectors do not resolve with infinite precision when and where photons are detected. We should therefore average over all possible microscopic states that agree with the detection record. We assumed in Sec. II C that the average number of photons detected per pulse in the APD is much smaller than one. Averaging over unresolved detection events is then equivalent to averaging over single-photon subtraction events.

The Wigner transformation of the density matrix is a linear transformation. Therefore, the averaged Wigner function W (x, p) is simply obtained by replacing A j ...D j in (27) by A...D, with the notation

X = P -1 tot j P j X j . (29)
Here P tot is the sum of the probabilities P j of microscopic states that agree with the detection record. From

Eq. ( 28) it follows that averaged quantities A...D involve sums like j |γ j | 2 , j |δ j | 2 or j γ * j δ j ; in the following we will write the respective averages as |γ| 2 , |δ| 2 or γ * δ.

The Wigner function [START_REF] Leonhardt | Measuring the Quantum state of Light[END_REF] of the photon-subtracted state and its detection-averaged version have a very general significance. Before going to Sec. III, devoted to a practical implementation of these results, let us finish with some reflections on the detection modes.

E. Detection modes

All the previous results were derived through the use of a set of detection modes âd (j). The coefficients entering in the averaged Wigner function W involve quantities like P tot , |γ| 2 , |δ| 2 or γ * δ, in which the detection operators only appear through their projection operator:

Π = j φ d (j) φ † d (j). ( 30
)
Hence one would find the same predicted averages if one would employ a different set of detection modes that has the same associated projector. This projector describes how the setup filters the signal before it enters the photon detector. For example, a time-domain filtering system will be described by Π = Π T , where Π T can be written in terms of a set of modes â(t) that are labeled by the time t:

Π T (t, t ′ ) = T (t)δ(t -t ′ ), (31)
where T (t) = 1 when the APD is switched on, and T (t) = 0 otherwise. To give another important example, a spectral slit can be described with a set of modes â(ω) with an associated projector

Π Ω (ω, ω ′ ) = T (ω)δ(ω -ω ′ ), (32)
where T (ω) = 0 if the frequency ω is filtered out, and T (ω) = 1 otherwise. Above we have assumed that the filtering of the signal after its production in the DOPA was included in the Bogoliubov transform U . Below we will give an alternative description, in which U is separated into the transformation due to the production of the squeezed light in the DOPA, and the subsequent filtering before detection. This alternative description will enable a more straightforward comparison with the empirical model in Sec. III.

So, instead of the input-output transform Eq. ( 6) for the photon detection operator, we now write âd,out (j

) = φ † d (j)u f (u â + v â * ), ( 33
)
where the Bogoliubov matrix u f accounts for filters (as the filters are passive, we have v f = 0); this matrix u f is of course unitary, even if the filters can present losses. In fact, losses will be modeled using beamsplitters, where the lost energy is reflected into auxiliary non-relevant modes. These modes do not interact with the rest of the experiment (i.e. they are unaffected by Bogoliubov transform u, v) and will not reach the APD. But all the other modes, referred to as relevant modes, should be considered as detection modes, and we then have:

j φ d (j) φ † d (j) = Π r , (34)
where Π r is the projector onto the subspace of relevant modes. Let Π r be the projector on the non-relevant modes. As the latter are unaffected by the transform u, v, we have Π r (u â + v â * ) = Π r â. Inserting the relation Π r + Π r = 1 1 into (33) then leads to âd,out (j

) = φ † d (j)u f Π r (u â + v â * ) + φ † d (j)u f Π r â. ( 35
)
Here the last term on the right annihilates vacuum, and commutes with annihilation operators like â0 or â1 : this term will add no contribution to the results of the previous subsection. The only change therefore consists in the substitution φ → Π r u † f φ, so that the operator Π in Eq. ( 30) should be replaced by

Π ′ = Π r u † f j φ d (j) φ † d (j)u f Π r = F † F, (36)
where F = Π r u f Π r and where we have used standard properties of projection operators (Π = Π † , Π 2 = Π). The operator F represents the action of the filters restricted to the subspace of relevant modes. If F is a projector, such as Π T of Eq. ( 31) or Π Ω of Eq. ( 32), then the effect of Π ′ is the same as of Π in Eq. ( 30). This can be easily understood: it is equivalent to say that the filtered modes are blocked, or that they are first redirected into auxiliary modes, and then blocked.

The operator Π ′ of Eq. ( 36) is a more general quantity than Π in Eq. (30), however, since Π ′ need not be a projection operator. It can for instance account for partial absorption of the modes. In that case, the spectral transmission T (ω) in Eq. (32) can assume any value between 0 and 1, to account for filtering systems more complex than a simple spectral slit.

Furthermore, the above expressions can simply be generalized to situations where several filters are used. For example, if a spectral slit Π Ω is followed by a time-domain filter Π T , then the above expressions for F and Π ′ become

F = Π T Π Ω , leading to Π ′ = Π Ω Π T Π Ω .

III. APPLICATION

At this stage, we have a complete description of the final state starting from the Bogoliubov transform (1). The results of the previous section are generally valid, since we started with a multimode model that was left unspecified. Our purpose now is to establish a concrete link with the photon-subtraction experiment as described in Ref. [10], and to improve its analysis.

A. Photon subtraction experiment

In the experiment by Ourjoumtsev et al. [10], pulses of squeezed light are produced. The setup is sketched in Figure 3. A squeezed vacuum, produced in a singlepass DOPA (a KNbO 3 crystal) by down-conversion of frequency-doubled femtosecond laser pulses, is sampled by a beam splitter with low reflectivity R = 1 -T . Two mode filters are placed in front of the APD: a spatial mode filter, that consists in a single-mode fiber, and a spectral slit of width Ω. If a photon is detected by the APD, then ideally it has been subtracted from the squeezed vacuum. This subtraction leads to a 1-photon squeezed state, which is very close to a 'Schrödinger kitten' state. Quantum state tomography with a balanced homodyne detector [START_REF] Leonhardt | Measuring the Quantum state of Light[END_REF] allows the complete reconstruction of this highly non-gaussian quantum state of light.

The Wigner function [START_REF] Leonhardt | Measuring the Quantum state of Light[END_REF] was derived assuming that the mode reduction was performed on the mode entering the homodyne detector. To relate our results to the empirical model discussed in the next subsection, we here choose to perform the mode reduction to the signal directly after the DOPA. We model the DOPA using the scheme presented in Fig. 2, where the parameters r and g are linked to the Bogoliubov transform through Eqs. (9)(10)(11)(12). For the calculation of the modes âd,out (j) detected by the APD, the sampling beamsplitter and the mode filters can be separately added to this transform, as explained in Sec. II E. As stated above, here we chose not to include the sampling beam splitter into the mode reduction. The Wigner function [START_REF] Leonhardt | Measuring the Quantum state of Light[END_REF] then describes the signal just after the DOPA. We therefore still need to account for this sampling beam splitter between the DOPA and the homodyne detection, as well as for other losses. For example, one usually accounts for imperfections of the homodyne detection by adding a fictitious beam splitter of transmission η hom just before the homodyne detection, where η hom is the homodyne detection efficiency. Both those beam splitters can be easily implemented by replacing the variances according to:

V (m) x/p -1 = η hom T [V x/p -1]. ( 37
)
and by multiplying A, B and D by η hom T . Before going into detailed calculations, let us first recall the empirical model that was proposed in Ref. [10] to account for experimental results.

B. Empirical model

It is useful to recall the empirical model proposed in Ref. [10] to explain the experiments, and to see by what assumptions our multimode model reduces to it. The DOPA is again modeled as in Fig. 2, producing the same squeezed vacuum. However, in the empirical model it is assumed that the detected photon is either in the homodyne mode with probability ξ, or in an orthogonal mode with probability (1 -ξ). In the latter case, the detection event is not correlated with the homodyne measurement, and one simply performs a homodyne measurement on squeezed vacuum.

The output density matrix obtained with the empirical model is similar to our ρ j in Eq. (19). In fact, the two would be identical if the detected photon was only due to photons in the mode ψ h or from H ⊥ . This is in general not the case, however, as there will be an admixture from â1,out in the photon detection operator. In fact, in order to completely account for the multimode nature of this experiment, the empirical model should be modified in the way depicted on Fig. 4, with the insertion of a beamsplitter of amplitude reflection and transmission coefficients ρ and τ that allows interference between âh,out and â1,out . A photon detection event in such a setup can indeed be equivalent to the application of â † θ (j) to the initial vacuum [see Eqs. (19,21)], provided

ρ τ cosh(r) = 1 tan θ j - 1 tan θ 0 , (38)
with tan(θ 0 ) = β/α. These angles θ 0 and θ j are mixing angles that fix the probability amplitudes of detection of a photon of mode 0 and of mode 1. Our angle θ j in general depends both on the squeezing properties of the light source, and on the filtering of the signal before the photon detector, whereas the empirical θ 0 only depends on the source. Within a narrow-filter approximation that will be detailed in next subsection, such a setup can also account for averaged quantities (29) with the use of an average angle θ instead of θ j (see Eq. ( 54) in the following). This possibility of interference between the homodyne signal and the â1,out signal is the crucial difference between the multimode and the empirical models: More interference makes the empirical model worse. The essential assumption of the empirical model is thus that the photon detection operator does not have a contribution from â1,out . Then γ j and δ j could be replaced by α and β, respectively, according to Eqs. ( 7) and ( 13), and the angle θ j in Eq. ( 19) by θ 0 (in which case, according to Eq. ( 38), the beamsplitter BS(ρ, τ ) in the equivalent model in Fig. 4 can be removed). The DOPA may in general be represented by an ideal single mode DOPA and a two-mode NDOPA as in Fig. 2, and the photon counting is a combination of dark counts as well as a coherent mixture of the two output ports of the NDOPA. Compared to the empirical model developed in Ref. [10], the only difference is the presence of the beam splitter BS(ρ, τ ) which was not present in the empirical model.

The empirical Wigner function can be easily deduced from Eq. ( 27) with the above replacements, and has the same form after the replacement of Ā • • • D by A emp • • • D emp . Coefficients A emp and B emp are obtained by multiplying Ā, B by ξ/ ξ, and replacing γ and δ by α and β, respectively. This gives

A emp = ξ (V x -1) 2 V x + V p -2 , B emp = ξ (V p -1) 2 V x + V p -2 . ( 39
)
The coefficient C emp is given by C emp = 1 -A emp /V x -B emp /V p , and D emp vanishes.

The empirical model produces intuitive results. However, it requires justification. If large spectral slits would be used, then the homodyne mode and many other orthogonal modes would hardly be affected by the slit. If the detected photon could have come from many modes orthogonal to ψ h , then the modal purity ξ would be unacceptably low, and also a large admixture of â1,out would enter the detection signal. Indeed, some of us found experimentally that the spectral slit should be as narrow as possible, while still allowing the detection of a signal, in order to find the highest modal purities (see also Sec. III C). Consequently, narrow slits have been used in the photon-subtraction experiment [10]. Although it is obvious that filtering is necessary, the use of a narrow spectral slit before the photon detector does not make the empirical model automatically valid. A quantitative comparison of both models is therefore needed to test the validity of the empirical model, as given below.

C. Concrete multimode model

Let us now develop a simple spatio-temporal multimode model for which the Bogoliubov transformation can be written explicitly. We assume that light propagation inside the DOPA is described by modes of the form

A( r, t) exp(iω 0 t -i k • r), (40)
where the plane wave is exactly phase-matched, and where the amplitude A satisfies the slowly-varying envelope approximation (SVEA). This approximation does not hold for all the light that exits the nonlinear crystal, but the homodyne mode is supposed to be phasematched, and we will assume that the filters before the APD block the modes that are not phase-matched. We will furthermore neglect diffraction effects within the DOPA. In the basis (x, y, t), where x and y are spatial variables running on the DOPA's transverse plane, the u and v of the Bogoliubov transform (1) then become diagonal (see appendix B), and we have âout (x, y, t) = u(x, y, t)â(x, y, t) + v(x, y, t)â † (x, y, t), (41) in terms of operators that we assume to have commutation relations

[â(x, y, t), â † (x ′ , y ′ , t ′ )] = δ(x -x ′ )δ(y -y ′ )δ(t -t ′ ). (42)
Coefficients in the transformation (41) have the form

u(x, y, t) = cosh[ ql E P (x, y, t) ], (43a) v(x, y, t) = sinh[ ql E P (x, y, t) ]. (43b)
Here E P (x, y, t) is the pump-beam amplitude, which we assume to be real-valued. The parameter q takes into account the nonlinearity of the crystal and l is its length. One can allow for Group Velocity Mismatch (GVM) in the crystal by convoluting E P by a rectangular unit gate of duration τ g , the time separation induced by the GVM after passing the crystal (see appendix B). More precisely, this convolution should be made twice, as we also have GVM for the Second Harmonic Generation (SHG) of the pump beam. The u and v are real-valued functions if the pump beam E P is so, which is a valid assumption if there is no frequency chirp. The homodyne mode ψ h (x, y, t) ≡ e h (x, y, t) will also be taken real-valued in the following. The homodyne signal is then given by âh,out = e h âout , where integration over x, y, t is implied. Mode reduction now starts with the identification

η â0 = dxdydt e h (x, y, t)u(x, y, t)â(x, y, t), (44)
which is Eq. ( 8) specified for our spatiotemporal model. The mode-reduction parameters are now given by spatiotemporal integrals, for example

η 2 = dxdydt e 2 h (x, y, t) u 2 (x, y, t) (45a)
ηα = dxdydt e 2 h (x, y, t) u(x, y, t) v(x, y, t). (45b)

Further parameters can be found analogously.

Averaging over photon detection events.-We have seen in Sec. II E that all averaged quantities can be obtained through the determination of the operator Π defined in (36) when the filters are separately added to the Bogoliubov transform. In the considered experiment two filters are used: a rectangular spectral slit, that can be described using (32); and a monomode fiber that selects a single spatial mode φ s (x, y), which we suppose to be real-valued, and therefore corresponds to the projector φ s (x, y)φ s (x ′ , y ′ ). Note that in this experiment detection times are unknown at the scale of pulses duration, so that there is no time-domain filtering. (In the analysis one should average over all possible photon detection times.) We therefore have to use

Π(x, y, ω; x ′ , y ′ , ω ′ ) = φ s (x, y)φ s (x ′ , y ′ )T (ω)δ(ω -ω ′ ),
(46) where T (ω) = η c R if ω enters into the spectral slit, 0 otherwise. Here R is the sampling beamsplitter reflectivity, and η c accounts for all other losses in the conditioning arm (APD efficiency, optics losses ...). As can be seen from Eq. ( 40), the field amplitudes are defined around a central frequency ω 0 (or 2ω 0 for the DOPA pump beam, see appendix B), so that the frequency ω = 0 for the amplitude Fourier transform corresponds in fact to this central frequency; in that way, a rectangular spectral slit well centered around this central frequency can be defined as ω ∈ [-Ω/2, Ω/2]. The operator Π defined in Eq. ( 46) should be applied to Fourier-transformed mode functions:

φ(ω) = dt φ(t)e -iωt . (47)
Using Eqs. ( 16), ( 46) and (47), the average total photon detection probability per pulse then becomes

P tot = η c RΩ 4π 2 dxdydt v 2 (x, y, t)φ 2 s (x, y), ( 48
)
where the average is taken over the detection modes φ d,j (x, y, t). The expression (48) increases linearly with the filter width Ω, but is valid only if Ω is small enough to warrant the SVEA. Furthermore, when conditioning upon a click in the detector in Sec. II C, we assumed that P tot ≪ 1, an assumption that can now be tested with the explicit formula (48).

In general, the parameters γ and δ appear in W as product sums like γγ * = j |γ j | 2 , or γδ * = j γ j δ * j . In our concrete model, the evaluation of detection-averaged coefficients in the Wigner function involves integrals of the type:

η c R 4π 2 dxdy dx ′ dy ′ Ω/2
-Ω/2 dω f (x, y, ω) h * (x ′ , y ′ , ω).

(49) For example, the average γγ * is found by substituting both f and h in (49) by the time-domain Fourier transform of the function e h u v φ s (x, y, t). [Here and in the following, we abbreviate products like f (x, y, t)h(x, y) by f h(x, y, t).] Since f and h are Fourier transforms of realvalued functions, the corresponding integrals (49) are real-valued as well. Hence, all mode-reduction parameters and coefficients in the Wigner functions are also real-valued. In particular, the coefficient D in W vanishes, see Eq. (28d). There is no difficulty to numerically evaluate integrals (49) and we will do that below, but let us first focus on an additional approximation that can considerably simplify these results, without becoming inaccurate.

Narrow-filter approximation.-We previously discussed the experimental observation that the spectral slit should be as narrow as possible. Another simplification is possible in that case, that simply consists in neglecting in the integrals (49) the frequency dependence of the mode profiles within the narrow width Ω, i.e. f (x, y, ω) ≃ f (x, y, 0) for |ω| < Ω/2. Let us recall that this value ω = 0 corresponds to the central frequency ω 0 of the pulses, where the real-valued amplitudes present a maximum. This has two consequences: first, the presence of this maximum justifies a zero-order Taylor approximation, provided the spectral width Ω is much smaller than the spectral width of the pulses (typically 2π/τ , where τ is the pulse duration); second, if all functions involved in (49) present a maximum at ω = 0, then the narrowfilter approximation generates an upper-bound for these integrals, and therefore for quantities like |γ 2 |, |δ 2 | or the modal purity ξ (see 20,29). This approximation will be applied and tested in Sec. III D, dedicated to the numerical results. This approximation brings the following simplification in the integrals (49):

= η c RΩ 4π 2 dxdy f (x, y, 0) dx ′ dy ′ h * (x ′ , y ′ , 0) (50) = η c RΩ 4π 2 dxdydt f (x, y, t) dx ′ dy ′ dt ′ h * (x ′ , y ′ , t ′ ).
Evidently, we end up with separate integrals over f and h, and using the definitions ( 14) we obtain the averages

γ = √ η c RΩ 2π 1 η dxdydt e h u v φ s (x, y, t) (51a) δ
= - αγ β + 1 β √ η c RΩ 2π dxdydt e h v 2 φ s (x, y, t). (51b)
In the narrow-filter approximation, averages of products are simply given by products of averages, γγ * = |γ| 2 , and γδ * = γδ * , etc. Essentially in the Ω → 0 limit the filter removes any temporal information about the time the photon was emitted from the DOPA. The photodetection then corresponds to a single mode with ω = 0, regardless of the average over detection times. We therefore find for the photon-subtracted squeezed state an average Wigner function of the form (28), with coefficients

Ā = P -1 tot |γ(η + α) + δβ| 2 (52a) B = P -1 tot |γ(η -α) -δβ| 2 , (52b)
and C = 1 -Ā/V x -B/V p and D = 0. Using the same substitution in (20,21) we can also introduce the averaged modal purity

ξ = |γ| 2 + | δ| 2 P tot (53)
and the average angle θ defined by:

tan θ = δ/γ. ( 54
)
Constant profiles.-Before dealing with a more realistic case, it is interesting to focus on the case of constant profiles. Let us assume a constant value for e h , E P , and φ s within a space-time support of volume Ξ. The normalization of the homodyne mode implies e 2 h = 1/Ξ, and equation (43,45) leads to η = u, α = v, and β = 0. As β vanishes, the mode â1 is no more defined, and Eqs. (14b,51b) cannot be used anymore. In fact the mode-reduction procedure now leads to an effective single-mode model rather than a two-mode model. The homodyne mode is now in the single-mode space H 1 spanned by â0 ,â † 0 . One can simply put δ j = 0 in (13), and hence δ = 0. This leads to the average angle θ = θ 0 = 0. The modal purity becomes ξ = 1, as it should for a single-mode model. Most importantly, we find C = -1, which according to Eq. ( 28) corresponds to the most negative value for the Wigner function at the origin, W (0, 0) = -1/π.

So, with constant profiles and a narrow filter slit, we recover from our multimode model the single-mode description for photon subtraction experiments. Since this limit leads to the most negative Wigner function, it represents the ideal limit for producing states for QIP applications, at least according to our simple multimode model. This shows that the multimode nature essentially appears through the mode distortions due to the non-constant space and time profiles of the pulses. The central role of gain-induced distortions is particularly clear with regard to the multimode nature of the squeezed vacuum produced by the DOPA: assuming a constant pump field, that is assuming no gain-induced distortions, is in fact enough to obtain β = 0. Let us now return to a more realistic model, taking into account these profiles.

Gaussian profiles.-As a more realistic simplification, we assume gaussian profiles for the various fields. For instance, we write the homodyne field e h as e h (x, y, t) = e h,0 exp -

x 2 + y 2 w 2 exp -2 t 2 τ 2 , ( 55
)
where w is the beam waist, τ the duration of the gaussian pulse, and e h,0 a normalization constant. The pump beam is usually obtained by SHG, in a crystal pumped by a beam identical to the homodyne beam. In the lowest order of the SHG process, the profiles of E P and e 2 h have the same shapes, so we can assume another gaussian profile:

E P (x, y, t) = E 0 e P (x, y, t) (56) = E 0 exp - x 2 + y 2 w 2 P exp -2 t 2 τ 2 P ,
where one expects the pump pulse duration to be τ P = τ / √ 2. However, if GVM is taken into account, this gaussian profile (56) must be convoluted by rectangular gates. In practice, such convolutions lead to beam profiles that are still very close to gaussians.

The final gaussian profile to be introduced here is the spatial mode φ s (x, y) of the filter in front of the APD. It is the LP 01 mode of a monomode fiber that is well approximated by a normalized gaussian of waist w f .

D. Numerical results

Here our goal is twofold: first, to compare our multimode analysis with the empirical model that was used before to analyze photon subtraction experiments. Second, by exploring our multi-parameter multimode model, we look for parameter regimes that are best suited for producing states with the most negative Wigner functions.

Expansions in pump field.-For our numerical work it is convenient to write all fields as Taylor expansions in the pump field E P . From Eq. ( 43) it follows directly that

u 2 (x, y, t) = 1 2 + 1 2 cosh[2qlE P ] = m b m e m P , (57a)
u v(x, y, t) = 1 2 sinh[2qlE P ] = m c m e m P , (57b)
v 2 (x, y, t) = - 1 2 + 1 2 cosh[2qlE P ] = m d m e m P , (57c)
with E P = E 0 e P (x, y, t) as in Eq. ( 56). This defines the constant coefficients b m , c m , and d m . For qlE 0 < 1 these expansions converge quite quickly. We then only have to insert relations (57) into the various integrals for an efficient numerical evaluation. For instance, we can rewrite (45a) as:

η 2 = m b m P m (58)
with

P m = dxdydt e 2 h e m P = 2 √ 2w 2 P τ P (mw 2 + 2w 2 P ) mτ 2 + 2τ 2 P . (59)
In the same way we have

ηα = m c m P m , (60a)
P tot = η c RΩ m d m Q m , (60b) γ
= √ η c RΩ η m c m R m , (60c) δ
= - αγ β + √ η c RΩ β m d m R m , (60d)
with

Q m = 1 4π 2 dxdydt φ 2 s e m P = π -3/2 τ P w 2 P 2 √ 2m(mw 2 f + 2w 2 P ) (61a) R m = 1 2π dxdydt φ s e m P e h (61b) = π -3/4 √ τ τ -2 + mτ -2 P 1 ww f (w -2 + mw -2 P + w -2 f )
.

After fixing parameters, these expansions in the pump field can be readily used for numerical evaluations.

Fixing basic parameters.-First we fix some parameters of our multimode model in order to present numerical results and to see how much our analysis differs from the one in Ref. [10], where filtering before the photon detection was not modeled explicitly. We take w = 1.2w P and a transmission T = 90% of the sampling beam splitter. Moreover, we fix w f = w/1.5, which is compatible with the coupling efficiency into the filtering monomode fiber (approximately 80%, see Ref. [10]).

Regarding efficiency of homodyne detection, the mode âh considered in Sec. II was defined as the mode that perfectly matches the local oscillator of the homodyne detection, in other words the matching efficiency equals unity by definition in our model. The transmission of the optics and the photodetection efficiency together lead to an overall efficiency of homodyne detection η hom . We put η hom = 0.93, in agreement with Ref. [10].

As stated above, if GVM is taken into account, the almost gaussian profile of the pump pulse is convoluted twice by a rectangular gate with time window τ g . A KNbO 3 crystal of length l = 100 µm has τ g = 120 fs. For an initial duration of the homodyne pulse τ ≈ 150 fs, the convolutions indeed lead to a nearly gaussian beam profile with τ P ≈ τ . We assume the identity τ P = τ in the following.

Negative Wigner functions.-As stated in the Introduction, the global minimum of a Wigner function is the standard figure of merit for the nonclassicality and 'nongaussianity' of the corresponding state. After subtraction of a single photon, the Wigner function W (x, p) is always most negative in the origin (since D = 0). the low-squeezing limit s → 1. This can be understood as there is less gain-induced distortions in that case. In Ref. [10], the best experimental results (highest modal purities) were obtained for s = 0.56. For this value of s, which can be selected by choosing the right value for the quantity qlE 0 , we obtain g = 0.50 and W (0, 0) = -0.034; the latter value is close to what was observed in [10], without correction for the detection efficiency. At this stage, it can be interesting to compare this result, obtained using the narrow-filter approximation, with a more accurate calculation based on a complete evaluation of integrals (49). Figure 6 presents the numerical results obtained for W (0, 0) at s = 0.56 as a function of the spectral slit transmission for the homodyne mode. (This transmission can be increased by making the spectral slit width Ω larger.) A minimal value is clearly reached for low transmissions, justifying a posteriori the use of narrow spectral slits in the experiment of Ref. [10]. Since for low transmissions, W (0, 0) does not differ much from its minimal value, the narrow-filter approximation that we made in Sec. III C gives accurate results. Figure 7 predicts the behavior of W (0, 0) when varying the size of the pump beam. Experimental values for the widths were related by ω = 1.2ω P [10]. Fig. 7 clearly shows that one can await a high increase of the negativity from a larger pump beam. This result was intuitive, as there is less gain-induced distortions in that case, but is here quantified. This can motivate the use of amplified pulses [START_REF] Dantan | [END_REF], in order to have a spatially broader pump beam (i.e. with larger w P ), but with the same intensity.

Comparison with empirical model.-Measured negative Wigner functions were interpreted in [10] using the empirical model as introduced in Sec. III B, where the photon is subtracted in the 'good' mode âh,out with probability ξ, and where the state is left in the squeezed vacuum with probability (1 -ξ). As explained before, the main difference between the empirical and our models is that the conditioned state in the former corresponds to the single-photon initial state â † θ0 with tan θ 0 = β/α, while it corresponds to â † θ in the latter, with tan θ = δ/γ. These angles θ 0 and θ are mixing angles that fix the probability amplitudes of detection of a photon of mode 0 and of mode 1. Figure 8 shows θ 0 and θ as a function of the squeezing parameter s. Clearly, θ 0 and θ do not differ too much, and by less than 10% for s = 0.56.

Another important difference between our model and the empirical model is that modal purities ξ j in our model are fixed by the relation (20), whereas the parameter ξ in the empirical model is a free parameter. This freedom can be used to fit the data, i.e. to have A emp = Ā or B emp = B. It is however not a priori possible to fit both parameters Ā and B (52a,52b) from (39) using only the fitting parameter ξ. In neither model should the variances V x and V p be considered as free fitting parameters of the photon-subtraction experiment, at least their values should agree with the values for V x,p obtained by homodyne measurements of the squeezed vacuum.

In our model the mixing angles θ j and their average θ take into account the filtering of the signal that is used for conditioning. In the empirical model, the corresponding angle θ 0 is independent of the filtering. Thus it is to be expected that this inaccuracy of the empirical model will lead to optimally fitted modal purities ξ opt in the empirical model that are systematically lower than the average modal purity ξ in our model. This is indeed what we find for the curves in Figure 9: the best fit in the present example is obtained for ξ opt = 0.87, a value that is indeed smaller but still close to ξ ≈ 0.91. The high quality of this fit (with an error less than 1.2%) is directly linked to the fact that in the present case θ 0 ≈ θ.

There is a possibility to improve this result if g is considered as a fitting parameter as well. We obtained an error of less than 0.4% between the Wigner functions for ξ opt = 0.89 and g opt = 0.53, i.e. for a value of g that differs by 6% from the value given by the multimode model. In other words, if g has a great influence on A emp and B emp in (39), it has a very low impact on V x , V p ; the change from g = 0.50 to g opt = 0.53 modifies the values of V x , V p by only a few 10 -3 , and for this reason it is very difficult to accurately measure g from squeezed vacuum [START_REF] Wenger | [END_REF]. These considerations explain why the empirical model can fit experimental data so successfully; even when θ 0 is not equal to θ, the parameter g gives a supplementary freedom for fitting.

IV. DISCUSSION AND CONCLUSIONS

We have introduced a straightforward and physically intuitive procedure that we call 'mode reduction' to simplify the multimode description of squeezed light to the bare essentials. For photon-subtraction experiments, this means that the homodyne signal is reduced to an effective two-mode description and the detector signal requires one extra orthogonal effective mode. We derived the Wigner function of the homodyne signal conditional upon the detection of a single photon, and we also showed how to average over possible measurement outcomes.

The general mode-reduction formalism was then applied to a detailed model describing photon subtraction of gaussian spatiotemporal pulses of squeezed light. This model features many experimental parameters such as beam waists and duration of the pulses that can be independently measured. Indeed, our model does not have free fitting parameters. This allows one to study in detail what are the crucial experimental parameters to produce optimally negative Wigner functions with pulses of squeezed light.

We compared our new model to the empirical model that was used before to analyze photon-subtraction experiments in [10]. In fact, the formulae for the output Wigner functions look similar. One crucial difference is that the empirical model does have a free parameter, namely the quantity called the modal purity. In our model modal purities also occur, be it with a slightly different meaning, but they are fixed quantities. A good agreement between our model and experiments therefore gives more understanding than an accurate fit with the empirical model.

We found that in the range of parameters of the measurements in [10], both our model and the empirical model are accurate. We reasoned that modal purities in our model would be systematically higher, and in our numerical example we found this to be the case. The accuracy of the empirical model strongly depends on the availability of the free parameter. It was nevertheless a surprise in the theoretical analysis that the mixing angles θ 0 and θ, describing the relative probability of measuring a photon in either one of two effective modes, differ at most 20% in a whole range of squeezing parameters.

Our mode-reduction procedure is closely related to the analysis of photon-subtraction experiments of Refs. [19,20,21]. One could express our mode-reduction parameters in terms of elements of the covariance matrix of Refs. [19,20,21]. Our output Wigner function in Eq. [START_REF] Leonhardt | Measuring the Quantum state of Light[END_REF] then reduces to the one in Ref. [21], but only in the special case that all our mode reduction parameters are real-valued so that D j in Eq. (28d) vanishes. This we assumed for simplicity in Sec. III C. Our mode-reduction procedure is carried out in the Heisenberg picture. We think that our approach has some advantages. In our approach it becomes quite intuitive in what sense it goes beyond the empirical model of Ref. [10]. In our concrete multimode analysis, we include effects not considered in Ref. [21], such as the transverse beam profile, for which we found that wider pump beams lead to more negative Wigner functions.

In conclusion, we presented a very concise model that can account for the multimode nature of projective photon-counting measurements. It gives an intuitive picture of photon-subtraction experiments, close to the empirical model previously published. This multimode model therefore gives consistent results, in agreement with previously published experiments where pulses of light with negative Wigner functions were produced conditionally. Our model can be used to predict the changes in the output upon variation of experimentally relevant parameters, and to optimize the setup design.

where the phase-matching condition k 2ω0 = 2 k ω0 is assumed to be satisfied. By using the SVEA in Maxwell's equations, neglecting diffraction terms and considering the first-order dispersion, we obtain ∂A ∂z + 1 v g,ω0 ∂A ∂t = qE P (x, y, t -z v g,2ω0

-δt)A * , (B3)

where A = A( r, t). The substitution of t -z/v g,ω0 by t then leads to ∂A ∂z ( r, t) = qE P (x, y, t -Dz -δt)A * ( r, t), (B4)

where D = v -1 g,2ω0 -v -1 g,ω0 is the GVM. With the assumption that E P is real-valued, the solution to Eq. (B4) becomes

A out = cosh[ qlF P ] A in + sinh[ qlF P ] A * in , (B5)
where the (x, y, t)-dependence was suppressed. The effective pump field F P is given by

F P = 1 l l 0 dzE P = 1 τ g τg+δt δt
dτ E P (x, y, t -τ ), (B6) with τ g = Dl the time separation induced by the GVM after crossing the crystal. Eq. (B6) shows that the effective pump field F P is a convolution of the pump field E P with a rectangular unit gate of duration τ g , which for δt = -τ g /2 is centered around the origin. As the quantized version of Eq. (B5), we then find Eq. (41) of the main text. The pump field E P of the main text is to be understood as the effective pump field F P derived here.

Evidently, F P → E P in the limit τ g → 0 (no GVM).

FIG. 3 :

 3 Photon subtracted state homodyne detection

FIG. 4 :

 4 FIG.4: When one applies the mode reduction procedure to the multimode model for photon subtraction as sketched in Fig.1, the resulting Wigner function is equivalent to the twomode model depicted in the figure. The DOPA may in general be represented by an ideal single mode DOPA and a two-mode NDOPA as in Fig.2, and the photon counting is a combination of dark counts as well as a coherent mixture of the two output ports of the NDOPA. Compared to the empirical model developed in Ref.[10], the only difference is the presence of the beam splitter BS(ρ, τ ) which was not present in the empirical model.

FIG. 5 :FIG. 6 :

 56 FIG.5:(color online). Minimal value of the Wigner function, W (0, 0), as a function of squeezing parameter s = exp(-2r), which is varied by changing the quantity qlE0. The narrowfilter approximation was made for the spectral slit. Fixed parameters: pulse parameters w = 1.2wP , τP = τ = 150 fs, transmission of sampling beam splitter T = 90%, efficiency of homodyne detection η hom = 0.93.

FIG. 7 :

 7 FIG. 7: (color online). Minimal value of the Wigner function,i.e. W (0, 0), as a function of wP /w, where wP is the waist of the pump field and w is the waist of the homodyne field. qlE0 is fixed such that squeezing parameter s = 0.56. Other parameter values as in Fig.5.

FIG. 8 :

 8 FIG.8:(color online). Mixing angles θ (solid line) and θ0 (dashed line) of modes 0 and 1, as a function of squeezing parameter s, which is varied by changing qlE0. Other parameter values as in Fig.5.

FIG. 9 :

 9 FIG.9:(color online). Sections along (x, 0) and (0, p) of the Wigner function W (x, p) (solid line), and corresponding best fits using the empirical model (dashed line) with ξopt = 0.87. Parameter values as in Fig.5. Results of our model and the empirical fits almost overlap.

Acknowledgments

We thank K. Mølmer for useful discussions. This work has been supported by the Danish Research Council through QUANTOP, by COMPAS, and by the Niels Bohr International Academy.

APPENDIX A: WIGNER FUNCTION

In order to derive Eq. ( 27) from Eq. ( 26), one has first to invert Eqs. (7,22), leading to â0 = η(â h,out -βâ † 1,out

One should then replace in Eq. [START_REF] Wenger | [END_REF] x 0,1 , p 0,1 by x h,out , p h,out , x 1,out , and p 1,out , and integrate over x 1,out , p 1,out . It is however convenient to make a change of variables so that the integral is over x 1 , p 1 instead of x 1,out , p 1,out . In this case the only transforms needed for this calculation is

as well as the transformation of the integral

One should then note that the Wigner function ( 26) is the product of a polynomial in x, p, and of a gaussian term exp (-R 2 ), with

) With Eq. (A3), the exponent can be rewritten as

The integral (A4) with Eq. ( 26) as its integrand can then be found by replacing in the integrand the squares x 2 1 and p 2 1 by

and by replacing the first-order terms according to

APPENDIX B: SLOWLY VARYING ENVELOPE APPROXIMATION

The goal of this appendix is the derivation of the local Bogoliubov transformation (41). We assume that inside the DOPA the pump pulse with an angular frequency 2ω 0 travels at a speed v g,2ω0 , with negligible absorption. This field can therefore be written as i

where δt is an arbitrary time delay and where 'i' is a purely conventional phase factor. Let us write the probe beam as E( r, t) = A( r, t) exp(iω 0 t -i k ω0 • r), (B2)