Subcycle squeezing of light from a time flow perspective
Matthias Kizmann Thiago Lucena de M. Guedes Denis V. Seletskiy Andrey S. Moskalenko Alfred Leitenstorfer Guido Burkard





Department of Physics and Center for Applied Photonics,
University of Konstanz, D-78457 Konstanz, Germany

Department of Physics and Center for Applied Photonics,
University of Konstanz, D-78457 Konstanz, Germany

Department of Engineering Physics, Polytechnique Montréal, H3T 1J4, Canada
Department of Physics and Center for Applied Photonics,
University of Konstanz, D-78457 Konstanz, Germany
 
 andrey.moskalenko@uni-konstanz.de
 Department of Physics and Center for Applied Photonics,
University of Konstanz, D-78457 Konstanz, Germany
 
 Department of Physics and Center for Applied Photonics,
University of Konstanz, D-78457 Konstanz, Germany
  
  guido.burkard@uni-konstanz.de
  Department of Physics and Center for Applied Photonics,
University of Konstanz, D-78457 Konstanz, Germany


Subcycle squeezing of light from a time flow perspective
    Guido Burkard
     December 30, 2023
========================================================



Light as a carrier of information and energy plays a fundamental role in both general relativity and quantum physics, linking these areas that are still not fully compliant with each other. Its quantum nature and spatio-temporal structure are exploited in many intriguing applications ranging from novel spectroscopy methods of complex many-body phenomena <cit.>  to quantum information processing <cit.> and subwavelength lithography <cit.>. Recent access to subcycle quantum features of electromagnetic radiation <cit.> promises a new class of time-dependent quantum states of light. Paralleled with the developments in attosecond science <cit.>, these advances motivate an urgent need for a theoretical framework that treats arbitrary wave packets of quantum light intrinsically in the time domain.
Here, we formulate a consistent time domain theory of the generation and sampling of few-cycle and subcycle pulsed squeezed states, allowing for a relativistic interpretation in terms of induced changes in the local flow of time. Our theory enables the use of such states as a resource for novel ultrafast applications in quantum optics and quantum information. 

 







In conventional optics, there was no control over the absolute phase of a light pulse with respect to its envelope. This situation changed abruptly with the advent of femtosecond frequency combs<cit.>. The carrier-envelope phase is especially relevant for few-cycle pulses which are at the heart of attosecond technology and extremely nonlinear optics<cit.>. So far, this area has exploited coherent states of light which come closest to the classical picture of an electromagnetic field with well-defined amplitude and phase. Only recently, quantum optics was carried to an analogous level with subcycle analysis of the noise properties of infrared electric fields by electro-optic sampling with femtosecond laser pulses <cit.>. Here, the concept of a carrier-envelope phase loses its meaning because highly nonclassical states of light may exhibit excessive phase fluctuations or no well-defined phase at all. Instead, it is the relative timing of a quantum noise pattern with respect to a subcycle probe which gains relevance. 
Interesting and unexpected phenomena may arise when the physics of such synchronal states of light is explored on a subcycle scale.


One of
the most fundamental nonclassical states of light is the squeezed vacuum state <cit.>. Employed in an interferometer, it can provide a detector sensitivity below the shot-noise limit <cit.> that was used in the measurements leading to the recent experimental discovery of gravitational waves<cit.>. 

Stationary squeezed vacuum states  carry a persisting flux of photons <cit.>, whereas there is a finite number of photons contained in a pulsed state. Pulsed squeezed light <cit.> with many optical oscillation cycles was realized based on parametric down-conversion <cit.>. The case of the generation of pulsed broadband ultrashort squeezing<cit.> was analyzed in terms of a large set of single-frequency or broadband shaped temporal modes <cit.>, where each mode is characterized
separately.


 The alternative time-domain
approach based on electro-optic sampling provided first insights into the temporally
resolved dynamics of few-cycle squeezing supported by a simplified theoretical description <cit.>.

This theory was restricted to the low squeezing regime or selected points in time. Furthermore, it has not yet incorporated a finite probe pulse duration for the detection of the temporal quantum noise patterns. Here, we present methods that overcome these limitations by considering extremely short (in terms of optical cycles) pulses of squeezed vacuum light and provide a general theoretical picture of their generation and subcycle-resolved detection. We predict that the expected asymmetry between the anti-squeezed and squeezed temporal noise, well known in conventional quantum optics, can in fact be reversed for ultrabroadband driving fields which are also within reach of state-of-the-art experiments.  Our results shed light on the interplay between the quantum nature of electromagnetic fields and general relativity, leading to creation of particles out of vacuum on ultrashort time scales.



We consider the setup depicted schematically in Fig. <ref>a. A coherent mid-infrared electric driving field with a classical amplitude ℰ(z,t) propagates along the z-axis and is polarized along the x-axis. It enters a thin transparent optical crystal, which has a length l, centered at z=0, and possesses a non-vanishing second-order nonlinearity. This leads to an interaction of the driving field with the co-propagating quantum electric field component ε̂(z,t), which is polarized along the y-axis, belongs to the same frequency range, and corresponds to the vacuum field when entering the crystal (see Methods for
details). The nonlinear interaction in the crystal subsequently modifies ε̂(z,t).  Far from material resonances the second-order susceptibility χ^(2) can be treated as a dispersionless tensor. Adopting the propagation of plane waves with perfect phase matching inside the crystal and a broadband version of the slowly varying amplitude approximation the evolution of the field operator ε̂(z,t) is given by (see Methods)

    ∂ε̂/∂ z   =d/nc_0[∂ℰ/∂ tε̂+(ℰ-n^2/d)∂ε̂/∂ t].

Here the effective nonlinear coefficient d is determined by the properties of χ^(2), n is the linear refractive index and c_0 is the vacuum speed of light. Equation (<ref>) can be solved analytically by using the method of characteristics, where the initial conditions of a partial differential equation are propagated along   characteristic curves  (see Methods for
details). In our case the characteristic curves are given by fixed values of an auxiliary time variable τ(z,t), which at any fixed z monotonically increases with t and may be called conformal time, in analogy to cosmology <cit.>. The gravitational analogues of the characteristic curves are world lines <cit.>. Denoting the quantum field at the entrance of the generation crystal as ε̂_in(t)≡ε̂(z=-l/2,t), the solution of equation (<ref>) can be found as

    ε̂(z,t)=∂τ(z,t)/∂ tε̂_in(τ(z,t)).

The form of the resulting quantum field
can be attributed to two different effects. Firstly, due to the driving field ℰ(z,t) the wave propagation of the incoming quantum field ε̂_in inside the crystal is governed by the conformal time τ(z,t). This effect is illustrated by the world lines in the horizontal plane of Fig. <ref>b.
The change of their density upon the propagation through the crystal determines the resulting difference in the flow of the conformal time with respect to the lab time t. Secondly, there is a modulation of the amplitude of the quantum field by the inverse conformal factor <cit.> ∂τ(z,t)/∂ t. This factor is shown in the vertical plane of Fig. <ref>b at the crystal exit and is essentially related to the behaviour of ∂ℰ(z,t)/∂ t (cf. ref. Riek2017), as discussed below. Note that the quantum field is enhanced (suppressed) in the same time segments where the density of the world lines is increased (decreased), see vertical lines in Fig. <ref>b. In a simplified picture, neglecting the effect of the modified density of the world lines, the amplitude modulation of the quantum field leads to a corresponding modulation in the temporal profile of its quantum noise. The quantum fluctuations of ε̂ are suppressed beneath the vacuum level in certain time segments while exceeding it in the neighbouring segments. Squeezing (anti-squeezing) can be connected to the deceleration (acceleration) of the local flow of the conformal time τ(z,t), determined by its derivative with respect to the lab time t. For a part of the discussion it is convenient to use a retarded reference frame with t'=t-n/c_0z, z'=z and τ'(z',t')=τ(z,t)+nl/2c_0. First we discuss the essence of the effect and its detection for the case of an idealized half-cycle driving pulse; afterwards we study the single-cycle case. Figure <ref>a depicts the evolution of the conformal time in the retarded reference frame through the crystal resulting from a half-cycle driving field of the form (see Supplementary Information)

    ℰ'(t')   =ℰ_0(Γ_0t'),

where ℰ_0 is
the amplitude of the field and Γ_0 determines its duration.

The squeezing strength r=|ℰ_0d|Γ_0 l/(nc_0) is a dimensionless positive parameter. At the entrance of the crystal the conformal time matches the lab time: τ(z=-l/2,t)=t and τ'(z'=-l/2,t')=t'. With the interaction turned on, while the quantum field propagates through the crystal its conformal time starts to deviate from the lab time. 

 Via Eq. (<ref>) the conformal time directly links the quantum field at the exit of the crystal at time t to the quantum field at the entrance at time τ_out(t)=τ(z=l/2,t), see Fig. <ref>b. According to Eq. (<ref>) the modification of the variance of the quantum electric field after passing the crystal is given by  (dτ_out/d t)^2 if considered at a fixed time moment. Experimentally,  the temporal dynamics of the quantum fluctuations of the electric field in dependence on the relative time delay of the probe pulse t_d can be traced utilizing the electro-optic sampling technique and statistical readout <cit.>. The time resolution and sensitivity are limited by the duration of the probe pulse t_p during which the information about the quantum field is collected. This is a significant aspect because apart from the change in its amplitude the quantum field of the incoming vacuum is effectively subjected to a modified flow of time at the crystal exit, as can be seen from Eq. (<ref>). Alternatively we can say that in the reference frame using the conformal time τ_out(t) in place of the lab time t the probe pulse shape and duration are modified according to the flow of τ_out(t) [cf. Eq. (<ref>) in Methods]. In order to obtain the time-resolved detected variance V(t_d)=⟨[ε̂^(d)_out(t_d)]^2⟩ for the outgoing quantum electric field ε̂_out(t)=ε̂(z=l/2,t), i.e. the temporal trace of its quantum fluctuations, we employ the quantum theory of electro-optic sampling <cit.> (see Methods for details). We define the relative detected variance (RDV) as

    RDV(t_d)   =V(t_d)
    -(Δε̂_vac)^2/(Δε̂_vac)^2.

Here (Δε̂_vac)^2=⟨[ε̂^(d)_in]^2⟩ denotes the detected variance of the incoming quantum vacuum field, which does not depend on the delay time t_d.



Figure <ref>b depicts the RDV for the case of the half-cycle pulse Eq. (<ref>), a very short detection time t_p=0.49 fs and different squeezing strengths r.  
For r=0.1,  RDV(t_d)
is an almost perfectly odd function (red line) with respect to the center of the driving pulse at t_d=0 fs. In this case, the temporal trace nearly coincides  with the waveform proportional to the third derivative of the driving field d^3 ℰ'(t_d)/d t_d^3 (dashed black line). The latter corresponds to the exact analytic solution in the limit of small r and t_p, which can be obtained via the squeezing operator in the frequency domain (see Supplementary Information). This dynamics differs from the simplified picture where for each fixed time moment the variance of the quantum field is determined by  [dτ_out/d t]^2 and in the low squeezing regime (r≪1) leads to temporal profiles following the first derivative dℰ'(t_d)/d t_d <cit.>. Notice that the deviation becomes less significant if multi-cycle, more narrowband driving fields are considered.


For small r, the RDV scales linearly with r and therefore remains symmetric. A build-up of asymmetry between the time segments with reduced and excess quantum noise can be clearly observed when r is increased. Firstly, we can see that the magnitude of the detected quantum fluctuations becomes more pronounced in the anti-squeezing period with respect to the squeezing period. This can be attributed to the fact that the conformal time τ_out(t) must always increase monotonically with t. As the slope of τ_out(t) approaches zero, the flow of time comes to a halt and according to  Eq. (<ref>) ε̂_out(t) must vanish. In the simplified picture, looking only at the prefactor in Eq. (<ref>), this would lead to an almost complete elimination of the quantum noise and squeezing limited by 100%. In contrast, acceleration in the flow of time in principle may lead to arbitrarily large prefactors determining the magnitude of ε̂ in Eq. (<ref>)  and therefore to an arbitrarily strong enhancement of the quantum noise. In the full picture including detection, the increase caused by the prefactor in Eq. (<ref>) is partly counteracted by the slower (quicker) flow of  the conformal time  in ε̂_in(τ_out(t)) in the case of squeezing (anti-squeezing). The duration of the probe pulse seen by the incoming vacuum field thus effectively  becomes smaller (larger) for squeezing (anti-squeezing) and the detected variance V(t_d) increases (decreases) due to this effect <cit.>. Equivalently, the effect of the lower (higher) local world line densities (see Fig. <ref>b) alone would enhance (suppress) the detected quantum noise.  Still, after the unperturbed vacuum contribution (Δε̂_vac)^2 is subtracted in Eq. (<ref>), the asymmetry in the RDV(t_d) is preserved to a large extent for sufficiently short detection times t_p. Secondly, with increasing r the time segments of squeezing become broader while the time segments of anti-squeezing narrow down.
This additional asymmetry can be comprehended looking at Fig. <ref>b. It is clear that the time intervals with steep slopes in the flow of the conformal time τ_out(t) (anti-squeezing) take less space on the horizontal axis, representing the lab time t, than the intervals with flat slopes (squeezing). This is what we observe for the case of r=2 in Fig. 3b (blue line). Uniquely to the strong squeezing in the time domain, polarity reversal of the driving field (Fig. 3a) does not change the position of the squeezing valley while reverses the arrival time of the anti-squeezing burst.


In Fig. <ref>d, the dynamics of the RDV is shown for different probe pulse durations t_p and a fixed squeezing strength r=2. As t_p grows (cf. Fig. <ref>c), higher frequencies in ε̂_out(t) are not captured anymore leading to a flattening of the detected temporal traces. The narrow anti-squeezing peaks are affected more strongly by such changes than the valley regions of squeezing. When the probe pulse duration approaches the oscillation time of the driving field, the asymmetry is almost completely lost.

We now switch our attention to the case of a single-cycle driving pulse (Fig. <ref>a).

Figure <ref>b shows the resulting dynamics of the RDV, also in terms of the degree of squeezing (see Supplementary Information). For the given amplitude of the driving field, the maximum positive degree of squeezing amounts to 49.6% while the minimum negative degree of squeezing (corresponding to anti-squeezing) constitutes only -35%. Therefore, the previously described asymmetry in favour of higher absolute values for anti-squeezing is reversed in this case. At first sight, one might wonder whether Heisenberg's uncertainty principle is violated in the described situation. Through the inspection of the RDV trace for the polarity-switched drive (inset in Fig. <ref>c), we discover that the amplitude of anti-squeezing is in fact nearly twice that of the squeezing (Fig. <ref>b). It is only in this temporally non-local sense that we are able to restore the familiar dependence of quantum noise on the squeezing strength (Fig. <ref>c and Supplementary Information), and hence directly recover Heisenberg's uncertainty  principle.


Our work provides a consistent time-domain theory of the generation and sub-cycle-resolved detection of ultrabroadband waveforms of pulsed squeezed radiation, linking these processes directly to a change in the local flow of time induced by the coherent driving field. This constitutes a non-perturbative analytical solution for the operator of the generated quantum field and the electro-optically sampled traces of its quantum noise. This solution is valid at any time and also in the high-driving regime, restricted only by conditions underlying the broadband version of the slowly varying amplitude approximation and causality. We applied our theory to predict time traces of the detected variance and corresponding degrees of squeezing belonging to ultrashort quantum fields created by half-cycle and single-cycle driving.
The results for the detected variance show that the use of electro-optic sampling for a time-resolved measurement of the noise pattern of a squeezed field inherently introduces an admixture of vacuum fluctuations. However, the asymmetries between squeezing and anti-squeezing are preserved by this measurement technique for sufficiently short probe pulses.
 Finally, we theoretically predict an effect that at first glance looks paradoxical: the conventionally observed asymmetry between squeezing and anti-squeezing can be reversed for specially designed driving fields.




§ METHODS



 §.§ Time-domain equation for the generation process

We consider a process where in a χ^(2) nonlinear crystal  a coherent driving field ℰ(z,t) polarized along the x-axis propagates in the z-direction (cf. Fig. <ref>a). ℰ(z,t) interacts with the y-component of a co-propagating quantum field ϵ̂(z,t). This gives rise to the y-component of the nonlinear polarization

 P̂^(2)(z,t)=-ϵ_0 dℰ(z,t)ϵ̂(z,t)

but generates no x-component. Such a situation can be realized for various nonlinear crystals, e.g., for ZnTe in a configuration where the present unit vectors e_x and e_z  would be aligned along the [001] and [110] crystallographic directions, respectively <cit.>. The driving field can be external but the situation is similar if this field is generated in the same crystal <cit.>.

We write the corresponding inhomogeneous wave equation for ϵ̂(z,t) with P̂^(2)(z,t) as the source term. Then we decompose ϵ̂(z,t)=∫_-∞^∞dΩ ε̂(z;Ω)e^i(k_Ω z-Ω t) and in the same way ℰ(z,t) and P̂^(2)(z,t) in plane waves with an almost ideally linear dispersion k_Ω=n/c_0Ω.
 With the latter, for a thin enough crystal perfect phase matching is a good approximation.
In the resulting frequency-domain integro-differential equation <cit.>
we make an assumption that
the term ∂^2 ε̂(z;Ω)/∂ z^2 can be neglected in comparison with k_Ω∂ε̂(z;Ω)/∂ z. This leads to

    ∂ε̂(z;Ω)/∂ z=-idΩ/n c_0∫^∞_-∞dΩ_1 ℰ^*(z;Ω_1-Ω)ε̂(z;Ω_1).

Transforming Eq. (<ref>) back into the time-domain we then obtain Eq. (<ref>).


 In essence, this step is similar to the well-known slowly varying amplitude approximation  (SVAA) <cit.>. Note, however, that in our case no amplitude or envelope was introduced in the time domain. Therefore, also few-cycle and subcycle fields with no unambiguously defined carrier frequency and phase can be described by this approach. It  may be considered as a broadband version of the SVAA<cit.>. An important prerequisite is the assumed linear dispersion holding in the frequency range of interest.
In a dispersive medium ultrashort pulses would experience shape distortion while propagating through the crystal <cit.>, requiring a more involved description of the generation process.




 §.§ Method of characteristics

The method of characteristics is used to solve a given partial differential equation (PDE) by reducing it to a set of ordinary differential equations (ODEs), which define parametric curves  in the space of the PDE variables (so-called characteristic curves). These curves allow then to construct the solution of the PDE.
Let us assume that the original variables of the PDE, in our case Eq. (<ref>), depend on the parameter ζ. Comparison of the total derivative of ϵ̂ with respect to ζ with Eq. (<ref>) gives us ODEs for z(ζ) and t(ζ) (with their respective initial conditions):

    dz/dζ=1,       z(ζ=0)=-l/2;
    
       dt/dζ=-d/nc[ℰ(z,t)-n^2/d]_z=z(ζ),t=t(ζ),
         t(ζ=0)=τ.

Here an additional parameter τ was introduced as the initial condition for t(ζ) when ζ=0. Solving Eq. (<ref>) reveals that ζ=z+l/2. Changing the variable from ζ to z in Eq. (<ref>) and solving it, we obtain t(z,τ), where now also an explicit dependence on τ appears. Finally, the function t(z,τ) can be inverted to obtain τ(z,t) with τ(z=-l/2,t)=t.
The functional form of τ(z,t) ensures that the total derivative of τ(z(ζ),t(ζ)) with respect to ζ vanishes,

    dτ(z(ζ),t(ζ))/dζ   =∂τ(z,t)/∂ z-d/nc_0[ℰ(z,t)-n^2/d]∂τ(z,t)/∂ t=0,

meaning that τ(z(ζ),t(ζ)) remains constant when moving along a characteristic curve {z(ζ),t(ζ)} parameterized by ζ.
This PDE for τ(z,t) can now be used to prove that Eq. (<ref>) indeed satisfies Eq. (<ref>). Therefore, the only remaining nontrivial task is to solve Eq. (<ref>) finding τ(z,t). Inserting the solution for τ(z,t) into Eq. (<ref>) then gives ε̂(z,t).

 For the purpose of solving Eq. (<ref>) the retarded reference frame with t'=t-n/c_0z, z'=z and ℰ(z,t)=ℰ'(t') is utilized. This transforms Eq. (<ref>) into a separable ODE. τ(z,t) can then be expressed as

    τ(z,t)   =f^-1(f(t-n/c_0z)+1/nc_0(z+l/2))-nl/2c_0

in the original reference frame, where f(t') is defined as solution of the following ODE

    df(t')/dt'=1/dℰ'(t').

Note that here no initial condition needs to be stated because it can be chosen arbitrarily. The necessary initial condition is already incorporated in the solution for τ(z,t) with τ(z=-l/2,t)=t, regardless of the specific form of f(t').




 §.§ Quantum electro-optic sampling

The quantum theory of electro-optic sampling <cit.> is employed in order to model the detection process. The variance of the outgoing quantum field ε̂'_out(t') is calculated as a function of the relative delay time t_d with respect to the center of an ultrashort higher-frequency probe pulse ℰ'_𝗉(t') (for convenience the retarded reference frame is used here). The resulting operator of the sampled field is given by

    ε̂^(d)_out(t_d)=∫_-∞^∞dt'  R(t_d-t')ε̂'_out(t').

The detector function R(t_d) is determined by the normalized intensity of the probe field
that for a fast oscillating probe field  can be replaced by the normalised intensity envelope:

    R(t_d)=|ℰ'_𝗉(t_d)|^2/∫_-∞^∞dt' |ℰ'_𝗉(t')|^2≈I'_𝗉(t_d)/∫_-∞^∞dt'  I'_𝗉(t').

Then for the outgoing quantum field, determined by Eq. (<ref>) at z=l/2, the operator of the sampled field can be expressed via the incoming quantum field ε̂_in(t) in the original reference frame as

    ε̂^(d)_out(t_d)=∫_-∞^∞dt  R(t_d-τ^-1_out(t))ε̂_in(t),

where x=τ^-1_out(y) denotes the inverse function of y=τ_out(x).
Let us consider the case where the incoming field corresponds to the vacuum state.
Then Eq. (<ref>) shows that sampling of the generated squeezed field by the probe pulse in an inertial reference frame
can be alternatively viewed as sampling of the bare vacuum field in a reference frame with a non-inertial time flow given by τ^-1_out(t). In the latter reference frame the shape of the used probe pulse is then effectively transformed.
In order to characterize the dynamics of the resulting quantum fluctuations we evaluate the relative detected variance (RDV) that is given by Eq. (<ref>). For the corresponding calculations in this paper we used
 a Gaussian shape of the probe intensity envelope with the FWHM duration t_p so that R(t_d)=2√(ln2)/√(π)t_pexp(-4ln2 t_d^2/t_p^2).  
 
The calculated values of the RDV typically are small compared to the values of RDV and respective degrees of squeezing that would follow directly from the simplified picture (see Fig. <ref>b).
The reason for this is that the RDV [cf. Eq. (<ref>)] is defined as the relative difference between the detected variance of the generated field and that of the unperturbed vacuum fluctuations.
In order to resolve the signal in the time-domain the probe duration has to be much shorter than the typical time scales of the signal. Thus the frequency spectrum of the probe has to be much broader than that of the coherent field driving the squeezing. This means that the probe pulse will not only sample modes at the frequencies that are affected by the driving field but also modes at higher frequencies which remain untouched. Therefore, the contribution of the unperturbed vacuum fluctuations to the RDV has to be significantly larger than the contribution from the frequency range affected by the generation process in order to appropriately resolve the temporal squeezing pattern.
 This can be interpreted as an admixture of vacuum fluctuations to the signal that can be modelled by losses. However, it is important to note that rather than being losses that can be avoided, these losses are inherent to the measurement process and are necessary in order to appropriately resolve the signal in time. Nevertheless,  losses do not necessarily represent an obstacle for the reconstruction of the degree of the generated squeezing <cit.>, which is shown in Fig. <ref>c.

Note that the relevant quantity measured experimentally, and connected to the calculated RDV, is the relative differential noise (RDN) <cit.>. It is given by


    RDN(t_d)   =(√(V(t_d))-Δε̂_vac)Δε̂_vac/(Δℰ_p,SN)^2 ,

[see Eq. (<ref>) for comparison].
Here (Δℰ_p,SN)^2 denotes the shot noise (SN) of the probe field ℰ_p.






§ ACKNOWLEDGEMENTS

We thank P. Sulzer and R. Haussmann for discussions. Support
by DFG via SFB767, by Baden-Württemberg Stiftung via
the Eliteprogramme for Postdocs
(project “Fundamental aspects of relativity and causality in time-resolved quantum optics ”)
and by Young Scholar Fund of the University of Konstanz is acknowledged.
M.K. is indebted to the LGFG PhD fellowship program of the University of Konstanz.



§ AUTHOR CONTRIBUTIONS

A.S.M., D.V.S. and G.B. conceived the idea. A.S.M. managed the project and supervised the research. M.K. found the exact analytical solution in the time domain, performed numerical calculations and prepared the figures. T.L.M.G. obtained the perturbative analytic solution via the squeezing operator in the frequency domain. M.K., T.L.M.G. and A.S.M.  wrote the first version of the paper.  D.V.S. and A.L. provided several important physical insights and interpretations. All authors discussed the results and contributed to the writing of the final manuscript.









§ SUPPLEMENTARY INFORMATION




 §.§ Slowly varying amplitude approximation and causality

Equation (<ref>) shows that the quantum field at the exit of the crystal at time t is influenced by the quantum field at the entrance at time τ_out(t). This fact raises the question whether in our calculations this can technically lead to an information exchange faster than
 the speed of light or even the conformal time preceding the lab time. Already the first case would contradict causality. Figure <ref>b shows τ_out(t)  as a function of t (full black line). Without the driving field, τ_out(t) is delayed with respect to t by the time nl/c_0 which the light needs to travel through the crystal (see dashed black line). The red-coloured area in the figure marks the region for τ_out(t) where the delay between τ_out(t) and t is less than l/c_0, which is the time required for light to travel the distance l in vacuum. Due to causality, the information about the incoming quantum field cannot propagate to the end of the crystal in a time shorter than l/c_0. Causality forbids this area for τ_out(t). Since its deviation from the lab time t grows with higher squeezing strengths r, the restriction that τ_out(t)-t+nl/c_0 (vertical difference between the full red and dashed black line in Fig. <ref>b) should not exceed (n-1)l/c_0 ultimately puts an upper bound to the maximally allowed r. However, the use of a broadband version of the SVAA, underlying our theoretical approach, imposes an even stricter upper bound to the squeezing strength that then already ensures causality, as we will show in the following. 
Within the utilized broadband SVAA we neglected ∂^2 ε̂(z;Ω)/∂ z^2 with respect to k_Ω∂ε̂(z;Ω)/∂ z. This is possible only if ε̂(z;Ω) does not change  too rapidly with z. The corresponding condition limiting the strength of the coherent field, which induces variations of ε̂(z;Ω), can be found from Eq. (<ref>) differentiating it with respect to z. Assuming that the driving field has a peak amplitude ℰ_0, we obtain ℰ_0 d≪ n^2. With the frequencies of the ultrabroadband driving field being on the order of Γ_0 this
restriction can be reformulated as a condition on the squeezing strength r,

    r≪nΓ_0 l	/c_0.



Solutions obtained within the broadband SVAA do not automatically preserve causality, i.e. the impossibility of information exchange faster than the speed of light. The condition 

    g(r):= Γ_0 max_t'[τ'(z'=l/2,t')-t')]≤Γ_0(n-1)l/c_0,

formulated in the retarded reference frame, ensures causality.
This also effectively sets a boundary for r. For r≪ 1, the left hand side of the expression above is proportional to r with a proportionality factor given by max ℰ'(t')/ℰ_0≲1. For higher r, the increase becomes even slower [see Fig. <ref>].

Since for nonlinear crystals typically n and n-1 are on the same order,
the limit imposed by the range of validity of the broadband SVAA, Eq. (<ref>), already ensures causality.
This means that it is sufficient to respect Eq. (<ref>). For squeezing strengths beyond the SVAA the changes in the refractive index induced by the driving field would be comparable to the initial refractive index n of the crystal. Furthermore, effects of higher order nonlinearities like χ^(3) would need to be taken into account.

Considering an experimentally realistic situation with Γ_0/(2π)=26 THz and a l=15 μm thick ZnTe crystal with n=2.57, Eq. (<ref>) leads to an estimation

    r ≪ 20.

In comparison, the limit given by causality in this case can be calculated to

    g(r)≤ 13.

In this work the refractive index of the nonlinear crystal is treated as a dispersionless constant. In doing so we neglect the phonon resonances, located for ZnTe at around 5 THz, and assume that the main contributions to the signal will be given by higher frequencies.
In principle, from Eq. (<ref>) we see that higher values of r become possible if the crystal thickness is increased. This possibility is, however, limited by the requirement to fulfill the ultrabroadband phase matching conditions. Overcoming this issue would open a way for ultrabroadband squeezing in the high gain regime.



 §.§ Half-cycle pulse

The shape of the half-cycle pulse ℰ'(t')=ℰ_0(Γ_0 t') used as the driving field in Fig. <ref>
is depicted in Fig. <ref> as the black line. Because of its non-vanishing integral over time, these type of pulses are
actually not supported in the far-field region with respect to the light source. However, for our particular problem the selected  temporal  profile can serve as a good model of a waveform capable to propagate in the far-field zone, which is relevant to the experimental situation. Let us consider the following allowed pulse shape (red line in  Fig. <ref>)

    ℰ'_r(t')=ℰ_0(Γ_0 t')-0.1ℰ_0(Γ_0 t'/10).

The correction introduced by the second term forces the integral over time to vanish.

In principle, one might want to use such or similar shapes for our calculations directly in place of a simpler profile ℰ'(t'). However, the latter has the advantage that the solution for the evolution of the conformal time inside the crystal could be found in an analytical form also for stronger driving fields. Compared with ℰ'(t'), the field ℰ'_r(t') has a somewhat lower central positive peak and possesses negative side wings, which are  low-amplitude and extended in time . An important point is that the difference between both discussed waveforms remains small for any moment in time.
Furthermore, in the calculations of the detected variance the difference is even less than it could be anticipated looking at Fig. <ref>.
Indeed, let us consider for example the case of the vanishing probe pulse duration t_p and weak driving when the dynamics of the detected variance is governed by d^3 ℰ'(t_d)/d t_d^3 or d^3 ℰ'_r(t_d)/d t_d^3, respectively. In the third derivative the contribution of the second term on the right hand side of Eq. (<ref>) is suppressed by an additional factor of 10^3 with respect to the first term. Thus the difference in the calculated variances may be neglected with a high degree of accuracy.



 §.§ Fitting procedure for the extraction of the degrees of squeezing

The red line in Fig. <ref>b depicts the squeezing pattern in the simplified picture.  Here the degree of squeezing is determined by a calculation similar to Eq. (<ref>), with a minus prefactor (following the convention to assign positive values of the degree of squeezing for squeezing and negative ones for anti-squeezing) and using the ratio between the variances of the outgoing field and vacuum  field taken at the relevant moment. At times where the driving field vanishes this ratio depends exponentially on the squeezing strength r.
For such a point inside a squeezing (anti-squeezing) region, it exponentially decays (grows) with r. In this case it is guaranteed that the squeezing turns into  anti-squeezing and vice versa by changing the polarity of the driving field. This exponential dependence is shown in Fig. <ref>c as the red lines, where the degree of squeezing at t_d=0 fs is plotted for both polarities of the driving field amplitude in dependence of the squeezing strength r.
The observed asymmetry between the maximal positive values of the degree of squeezing (corresponding to squeezing) and its minimal negative values (corresponding to anti-squeezing) can serve as a direct measure of the degree of squeezing.

Figure <ref>c demonstrates that a similar measure can be used to extract the degree of squeezing from the RDV within the full description accounting for a finite detection window. The two grey lines reflect the corresponding maximal squeezing and anti-squeezing values obtained for the RDV for different polarities of the driving field amplitude. In order to extract the degree of squeezing from these lines an exponential fit

    h_±(r)=A_1[exp(± A_2r)-1]

is applied. The degree of squeezing is then given by -RDV/A_1. Since the two fitting functions only differ in the sign of the exponent, the parameters A_1 and A_2 can be fitted well either to the squeezing or the anti-squeezing values, where the respective differences decrease for shorter probe pulse durations. One should select rather the squeezing branch, since the corresponding dynamics is better resolved than for the anti-squeezing, due to the slower flow of the conformal time in the former case. The result of this fit is shown in Fig. <ref>c as the dashed red lines. For the case of a much shorter probe pulse duration the extracted squeezing and anti-squeezing curves would coincide with the red curves of the simplified picture.

The right axes of Figs. <ref>b and <ref>c show the degrees of squeezing obtained in this way from the RDV and within the simplified picture. 






 §.§ Single-cycle RDV for different squeezing strengths

Fig. <ref> depicts the RDV for the single-cycle driving field for different squeezing strengths. For r=0.5 the maximal absolute degree of squeezing exceeds the maximal absolute degree of anti-squeezing. Increasing the squeezing strength however restores the conventionally expected asymmetry of more pronounced anti-squeezing than squeezing.






 §.§ Squeezing operator approach

It can be shown<cit.> that up to the first order in the driving coherent field
Eq. (<ref>) results in the following transformation of the annihilation operator for the quantum field ε̂:

    â(z;Ω)=  â(Ω)+(z/l+1/2)sign (Ω)∫^∞_-∞dΩ_1 Ξ_sym(Ω,Ω_1)â^†(Ω_1),

where Ξ_sym(Ω,Ω_1)=-idl/n c_0sign(ΩΩ_1)√(|ΩΩ_1|)ℰ_Ω_1+Ω, ℰ_Ω denotes the Fourier transform of the driving field in the retarded reference frame ℰ'(t'), z∈[-l/2,l/2] and â(Ω)≡â(z=-l/2;Ω) ≡â_in(Ω) corresponds to the incoming field ε̂_in. For the outgoing field ε̂_out at the end of the crystal we have â(z=l/2;Ω) ≡â_out(Ω). The transformation of the creation operator is given by the Hermitian conjugate of (<ref>).

The linear relations between â_out(Ω), â^†_out(Ω) and â_in(Ω), â^†_in(Ω) can be viewed as a generalized  Bogolyubov transformation <cit.>. It can be induced by a unitary squeezing operator Ŝ as â_out(Ω)=
Ŝâ_in(Ω)Ŝ^† [â^†_out(Ω)=
Ŝâ^†_in(Ω)Ŝ^†]. One finds<cit.>

    Ŝ=exp{1/2∬^∞_0dΩ_1dΩ_2[Ξ^*_sym(Ω_1,Ω_2)â_(Ω_1)â_(Ω_2)-Ξ_sym(Ω_1,Ω_2)â^†_(Ω_1)â^†_(Ω_2)
    
    +Ξ^*_sym(-Ω_1,Ω_2)â^†_(Ω_1)â_(Ω_2)-Ξ_sym(-Ω_1,Ω_2)â_(Ω_1)â^†_(Ω_2)]}.

This continuous-mode form of the operator Ŝ extends <cit.> the definition of a discrete multimode squeezing operator <cit.>. Moreover, since in a thin crystal the generated photon pairs copropagate along the same axis, both frequency conversion and parametric down conversion <cit.> happen simultaneously.
This leads to the creation of photons in an even broader frequency range.
The resulting broadband quantum state of light in the Schrödinger picture is obtained by applying the squeezing operator (<ref>) to the vacuum state |0⟩: |{ξ_Ω}⟩=Ŝ|0⟩.


The outgoing quantum field in the retarded reference frame reads

    ε̂'_out(t') =i∫^∞_0dΩ √(ħΩ/4π nϵ_0 c_0 A)[â_out(Ω)e^-iΩ t' - â_out^†(Ω)e^iΩ t'],

where A is the normalization area and we have taken into account the relation k_Ω=Ω n/c_0. In this work we are interested in the difference in the detected variances for the outgoing and incoming fields ⟨ 0|[ε̂^(d)_out(t_d)]^2|0⟩
-⟨ 0|[ε̂^(d)_in]^2|0⟩, which occurs in Eq. (<ref>).
For the incoming vacuum field the detected variance does not depend on the delay time t_d of the probe field. In the limit case of the vanishing probe duration t_p, the considered difference converges to the normally ordered <cit.> variance of the outgoing quantum field
⟨ 0|[ε̂_out(t_d)]^2| 0⟩ at the time t=t_d. Using  Eq. (<ref>) and Eq. (<ref>) we can compute this quantity up to the first order in the driving field. In principle, it is possible to obtain also the higher order terms <cit.>. For that, the further expansion terms should be calculated in Eq. (<ref>) from the non-perturbative result for the operator Ŝ provided by Eq. (<ref>).

Let us consider a driving field with a half-cycle shape (<ref>). Its Fourier transform reads ℰ_Ω=ℰ_0/2Γ_0sech(πΩ/2Γ_0). Then up to the first order in perturbation theory we can calculate

    ⟨ 0|[ε̂_out(t_d)]^2| 0⟩ = ħΓ^2_0/4π nϵ_0 c_0 Ar/6[tanh^3(Γ_0 t_d)sech(Γ_0 t_d)-5tanh(Γ_0 t_d)sech^3(Γ_0 t_d) ].

This result holds when r=|ℰ_0d|Γ_0 l/(nc_0) ≪ 1. It can be seen that under these conditions the variance is an odd function of time (cf. dashed black line in Fig. <ref>b) and scales linearly with the squeezing strength r. Thus the temporal profiles of the squeezing and anti-squeezing remain symmetric with respect to each other until the limit of validity of the utilized approximation is reached.












