The surprising persistence of time-dependent quantum entanglement

The mismatch between elegant theoretical models and the detailed experimental reality is particularly pronounced in quantum nonlinear interferometry (QNI). In stark contrast to theory, experiments contain pump beams that start in impure states and that are depleted, quantum noise that affects -- and drives -- any otherwise gradual build up of the signal and idler fields, and nonlinear materials that are far from ideal and have a complicated time-dependent dispersive response. Notably, we would normally expect group velocity mismatches to destroy any possibility of measurable or visible entanglement, even though it remains intact -- the mismatches change the relative timings of induced signal-idler entanglements, thus generating"which path"information. Using a"positive-P"approach ideally suited to such problems, we show how the time-domain entanglement crucial for QNI can be -- and is -- recoverable despite the obscuring effects of real-world complications.

Quantum entanglement is important because it plays a key role in a range of quantum devices, notably in induced coherence [1][2][3] based quantum imaging / QNI schemes [4][5][6][7][8][9][10][11]; and in the time domain is a subject of wide-ranging and active study [12][13][14][15][16]. However, the complicating effects of material dispersion in the entanglement-generating nonlinear medium, or during subsequent propagation, are typically not considered [17][18][19]. This has wider relevance, not only for quantum interference in general, but also e.g. in quantum data transmission [20][21][22] and QNI. We demonstrate in this Letter how and why, despite the complete de-synchronization of entangled fields caused by material dispersion, a slow detection process can perform an unexpected "entanglement recovery", so that time resolved QNI can, after all, unexpectedly succeed.
Our testbed for examining the limitations on time-resolved quantum measurements is a pulsed QNI system where time dependence is relevant for all field interactions, in particular with the material dispersion (i.e. group velocity, and group velocity dispersion (GVD)) present alongside the entanglementgenerating nonlinearity. In ghost imaging, for example, one can imagine a clear separation between standard (spatial) schemes [23,24] and temporal schemes [25,26], but if material dispersion was present during propagation, such simplicity would be disrupted. To address such intrinsic complications requires a shift in both theoretical methods and mindset; a description can no longer rely on using only a small number of possible states (typically Fock states), and judgements based on path indistinguishability or phase shifts. Instead, a set of time dependent states is required, and consequently it is not path lengths but relative timings that matter.
How the fields propagate and are transported through the QNI system is shown on figure 1; the layout is very similar to that of Kolobov et al [27]. The key feature of the nonlinearity is that it produces correlated signal and idler photon pairs from an incoming pump field; this is often achieved using a χ (2) interaction, but here we use a degenerate-pump four wave mixing (FWM) (see e.g. [28] 1. A representation of the QNI scheme (e.g. [27]). The two identical pump fields are in blue, signal fields in green, and idler in red; here the FWM process combines two pump photons to create a signal-idler pair. The resulting signal fields then interfere at the beamsplitter (BS) and detected at D A and D B . In an Imperial experiment [28], the nonlinear stages (NL1 and NL2) take place in opposite directions through the same length of fibre.
first nonlinear stage (NL1) and generates entangled signal and idler fields; and whilst the signal is diverted to the final beamsplitter, the idler instead interacts with the to-be-imaged object, and then serves as a co-input, with a copy of the original pump field, to the second nonlinear stage (NL2). The signal field departing NL2 is thus influenced by an idler field entangled with the first signal field, and this information is extracted by interference at the beamsplitter, before photon detection. The nonlinear propagation model in our simulations is a well-established one derived originally for the prediction of squeezing generation in optical fibres [29][30][31][32]36], has also been used to model multi-field parametric processes [37], and here centres primarily on multi-field nonlinear propagation through a dispersive material. It also includes stages representing setting up the initial conditions, interaction with an object to be imaged, and mixing at a beamsplitter at the interferometer output. In this work we use an established off-diagonal coherent state basis positive-P approach [38][39][40], that enables both group velocity and dispersion to be easily implemented in a numerical scheme [29]; crucially, its off-diagonal nature allows a complete representation of the full quantum mechanical density matrix of the system, and of its dynamics.
Our description uses time-dependent quasi-conjugate pump field amplitudes α u (t), α † u (t), which are only complex-conjugate on average; likewise for signal α s (t), α † s (t), and idler α i (t), α † i (t). There is also a set of independent timedependent noise increments {dW a (t)}. In simulation, time is discretized and labels a set of finely divided, sub-pulse length modes where each field is held as an array of sequential "timebins" at t ∈ {t j }, each of which contains pairs of complex field amplitudes; these interact and evolve as the fields propagate step-by-step through space. Since these time-bins have a simulation-specific duration, field averages such as <α † s α s > represent intensity (in photons per second), not photon number, and normalizations and parameter scalings reflect this.
This field information is held along with fixed parameters for the nonlinear coupling κ, losses γ m , and dispersive properties D m (ω), where the field subscripts, as above, are m ∈ {u, s, i, r}. Since this is a stochastic technique, very many independent copies of the evolution need to be run (for the simulations here, typically in the range 10 5 to 10 7 ), and the necessary ensemble averages taken.
This propagation model is very close to previous ones (e.g. [29,32,37]), but here we have three co-propagating fields and a different nonlinear interaction; namely one for the degenerate FWM process with only resonant wave mixing terms where 2ω u = ω s + ω i . Since self-phase modulation (SPM) terms are not significant in the low-power regime (see e.g. [28], the interaction Hamiltonian we need here is simplŷ This interaction term results in propagation equations which are best expressed in the incremental form suited to such stochastic differential equations (SDEs) [38]. In an appropriate co-moving frame, including loss along with nonlinear effects [32], but suppressing the time argument common to all fields α m , α † m and quantum noise increments dW a , we have These equations are used to update each time-bin in the temporal profile of the fields as they propagate (step) forward in space. Deterministic evolution terms are in square brackets [...], and prefactors for stochastic (noisy) terms in braces {...}.
Here we see that there are both coherent interactions between the three fields, and correlated nonlinear quantum noise terms. The noise increment dW 1 drives both α s and α i , whilst dW 2 drives α † s and α † i , ensuring the pairs are correlated but not complex conjugate. The comparable classical model (or even a semi-classical model, see e.g. [33][34][35]) would have only three equations and no noise terms.
Material dispersion is the other key feature, and we interleave it with the nonlinearity in a split-step scheme (see e.g. Carter et al. [29,32]), using linear phase shifts in the spectral domain for group velocity, and quadratic shifts for GVD.
Material parameters in the simulations are chosen to be compatible with fiber-based photon pair sources based on spontaneous four-wave mixing [28], which use about 100cm of Thorlabs PM780-HP fibre (see Table I). For the simulation, we convert parameters into units based on meters (m), picoseconds (ps), and field excitation amplitudes referenced back to photon numbers per picosecond. A crucial step here is to consider pulsed operation on a picosecond scale, where the group velocity mismatches are significant. This regime is where the time-dependent nature of the entanglement will be most exposed to the disruptive effects of material dispersion.
Since the effect of group velocity mismatches and GVD turns localised entanglements into temporally distributed ones in our simulations we see a fan of signal amplitudes that spreads out behind the pump pulse, whilst a fan of idler amplitudes spreads out before. Thus, crucially, the signal-idler entanglement is not only distributed over a range of times, it is also between different times; i.e. it is a multi-time correlation. For best imaging visibility, this requires careful synchronisation at NL2, where the first-generated part of the idler in NL2 should be coincident with the pump pulse as it enters.  [28]. Since the fibre is weakly guiding, these are based on those for bulk silica. The nonlinearity in SI units is n 2 = 3 × 10 −20 m 2 /W, and the loss is γ = 0.004/m. In simulation units, the PM780-HP's stated transverse field mode area of about 25µm 2 and the pump photon energy of 25 × 10 −20 J, means that the rescaled nonlinearity is n 2 = 0.3 × 10 −15 per photon-picosecond. In a time-bin T ps long, the pump photon has a power ∼ (0.25/T ) µW, so that a power of 100W corresponds to a flux of 4 × 10 8 photons/ps. Estimates for signal and idler fields are similar. Note that the pump powers used in the simulations are increased to enable good simulation statistics (see Supplementary information).
Objects placed in the idler beam leaving the NL1 stage disrupt the entanglement with the signal beam, and that disruption changes the numbers of detected photons, enabling the object's presence or properties to be inferred. However, in non-imaging contexts, we can view them as representatives of further disruptive unwanted real-world effects: loss, phase shifts through optical elements, or extraneous couplings. Here we consider passive objects that impart (a) a phase modulation ∆φ of the idler field, or have (b) a reflectivity r that reduces the idler amplitudes as they pass (as in [27]); so r = 1 removes all entanglement. We also consider (c) imaging of dynamic objects, a key feature since an object's time dependence will affect visibility, just as timing, group velocity, and GVD do.
Our linearly coupled dynamic objects with amplitudes at the position z specific to the object; and a field-object interaction strength η. At the object we set the initial conditions at Now that the incident field α m , α † m has excited the dynamic object, this excitation acts back on the field and modifies it. We therefore then update α(t), α † (t) according to the same linear interaction but here integrated forward in space, using and keeping the SDE notation for consistency. Note that this could be extended to allow for nonlinear couplings or dynamics with the object, or to even use (e.g.) a two-level-atom or Raman models (see e.g. a classical counterpart [41]).
Detection and Visibility: The photon number rate measured at the detectors is taken to be the ensemble-averaged photon number of the relevant field at that point. An improved detector model could be implemented using similar approach to time-dependent objects. An important quantity is the visibility of the entanglement, which is the difference of the detector counts ("signal") divided by the sum of the detector counts ("background"); i.e. |n A − n B |/(n A + n B ). To suppress sampling artifacts in the pulse wings we add to the background an offset of 0.05% of its maximum.
We test the basics of the simulations with a simple CWequivalent parameter set with no group velocity or GVD effects, and look at how the visibility varies with object phaseshifts and object reflectivities. Here, correlations between the first signal's time-bins and the second signal's ones will always be synchronised, maximizing the visibility. For sufficiently low signal-idler generation efficiencies, we should see equal photon number intensities in the signal's field-modes, but distinct photon number intensities after the beamsplitter, i.e. at the two detectors; this is clearly shown on fig. 2.
Time averaging is a crucial part of any detection model used here, since real detectors are very slow (typically 100ps) when compared to the temporal resolution of our simulations (∼ps). Although more sophisticated models can easily be imagined, here we simply sum the time-binned amplitudes of each field over some chosen m-bin detector response time, before ensemble averaging to get the detected photons: . (11) Crucially, this averaging process in the detector helps expose correlations between time-bins that have become offset due to dispersion mismatches. Note that this summation is essentially the same as the process for combining multiple short time bins into a longer one; just as we need to do in numerical convergence checks. On fig. 3, the entanglement is always fully present, but has its visibility reduced by v g or GVF mismatch. However visibility can be recovered using longer averaging intervals; at least up to a (v g ) cut-off when the time difference exceeds the averaging windows. Pulsed simulations: Here we standardise on an input 40ps pump pulse and a 512ps window divided into 2ps bins. Simulation pulse intensities were chosen as low as practicable, given the constraints of computation time and the requirement for good simulation statistics. Further, the pump-idler pulses were ideally synchronised as they entered the NL2 stage. However, the signal fields are also mis-timed at the output beamsplitter by δ τ s = 20ps (i.e. 10 time bins). This not only mimics imperfect experimental setup, it is also useful in providing an example where the detector averaging has more to recover. We also use a technique, described in the Supplementary Material, to reduce the effect of sampling error -something which would otherwise make positive-P simulations either problematic or computationally prohibitive.
Results for no object and standard material parameters, but also considering artificially reduced group velocity offsets, are shown on fig. 4(a-f), where they are compared with full parameter results. The displacement to negative t of the visibility peak is a result of the group velocity walk-off of the signal field. We see that the detected photon rates n decrease at larger group velocity mismatches -this is due to the increased spreading of the generated fields, and hence less effective nonlinear generation. Despite the decrease in detected n, we see that the reduction in detector-averaged entanglement visibility is relatively minor; whilst in stark contrast, the drop in unaveraged visibility is significant. Thus we see that sufficiently long averaging interval allows us to recover most of the max- imum possible visibility, albeit not all; and as we would expect the averaging also helps reduce the significant statistical variation visible on the un-averaged data. Thus fig. 4 shows that group velocity mismatches are not as problematic as they might at first appear, since the generated entanglements, however scrambled they might be by the gradual nonlinear generation and significantly dispersive propagation, can be recovered to a significant extent by time-averaging at the detector. In fig. 4 the simulations ensured the idler pulse arrived at the NL2 stage in correct synchronization with the pump pulse. We can see the effect of mis-timing the idler pulse at this point on fig. 5(a), where at least for these parameters -notably 40ps pulse widths -the fall off in averaged visibility is gradual. This is due to a combination of the pump pulse length (40ps), the averaging time (64ps), and the group velocity spreading (∼ 80ps). On fig. 5(b,c), and in broad agreement with trends in fig. 2, we see a fall-off of the time-dependent entanglement visibility with object phase depth and object reflectivity.
Dynamic objects further emphasise the potential role of time-dependence. Fig. 5(d) shows recovered entanglement visibility values for a range of interaction strengths η . At low η , idler field excitations are coupled into the object but not out again, leading to reduced visibility. However, as the η increases even further, those excitations can also start being coupled back out, leading to a partial recovery. This nontrivial behaviour suggests the possibility of interesting tradeoffs when considering the imaging of dynamic objects. In conclusion, we have shown how the time-averaging process inherent in slow detector response times enables recovery of the entanglement necessary quantum nonlinear interferometry, even when confounding effects such as group velocity differences, dispersion, mis-timing of pulse arrivals, or objects interposed in the idler beam are allowed for. Indeed, group velocity mismatches would normally be expected to play a critical role, seeing as they can rapidly de-synchronise mutually entangled time-slices of the light field. This can be seen in our simulation results (e.g. in figs. 3(a) and 4(f)) which show a significant improvement when changing from un-averaged to averaged visibilities. In contrast to visibilities, which are a ratio, nonlinear generation efficiencies suffer penalties from the effect of material dispersion regardless of averaging. Finally, our linearly coupled dynamic object model acts as a starting point for more sophisticated interactions; a feature likely to be important in systems relying on short pulses, resonant behaviour, and time-domain entanglement.

A. Representation of entanglement
In the model used here, entanglement is represented as classical statistical correlations between complex field amplitudes which are capable of reproducing the off-diagonal elements of the density matrix -i.e. they can represent all the necessary quantum properties [38]. This is possible because each field is represented by two independent amplitudes α m and α † m ; and although they are complex conjugate on average, i.e. < α m >=< α † m > * , in any given trajectory making up part of the large ensemble, they need not be.
By looking at the SDE's for the nonlinearity (2) - (7), we can see that the non-daggered and the daggered amplitudes are driven by different noises. Thus both the mean photon numbers < α † s α s > and < α † i α i > could even remain nearly zero even whilst a significant quantum statistical correlation (i.e. entanglement) is being created between the signal and idler fields; i.e. between α s and α † s , and between α † s and α † i .

B. Simulation statistics: reducing sampling error
It has long been known that that getting good simulation statistics with the fully quantum mechanical positive-P representation can be challenging [32], so that it is very common to resort to the much simpler, but approximate, truncated Wigner representation [40] in simulations. In cantrast to the positive-P, the truncated Wigner representation is essentially a semiclassical hidden variable model that represents the quantum uncertainty as simple statistical fluctuations in the field amplitudes [33,34]. This halves the number of equations required by the simulation (needing only a single complex α rather than the double α and α † ), reducing the state space, and resulting in an effective and sufficiently accurate method when e.g. studying quadrature squeezing [42], and its generation using optical pulses in nonlinear materials [29,30,32,36].
However, in this work, here we want to ensure an accurate representation of the quantum effects, and so we stay with the full positive-P model (cf the case of quantum tunnelling [33,34] and nonlinearity and the quantum vacuum [35]). Since the subtler effects of quantum entanglement is addressed here, rather than the simpler quadrature moments, a truncated Wigner representation would not suffice, since it imposes an unavoidable linkage between correlations and photon number. This leaves us requiring the use of a full positive-P description and concommitant extremely long run times. This is especially problematic since we may need to resolve very low average photon numbers inside an extremely noisy background.
We address this using a hybrid strategy which allows us to use just one very high ensemble number M B simulation to get a good estimate of the background for some particular case, and then adjusting this using two more simulations with lower ensemble numbers M R but perfectly matched random noise generation. We call the M B simulation the "background", and the other two the "reference" and "target" simulations. The reference simulation has identical parameters to the background simulations, and the target simulation has the parameters coresponding to the particular result we are interested in. The difference between the reference and target simulations only depends on the differences between the system parameters, with only "second-order" noise effects -resulting from how the noise influences propagation differently, and not ii from different random number sequences.
When trying to evaluate a photon number n for some chosen target parameters, we proceed as follows, using the usual notation where a statistical average is denoted < n >, but additionally adding a subscript to denote the ensemble size, with ∞ to denote the idea infinite-ensemble case. Each trajectory in the background or reference ensembles returns a value n j , and each in the target ensemble returns n j .
Thus for n we can have either a low sampling error, or a larger sampling error, as follows with the sampling error reducing for each as M B and M R are increased; thus for large enough M values we have i.e.
but noting that the noise-sampling error in this approximate equality (15) is dominated by the larger variation in the smaller reference ensemble.
Similarly, for n we we have Since (15) should average to zero, we can now write which, given its dependence on both n M R and n M R , would at first appear to suffer a larger sampling error based on the smaller M R , not a small sampling error based on the large M B . However, the key point is that since we have exactly matched the simulation noise values between the reference simulations of n j and target simulations of n j then their difference is due to differences of the trajectory dynamics between the two simulations, and is only weakly dependent on the specific noise values. Since this cancels out the bulk of the sampling error due to the noise in these smaller-ensemble simulations, we can now write where now it is the sampling error from the background ensemble that dominates. Thus for any parameter set, we only need to do one large ensemble M B background simulation, one smaller ensemble M R reference simulation otherwise identical to the background one, and then then many small M R target simulations which do have a parameter variation compared to the background simulation. Background parameters were of systems with perfectly transparent objects, and our target simulations varied only the object properties. This meant that all the results shown on fig. 5 could be generated from the same background simulation, a large and -in our case -very necessary reduction in total simulation time.
Here we typically did background simulations with M B sizes from ∼ 10 6 ( fig. 4(a,d)) to ∼ 16 × 10 6 (figs. 4(c,f) and 5), and reference and target simulations with M B = 128 × 10 3 . In the most extreme case, this gave us a factor of 125 improvement in effective simulation speed, as well as the crucial decrease in statistical error. Whilst it may be possible to vary other parameters, or perhaps even several parameters at once, this will induce a greater divergence between propagation in the reference and target systems, and so affect the extent of any improvement.

C. Nonlinearity
As described in the main text, the key feature of the degenerate-pump spontaneous FWM interaction used in our simulations is that it produces pairs of entangled signal and idler fields (photons). It shares this feature with the more commonly used second order parametric nonlinearity [40] (usually denoted χ (2) ), which generates the same kind of entangled pairs but from single pump photons, not pairs. Although this difference is not trivial, the SDE equations for a χ (2) nonlinearity are similar in form, although with no quantum noise term applied to the pump field. To test what differences might appear between the two models, some simulations were also done with this χ (2) model and system parameters rescaled to match the nonlinear effects; the results were remarkably similar in character, indicating that our conclusions are not specific to the FWM model we present here.
Note that the nonlinear response is treated here as if it were instantaneous. This is a reasonable approximation since typical nonlinear response times in dielectrics are of the order of femtoseconds or less [43].

D. Evolution: from density matrix to SDEs
Since to our knowledge the positive-P Fokker-Planck equation for the degenerate-pump FWM interaction Hamiltonian are not currently in the literature, so we now sumarize the derivation. This interaction Hamiltonian (1) gives a contribution to the density matrix evolutioṅ The expansion of the density matrix in terms of coherent states, albeit for just a single mode, is ρ = |α α † α|α † P(α, α † ;t) d 2 α d 2 α † .
Using the standard positive-P operator correspondences, we can convert the density matrix dynamics into a dynamics for