Nonlinear atomic tunnelling boosted by bright squeezed vacuum

Experimental details

Both coherent and quantum light sources are pumped with the same femtosecond laser pulse (790 nm, 28 fs, 10 kHz) generated by a Ti:sapphire multipass amplifier laser system (Femtolaser). The coherent light source centred at 1,580 nm and with a pulse duration of 70 fs, is produced by a commercial optical parametric amplifier (Light Conversion, TOPAS-Prime). For the quantum BSV light source, the pump beam collimated to a 4-mm diameter is propagated through two cascaded 3-mm BBO crystals. Both BBO crystals are cut for type-I collinear frequency-degenerate phase matching to generate high-gain parametric down-conversion. Here, the optical axes are oriented oppositely in the horizontal plane to minimize the spatial walk-off. The distance between two crystals is set at 80 cm so that only the spatial mode with the lowest diffraction undergoes the phase-sensitive amplification. The pulse duration of the BSV light is measured to be approximately 150 fs using the technique of cross-correlation frequency-resolved optical gating based on sum-frequency generation of the to-be-calibrated BSV pulse and a reference infrared pulse. The effect of different pulse durations for coherent and BSV lights is discussed in the section ‘Effect of pulse duration’.

The second-order correlation function g⁽²⁾ of the generated BSV is measured using the standard Hanbury Brown–Twiss technique. The BSV source is split at a nonpolarizing 50/50 beam splitter into two channels, the photon number distributions of which are measured with two fast InGaAs photodiodes. The signals from the photodiodes are recorded by a multichannel oscilloscope and then integrated over time to present the total photon numbers. Finally, the second-order correlation function is calculated by \(g^(2)=\langle n_1n_2\rangle /(\langle n_1\rangle \langle n_2\rangle )\), where n₁ and n₂ are photon numbers measured for two channels. Considering the typical measured g⁽²⁾ value of 1.5, there are multiple frequency modes amplified in the BSV generation process, which has been further confirmed by a spectral filtering of the BSV light using a standard 4-f monochromator (see section ‘Mode analysis of photon number statistics’ for details). It is worth noting that to preserve the peak intensity required for tunnelling ionization, no spectral filtering is applied in our photoionization experiments. The details of controlling the correlation function g⁽²⁾ and validating the multi-mode BSV light are demonstrated in the section ‘Mode analysis of photon number statistics’.

The beam diameters for both the coherent light and the BSV light are approximately 3 mm and are tightly focused by a concave silver mirror (focal length f = 75 mm) inside an ultrahigh-vacuum chamber onto the sodium vapour jet. The sodium vapour is produced in a resistively heated crucible at 170 ^°C and collimated by a 2-mm skimmer. The peak intensity of the coherent light in the interaction region is estimated to be 1 × 10¹³ W cm⁻². The corresponding Keldysh parameter is calculated to be 1.48, which places the ionization process in the non-adiabatic tunnelling regime. The photoelectron energy spectrum (Fig. 2b) featuring a smooth structure and the absence of resonant peaks⁶¹ also indicates that ionization proceeds through a tunnelling pathway, without significant population of intermediate states such as Na(4s). This distinguishes our results from regimes in which atomic saturation and resonances suppress quantum-optical enhancement effects⁶² and allows for an observation of the BSV-induced enhancement.

The strong-field tunnelling ionization of a single atom generates both photoelectron and photoion, whose momenta are coincidently measured using the cold-target recoil ion momentum spectroscopy reaction microscope^52,53. These charged particles are accelerated by a homogeneous electric field (7 V cm⁻¹) with the assistance of a weak magnetic field (11 G) and finally strike the detectors at the opposite ends of the spectrometer. The momenta of the emitted electron p^(ele) and the corresponding parent ion p^(ion) for each laser–atom interaction event are reconstructed from the measured time of flight and position of impact during the offline analysis. The principle of momentum conservation in the centre-of-mass frame dictates that the momentum sum of the electron and ion from the same interacting atom should vanish, that is, p^(ele) + p^(ion) = 0. In practice, the coincidence condition of \(| p_z^(\rme\rml\rme)+p_z^(\rmi\rmo\rmn)| < 0.2\,\rma.u.\) along the time-of-flight axis of the spectrometer, which has the highest momentum resolution, is used to select the right coincidence events from other background signals or false coincidences.

Both coherent and quantum light sources are converted from linear to elliptical polarization by a broadband quarter-wave plate. By rotating the fast axis of the plate relative to the initial linear polarization, the relative intensity between the two orthogonal components is controlled while introducing a fixed π/2 phase shift, thereby producing the desired ellipticity of 0.7, which is appropriate for performing angular streaking analysis using the accessible BSV light (see section ‘Choice of elliptical polarization’ for details).

QADK theory

The semiclassical ADK theory^57,58,59, commonly used in strong-field ionization, treats the electron quantum mechanically while describing the light as a classical field. It represents the adiabatic limit of the strong-field approximation^12,63,64 and the Perelomov–Popov–Terent’ev theory^65,66,67,68. In this framework, the doubly differential ionization rate is given by

where Coulomb-related prefactors are not included here for consistency. For circular polarization, the time dependence factorizes out, allowing us to directly replace A(t) with A₀ and F(t) with F₀. For elliptical polarization, we replace A(t) with εA₀ and F(t) with F₀ for ionization along the major axis or replace A(t) with A₀ and F(t) with εF₀ for that along the minor axis. For simplicity of the expression, we adopt the substitution A(t) → A₀ and F(t) → F₀, and the singly differential rate is expressed as

This expression is often further expanded for small initial transverse momentum, leading to an alternative form

The quantum enhancement observed in our experiment is fundamentally rooted in the quantum statistical property of the BSV light, which is imprinted onto the emitted electron by light–electron entanglement inherent to the quantum-optical description of the tunnelling ionization process. To rigorously capture this physics, we have developed a QADK theory based on the entangled light–electron Hamiltonian, given by

where \(\hat\bfA=\fracg\omega (\hata\boldsymbol\varepsilon +\hata^\dagger \boldsymbol\varepsilon ^* )\) is the vector potential operator, ε is the polarization vector and g is the coupling strength. The electron is initially in the ground state |g⟩, and the light is initially in a BSV state \(\hatS|0\rangle \), with \(\hatS=\exp \left(\frac12\xi ^* \hata^2-\frac12\xi \hata^\dagger 2\right)\) the squeezing operator, where the squeezing parameter ξ = re^iϕ, with r and ϕ the squeezing amplitude and squeezing phase, respectively. The initial light–electron state is, therefore, \(|\Psi _\rangle =|g\rangle \otimes \hatS(\xi )|0\rangle \), which will evolve into an entangled light–electron state as their interaction initiates.

To facilitate calculations of entangled light–electron dynamics, we apply a unitary transformation \(\widetilde\hatO=\hatS^\dagger \hatU_\rmF^\dagger \hatO\hatU_\rmF\hatS\) with \(\hatU_\rmF=\exp (-\rmi\hatH_\rmFt)\), which turns the Hamiltonian into

where the transformed vector potential \(\hat\bfA_\rmS=\hatS^\dagger \hatU_\rmF^\dagger \hat\bfA\hatU_\rmF\hatS\), and the initial state becomes \(|\widetilde\varPsi _\rangle =|g\rangle \otimes |0\rangle \), so that the light part of the initial state is simply vacuum in the transformed gauge.

For a BSV state that has a large squeezing amplitude r, \(\hat\bfA_\rmS\) can be approximated by

where \(\hatp_\Lambda =(\hata-\hata^\dagger )/\sqrt2\rmi\) is the momentum quadrature of BSV light, \(\boldsymbol\varepsilon (t)=\frac\rmi2(\rme^-\rmi\omega t\boldsymbol\varepsilon -\rme^\rmi\omega t\boldsymbol\varepsilon ^* )\) with ε = e_x + iεe_y, \(A_\equiv 2\fracg\omega \sqrt\fracn1+\varepsilon ^2\) with n the average photon number is the vector potential amplitude defined in such a way that it matches its classical counterpart with the same actual intensity because \(I_\rmB\rmS\rmV=\fracc\varepsilon _\omega ^2\overline\langle 0=\fracc\varepsilon _2\omega ^2(1+\varepsilon ^2)A_^2=I_\rmc\rmo\rmh\) with ε₀ the vacuum permittivity and c the vacuum light speed, and \(\widehatA_\rme\rmff\equiv \sqrt2A_\widehatp_\Lambda \) is the effective vector potential under BSV quantum light. Correspondingly, the effective electric field \(\widehatF_\rme\rmff\equiv \sqrt2F_\widehatp_\Lambda \). We note that the stretched quadrature has been rotated to the momentum direction in equation (6). In this transformed gauge, the classical and quantum aspects of the light–electron interaction factorize, which greatly simplifies the calculation of the entangled dynamics. In particular, the quantum-optical description proceeds in close analogy to its classical counterpart: we replace the classical vector potential A₀ with the effective vector potential \(\widehatA_\rme\rmff\) and the classical electric field F₀ with \(\widehatF_\rmeff\), whereas all other elements of the theory remain unchanged. It is important to note, however, that the light momentum \(\widehatp_\Lambda \) enters the effective fields \(\widehatA_\rmeff\) and \(\widehatF_\rmeff\), and the distribution of p_Λ encodes the information of the quantum statistical property of the quantum light. For BSV, whose initial state is vacuum |0⟩ in the transformed gauge, it can be expanded in the light momentum basis:

Therefore, the quantum statistical weight of light momentum p_Λ corresponds to \(\chi (p_\Lambda )=\rme^-p_\Lambda ^2\) in the normalized probability distribution. The initial vacuum state |0⟩ is expanded in the momentum quadrature basis p_Λ because \(\hat\bfA_\rmS\) is diagonal in this representation. This simplification relies on the fact that, under large squeezing, the contribution from the conjugate position quadrature becomes negligible. This property allows us to reduce the dimensionality of the expansion. By contrast, other studies expand the initial state in a coherent-state basis^2,40 or a number-state basis⁶⁹. Although all expansions are, in principle, equivalent, it is important to note that the coherent-state basis is not orthogonal; therefore, we cannot ignore off-diagonal elements and expect to obtain identical results.

In this framework, the doubly differential ionization rate under a single-mode BSV light is given by

where \(\chi (p_\Lambda )=\rme^-p_\Lambda ^2\) stems from the quantum statistical property of the BSV light. The QADK theory is a generalization of the semiclassical ADK theory to the quantum domain. It explicitly accounts for the entanglement between the photoelectron momentum p_e and the light momentum p_Λ, which imprints the quantum statistical property of BSV on the emitted electrons.

The derivation above for the single-mode BSV can be extended to the multi-mode case⁶⁰ by replacing the field momentum p_Λ by \(P_\Lambda =\sqrt\sum _jp_\Lambda j^2\). This substitution follows from the statistical independence of the Schmidt modes. Because the relative squeezing angles between different modes are random, the momenta add incoherently. This assumption is implicit in the multi-mode nonlinear response derived in refs. ^3,39, in which all cross-terms between modes vanish. This treatment is physically justified, as it also ensures gauge invariance of the resulting observables. Meanwhile, the Gaussian statistical distribution for the field momentum should be accordingly changed to be \(\chi (P_\Lambda )=\frac1\varGamma (N/2)P_\Lambda ^N-1\rme^-P_\Lambda ^2\) with N the mode number of the field derived using convolution across different field modes. This leads to the doubly differential ionization rate under an N-mode BSV light

The quantum statistics of the BSV light redistributes the entire weight of the QADK differential rate, which directly manifests as the extended photoelectron kinetic energy spectrum we observe, as visually confirmed in Extended Data Fig. 1. Extended Data Fig. 1a shows the correlated light–electron energy distribution obtained from the QADK theory, in which light energy \(E_\Lambda \equiv A_^2P_\Lambda ^2\) and electron energy \(E_\rme=p_\rme^2/2\), whereas Extended Data Fig. 1c corresponds to the projection onto the electron energy axis E_e. The semiclassical ADK ionization rate represents a single slice of the QADK distribution along a fixed light momentum at \(P_\Lambda =1/\sqrt2\) (or \(E_\Lambda =A_^2/2\)):

As the classical slice locates below the dominate region of the positively tilted quantum distribution, two key consequences emerge: (1) the BSV light yields a much higher ionization rate than the coherent light; and (2) the BSV light produces a much higher and broader photoelectron energy distribution.

To quantify the quantum enhancement, we attenuate the BSV field until its peak photoelectron momentum matches that produced by the coherent light, as shown in Extended Data Fig. 1b,d. This condition is met when the average photon number of the BSV light is reduced by a factor of 14 relative to its original value, corresponding to a quantum enhancement factor of 14 based on peak photoelectron momentum matching. The residual discrepancy from the observed 20-fold quantum enhancement might be attributed to quantum correlations across different field modes, which are not included in our current model.

Analysis of electron number statistics

Following the development of the QADK theory above, the number of electron n_e can be calculated as

where K is a scaling parameter depending on the density of the atom and conversion efficiency, γ(P_Λ) is the singly differential ionization rate for a specific field momentum P_Λ:

According to the probability conservation relationship

the distribution of electron number can be calculated as

The resulting fit to the BSV-driven electron number distribution is presented in Fig. 2a (orange solid line). It is necessary to note that the logic to calculate the electron number distribution is the same as in refs. ^3,39, except that we operate in the field momentum space.

Linear scaling between effective intensity and g
⁽²⁾

The second-order correlation function g⁽²⁾ is a fundamental quantity characterizing the photon number statistics and quantum properties of light fields. We note that although the rate of a perturbative n-photon process generally scales with the nth order correlation function g⁽ⁿ⁾ (refs. ^{70,71,72,73,74,75,76,77,78,79}), the characterization of the non-classical photon number statistics of the BSV light itself is primarily captured by the standard and widely adopted metric of second-order correlation function g⁽²⁾. This is because g⁽²⁾ serves as the primary measure of non-classical statistical properties of the quantum light, directly reflecting the super-Poissonian character inherent to BSV.

For the multi-mode BSV light generated in the present study, the total photon number is distributed across multiple frequency modes⁶⁰. Understanding how g⁽²⁾ scales with the effective intensity (photon number per mode) is essential for designing quantum light sources with tailored statistical properties. The multi-mode BSV state is generated by applying a multi-mode squeezing operator to the vacuum state, assuming the modes are independent:

where N is the number of modes⁶⁰ and \(\hatS_k(\xi _k)=\exp \left(\frac12\xi _k^* \hata_k^2-\frac12\xi _k\hata_k^\dagger 2\right)\) is the squeezing operator for the kth mode, with \(\xi _k=r_k\rme^\rmi\phi _k\). The total photon number operator is \(\hatn=\sum _k=1^N\hata_k^\dagger \hata_k\). The zero-delay second-order correlation function is defined as

where \(\langle \widehatn\rangle \) and \(\langle \widehatn^2\rangle \) are the first and second moments of the total photon number distribution.

For a single-mode BSV, the moments are

For N independent modes, the total moments are

Assuming all modes are equally squeezed (r_k = r), the moments simplify to

Substituting into g⁽²⁾, we obtain

in the limit of large photon number \(\langle \widehatn\rangle \).

The above relationship, g⁽²⁾ = 1 + 2/N, links the second-order correlation function to the effective number of independently amplified modes N (refs. ^80,81). It is consistent with known properties of BSV, in which the single-mode limit yields g⁽²⁾ = 3, whereas the expression converges to the coherent state value g⁽²⁾ = 1 as N → ∞. Moreover, it demonstrates that the measured g⁽²⁾ value of the generated BSV light inherently reflects an effective number of independent, equally squeezed modes. In the case of unequal gain for different modes, it only breaks down to a refinement in the interpretation of N from an exact mode count to the Schmidt number⁶⁰, without altering the fundamental conclusions of our analysis. Therefore, we may simply refer to N as the mode number.

For a multi-mode BSV state, the photoelectron momentum spectra are generated through interactions with individual modes and summed across all modes. Consequently, the peak momentum distribution directly reflect the effective intensity per mode, which is determined by the average photon number \(\langle \widehatn_k\rangle \) in each mode:

The correlation function g⁽²⁾ can then be directly linked to \(\langle \widehatn_k\rangle \) as

In this expression, the average photon number \(\langle \widehatn_k\rangle \) per mode is proportional to the effective intensity I_eff, whereas the total photon number \(\langle \widehatn\rangle \) is proportional to the average pulse energy P. The effective intensity scales with the mean photon number per mode, \(\langle \widehatn_k\rangle \), reflecting the fact that the total ionization signal results from the contributions of interactions with each individual field mode. Therefore, a linear scaling between the effective intensity I_eff and the second-order correlation function g⁽²⁾ results:

Effect of pulse duration

In our experiments, although the pulse durations of the 70-fs coherent pulse and 150-fs BSV pulse are different, they have minor influence on our comparison for the two key reasons.

First, in our measurement, the comparison between coherent and BSV light sources focuses on the effective peak intensity. Tunnelling ionization, as a highly nonlinear process, is predominantly driven by the peak electric field of the pulse. Our angular streaking technique directly measures this effective peak intensity, enabling a controlled comparison at the intensity level relevant to the process.

Second, to confirm that the observed photoelectron kinetic energy spectrum is insensitive to pulse duration, we performed a simulation using coherent pulses of varying durations. As shown in Extended Data Fig. 2, the photoelectron kinetic energy spectra from tunnelling ionization driven by 70 fs and 150 fs coherent pulses are virtually identical. This demonstrates that over this duration range, the temporal envelope does not introduce qualitatively different physics that could account for the observed quantum enhancement.

Therefore, although identical pulse durations for coherent and BSV pulses are, in principle, ideal for comparison, the evidence supports that our comparison is physically sound, and the conclusions about quantum enhancement remain robust.

Choice of elliptical polarization

The choice of elliptical polarization for the driving field, as opposed to linear, was a deliberate and important aspect of our experimental design. In linear polarization, the photoelectron momentum distribution is affected by rescattering and interference effects, resulting in a complex structure that cannot be uniquely mapped to a single value of the vector potential amplitude. The elliptical streaking field, by contrast, enables angular streaking, which generates a photoelectron momentum distribution that is centred around an ellipse corresponding to the vector potential amplitude.

The use of elliptical polarization was further necessitated by signal-to-noise considerations. Although circular polarization provides an ideal rotating field for angular streaking, it distributes pulse energy equally between two orthogonal axes, effectively halving the peak intensity along any direction. Given the experimental challenge of generating intense BSV, circular polarization would have resulted in ionization yields too low for statistically robust coincidence measurements. Elliptical polarization preserves a significant rotating field component for angular streaking while maintaining a strong field amplitude along the major axis, thereby ensuring sufficient ionization rates without compromising the physical interpretation of the streaking signal. The selected ellipticity thus represents an optimal trade-off between streaking interpretation and ionization efficiency, in which the peak momentum along the direction of maximum emission probability sits at εA₀.

Calibration of effective intensity

In both the QADK and semiclassical ADK theories, electron momentum, rather than kinetic energy, is directly involved as shown in equations (8) and (10). Therefore, in our study, we calibrate the intensity using photoelectron momentum distributions rather than energy. As shown in Extended Data Fig. 3, these distributions maintain a Gaussian-like profile across all g⁽²⁾ values, enabling a reliable and accurate extraction of the peak momentum. This photoelectron peak momentum directly reflects the effective intensity of the BSV field. Here, the enhanced tails of the momentum distribution are a demonstration of the quantum statistical nature of the BSV light.

Nevertheless, in Fig. 3a of the main text, we present the results in terms of kinetic energy in eV rather than momentum, because the energy unit of eV is a more intuitive and comprehensible unit, leading to an easier comparison with existing literature⁴⁰. It should be noted that the calibrated value of the effective intensity differs slightly depending on whether the matching is performed using kinetic energy or momentum, as \(\langle p_\rme^2\rangle \ne \langle p_\rme\rangle ^2\). That said, the key conclusions of our work, particularly the quantum enhancement and linear scaling of the effective intensity with g⁽²⁾, remain invariant when we carry out the matching using the same physical observable, be it momentum or kinetic energy.

Mode analysis of photon number statistics

Through singular value decomposition, a multi-mode field can be expressed on the basis of independent effective modes. The effective number of these modes, generally quantified by the Schmidt number⁶⁰, is a measure of the number of independent, uncorrelated components within the BSV field, as widely used in both theory and experiments. In the present study, the generated BSV light contains a single spatial mode and multiple frequency modes. To confirm this, we performed spectral filtering of the BSV light using a standard 4-f monochromator. The slit in the monochromator can change the spectral bandwidth and thus alter the frequency mode number of the BSV source. Extended Data Fig. 4a demonstrates the measured second-order correlation function g⁽²⁾ of the BSV light after the 4-f monochromator as a function of the slit width in the 4-f monochromator, which manipulate the spectral bandwidth and the frequency mode of the BSV light. After spectral filtering, we achieved g⁽²⁾ values up to 2.6, which confirms the contribution of multiple frequency modes to our BSV source. We note that, to maintain the high peak intensity required for atomic ionization in our experiment, we did not use spectral filtering to the generated BSV that enters the vacuum chamber for the ionization of the atoms as a multi-mode light.

In our experiments, we primarily control the g⁽²⁾ parameter of the BSV source by varying the pump power incident on the nonlinear crystal. Extended Data Fig. 4b shows the measured dependence of g⁽²⁾ on pump power, demonstrating that we can systematically tune this parameter. To validate the independent, equally squeezed multi-mode characterization of the BSV light, we measured photon number distributions at distinct g⁽²⁾ values. Extended Data Fig. 5 presents six typical photon number distributions (blue shades), in which lower g⁽²⁾ values exhibit multi-mode characteristics: the peak is shifted towards higher photon numbers, and the overall distribution shows a progressively closer resemblance to a Poissonian distribution. For quantitative comparison, we calculated the expected photon number distributions, shown in Extended Data Fig. 5 (blue curves), according to the multi-mode photon number distribution³:

where Γ is the gamma function, N is the number of modes, and n denotes the photon number. The effective mode number is estimated as \(N=2/[g^(2)-1]\) (equation (20)). The good agreement between measurement and calculation provides strong evidence for the multi-mode characterization of our BSV source and the validity of equation (20) derived under the assumption of multiple independent equally squeezed modes. This treatment is a common practice for characterizing sources based on parametric down-conversion and is consistent with the methodology used in established works^1,3,60.