interferometer which introduces a controllable delay $\Delta\tau$ (phase shift
$\varphi$) between the Er and Ein fields. An important requirement for the
applicability of the method is that the photon number of the local oscillator
field should be much larger than that of Ein to be characterized (for details
see Ref. [141]). The fields after the BS are detected by a balanced
differential photodetection system consisting of two identical photodiodes
(PD), which provides at each value of $\varphi$ the photocurrent difference
$i_{\varphi}$. Setting the delay stage around $\Delta\tau\approx 0$, the
characterization of the quantum state of light can be achieved by recording
for each shot the value of $i_{\varphi}$ as a function of $\varphi$. The
beauty of the homodyne detection scheme [140] is that it provides access to
the measurement of the values of the field quadrature. The values of the
photocurrent difference $i_{\varphi}$ are directly proportional to the
measurement of the electric field operator i.e.,

$\displaystyle\hat{E}_{in}(\varphi)\propto\hat{x}_{\varphi}=\cos(\varphi)\hat{x}+\sin(\varphi)\hat{p}.$
(93)

The unknown proportionality factor between the measured photocurrent and the
quadrature values in Eqs. (93) can be obtained by using the condition that the
variance of the vacuum state is $1/2$, as for the vacuum state
$\left\langle\hat{x}_{\varphi}\right\rangle=0$, and $(\Delta
x_{\varphi})^{2}=1/2$. In other words, the homodyne data are scaled according
to the measured vacuum state quadrature noise. Experimentally, that is all
that someone has to do in order to characterize the light state. This is
because repeated measurements of $\hat{x}_{\varphi}$ at each $\varphi$
provides the probability distribution
$P_{\varphi}(x_{\varphi})=\langle{x_{\varphi}}|{\hat{\rho}}|{x_{\varphi}}\rangle$
of its eigenvalues $x_{\varphi}$, where
$\hat{\rho}\equiv\outerproduct{\phi}{\phi}$ is the density operator of the
light state to be characterized and $|{x_{\varphi}}\rangle$ the eigenstate
with eigenvalue $x_{\varphi}$. The density matrix $\hat{\rho}$, which provides
complete information about the light state, can be obtained in the Fock basis
by calculating the matrix elements $\rho_{nm}$ using an iterative
Maximum–Likelihood procedure beautifully described in Ref. [142]. Having these
values, the mean photon number of the light state can be obtained by the
diagonal elements $\rho_{nn}$ of the density matrix $\hat{\rho}$, and the
relation $\langle{n}\rangle=\sum n\rho_{nn}$.

One of the most complete and commonly used ways to visualize the quantum
character of light states is via the quasi–probability distribution in phase
space $(x,p)$, namely, the Wigner function $W(x,p)$, which has been
extensively addressed over the years by an enormous number of research
articles and books (see for instance refs. [140, 143, 144, 145, 146]). What is
important to be mentioned here is that: a) the marginal distributions
$P(x)=\langle{x}|{\hat{\rho}}|{x}\rangle=\int^{+\infty}_{-\infty}\differential
p\,W(x,p)$ and
$P(p)=\langle{p}|{\hat{\rho}}|{p}\rangle=\int^{+\infty}_{-\infty}\differential
x\,W(x,p)$ yield the position and momentum distributions, respectively, and b)
if we introduce a phase shift $\varphi$, the $x$ and $p$ components rotate by
$\varphi$ in phase–space via the operator ${\hat{U}}(\varphi)$ shown in Eq.
(86). Taking this into account, the equation which connects the
$P_{\varphi}(x_{\varphi})$ with $W(x,p)$ for each $\varphi$ is

$\displaystyle
P_{\varphi}(x_{\varphi})=\langle{x}|{\hat{U}}(\varphi){\hat{\rho}}{\hat{U}^{\dagger}}(\varphi)|{x}\rangle$
(94) $\displaystyle=\int^{+\infty}_{-\infty}\differential
p_{\varphi}W(x_{\varphi}\cos{\varphi}-p_{\varphi}\sin{\varphi},x_{\varphi}\sin{\varphi}+p_{\varphi}\cos{\varphi}),$

where $x_{\varphi}=x\cos{\varphi}+p\sin{\varphi}$ and
$p_{\varphi}=-x\sin{\varphi}+p\cos{\varphi}$. The above equation is called
Radon transformation and has been extensively studied in the theory of
tomographic imaging [147], while the inverse Radon transformation, implemented
via the standard filtered back-projection algorithm, is the key tool in
quantum state tomography and quantum state reconstruction [128, 145, and
references therein]. Hence, the Wigner function can be reconstructed by
applying the algorithm directly to the quadrature values $x_{\varphi,k}$,
using the formula
$W(x,p)\simeq\frac{1}{2\pi^{2}N}\sum_{k=1}^{N}K(x\cdot\cos(\varphi_{k})+p\cdot\sin(\varphi_{k})-x_{\varphi,k})$
[128], where $N$ is the number of the recorded quadrature–phase pairs
$(x_{\varphi,k},\varphi_{k})$. In the previous expression, $k$ is the index of
each value and
$K(z)=\frac{1}{2}\int_{-\infty}^{\infty}\differential\xi\,|\xi|\exp(i\xi z)$
is called integration kernel with
$z=x\cdot\cos(\varphi_{k})+p\cdot\sin(\varphi_{k})-x_{\varphi,k}$. The
numerical implementation of the integration kernel requires the replacement of
the infinite integration limits with a finite cutoff $k_{c}$. In order to
reduce the numerical artifacts (rapid oscillations) and allow the details of
the Wigner function to be resolved, the value of $k_{c}$ is typically set to
be around $k_{c}\sim\\!4$. Alternatively, the $W(x,p)$ can be also obtained
using the density matrix formalism [142, 148].

In the following, we provide three characteristic examples of the homodyne
traces, and their corresponding Wigner functions as have been calculated for
an optical coherent, “cat”, and “kitten-like” states. These cases have been
experimentally measured in [91, 93] by means of the arrangement shown in Fig.
9. In Fig. 15 we show the calculated homodyne trace (left panel) and the
Wigner function (right panel) of a coherent state. This corresponds to the
case where the HHG process is switched off.

Figure 15: Quantum tomography of an optical coherent state $\ket{\alpha}$. In
the left panel we show the calculated homodyne detection signal $x_{\phi}$,
while the corresponding Wigner function centered at $|\alpha|=2$ is shown in
the right panel.

The calculation has been performed using a coherent state $\ket{\alpha}$ of
amplitude $|\alpha|=2$. In this case,
$\hat{\rho}\equiv\outerproduct{\alpha}{\alpha}$ and therefore the mean value
$\left\langle\hat{x}_{\varphi}\right\rangle=\langle{\alpha}|\hat{x}_{\varphi}|{\alpha}\rangle$
exhibits an oscillation with
$\left\langle\hat{x}_{\varphi}\right\rangle=|\alpha|\cos{\varphi}$.
Furthermore, the $P_{\varphi}$ defined in Eq. (94) and the quantum fluctuation
(which is the same as the vacuum state $(\Delta x_{\varphi})^{2}=1/2$) remains
constant along the cycle. All these characteristic features are shown by the
Wigner function which depicts a Gaussian distribution given by

$\displaystyle W(\beta)=\frac{2}{\pi}e^{-2|\beta-\alpha|^{2}}.$ (95)

In the above equation we have used the transformation $\it
x\equiv\real[\beta-\alpha]$ and $\it p\equiv\imaginary[\beta-\alpha]$ (where
$\beta$ is a variable) which centers the Wigner function at the origin as for
$\beta=\alpha$, $x=0$ and $p=0$. This is a convenient way for plotting a
Wigner function, as well as writing its expression in a more compact form.

The situation changes dramatically when the HHG process and the QS are
switched on. In Fig. 16 we show the calculated homodyne traces (left panels)
and the corresponding Wigner functions (right panels) of an optical “cat-like”
(Figs. 16 (a) and (b)) and a “kitten-like” (Figs. 16 (c) and (d)) state
created when we introduce the conditioning on HHG.

Figure 16: Quantum tomography of optical “cat-like” CSS ((a) and (b)) and
“kitten-like” CSS ((c) and (d)) states created when the HHG and the QS are
switched on. The left panels show the calculated homodyne detection signal
$x_{\phi}$ and the right panels the corresponding Wigner functions centered at
$|\alpha|$. The calculations have been preformed using the state
$\ket{\phi_{c}}=\ket{\alpha+\chi}-\xi_{(IR)}\ket{\alpha}$ (with $|\alpha|=2$)
created by the conditioning on HHG process. For the optical “cat-like” state
((a) and (b)) the $|\chi|$ was set at 0.8, while for the “kitten-like” state
((c) and (d)) is reduced to 0.1.

In Figs. 16 (a) and (b), the calculation has been performed using the optical
“cat” state $\ket{\phi_{c}}=\ket{\alpha+\chi}-\xi_{(IR)}\ket{\alpha}$ shown in
Eq. (81) with $|\alpha|=2$, and $|\chi|=0.8$. In Figs. 16 (c) and (d), the
value of $|\chi|$ has been reduced to 0.1. In both cases,
$\hat{\rho}=\outerproduct{{\phi_{c}}}{{\phi_{c}}}$ and the mean value
$\left\langle\hat{x}_{\varphi}\right\rangle=\bra{\phi_{c}}\hat{x}_{\varphi}\ket{\phi_{c}}$
exhibits an oscillation with
$\left\langle\hat{x}_{\varphi}\right\rangle=\sqrt{\langle{n}\rangle}\cos{\varphi}$.
On the other hand, and contrarily to what we saw in the case of a coherent
state, $P_{\varphi}$ changes along the cycle and depicts a minimum around the
center of the quantum fluctuation which is larger than the vacuum state. All
these characteristic features are shown in the equation of Wigner function
which reads,

$\displaystyle W(\beta)=$ $\displaystyle\dfrac{2}{\pi
N}\Big{[}e^{-2\lvert\beta-\alpha-\chi\rvert^{2}}+e^{-\lvert\chi\rvert^{2}}e^{-2\lvert\beta-\alpha\rvert^{2}}$
(96)
$\displaystyle-\big{(}e^{2(\beta-\alpha)\chi^{*}}+e^{2(\beta-\alpha)^{*}\chi}\big{)}e^{-\lvert\chi\rvert^{2}}e^{-2\lvert\beta-\alpha\rvert^{2}}\Big{]},$

where $N=1-e^{-\lvert\chi\rvert^{2}}$ is the normalization factor for Eq.
(81). The Wigner function depicts the features of a genuine optical “cat”
(“kitten” for $\chi=0.1$) state as it has a ring-like shape with a minimum
around the center, and shows pronounced non-classical signatures by means of
its negative values.

####  V.2.5 Analysis of an experimentally reconstructed Wigner function

Obtaining quantitative information about the non-classicality of the optical
coherent state superposition measured by a quantum tomography method, requires
the optimization $k_{c}$ used to reconstruct the Wigner function and an error
analysis. Both are associated with the number of points of the measured
homodyne detection signal $x_{\phi}$, and the error introduced by the balanced
detector of the homodyne detection system in the measurement of the photon
number.

A convenient way to estimate the error is to use the Wigner function of a
coherent state [93]. In this case the error can be obtained by comparing
(subtracting) the ideal Wigner function of a coherent state with the Wigner
function of the state reconstructed from the experimental data. Figure 17 (a)
shows how this procedure can be used to find the optimum value of $k_{c}$ and
minimize the error of the amplitude of the reconstructed Wigner function
$W(x,p)$. The same procedure can be used for obtaining the error in the mean
photon number. This can be done using the density matrix elements $\rho_{nm}$
in the Fock basis $\ket{n}$, where the mean photon number is obtained from the
diagonal elements $\rho_{nn}$, and the relation $\langle{n}\rangle=\sum
n\rho_{nn}$. Figure 17 (b) shows the dependence of the accuracy of measuring
the photon number $\langle n\rangle$ on the mean photon number $\langle
n\rangle$ of the light state (Error (%) $=|\langle n_{rec}\rangle-\langle
n_{th}\rangle|/\langle n_{th}\rangle$). Where $\langle n_{rec}\rangle$ is the
mean photon number of a coherent state calculated using the number of data
points of the measured homodyne detection signal $x_{\phi}$, and $\langle
n_{th}\rangle$ is the mean photon number value of an ideal coherent state. As
shown in Fig. 17 (b), this procedure can been repeated for coherent states of
different photon number values. The dependence of the Wigner function
amplitude and mean photon number errors on the number of points of the
measured homodyne detection signal $x_{\phi}$ is shown in Figs. 17 (c) and (d)
for $k_{c}\approx 3.7$ and $\langle n\rangle=3$, respectively. A homodyne
trace which contains $\geq 5\times 10^{3}$ measured $x_{\phi}$ values (for
$0\leq\phi\leq\pi$) results in a Wigner function amplitude and mean photon
number error $\leq 1.5$ %.

Figure 17: (a) Dependence of the error amplitude of a reconstructed Winger
function of a coherent state on $k_{c}$, using a homodyne trace with $\sim
10^{4}$ points for $0\leq\phi\leq\pi$. (b) Dependence of the accuracy of
measuring the mean photon number on the mean photon number of the light state.
(c) Dependence of the Wigner function amplitude error on the number of points
of the measured homodyne detection signal $x_{\phi}$, for $k_{c}\approx 3.7$.
(d) Dependence of the accuracy of measuring the mean photon number of the
light state on the number of points of the measured homodyne detection signal
$x_{\phi}$, for $\langle n\rangle=3$. In all graphs, the red solid line is a
15-points running average of the data (gray points). Figs (a) and (b)
reproduced from Ref. [93].

After the error analysis, the non-classicality of the measured light state can
be obtained by quantitatively evaluating its quantum features (such as its
negativity) using the reconstructed Wigner function. The negativity is
associated with the contrast of the fringe pattern created by quantum
interference effects between the coherent states participating in the coherent
state superposition. One of the reasons why the negativity can be reduced
compared to the ideal case is the presence of decoherence effects on the
coherent state superposition after its creation. Although the optical coherent
state superpositions can be considered as (almost) decoherence free compared
to the particle state superposition, the presence of optical losses during the
propagation in optical elements affects the non-classicality of these light
states. This leads to a Wigner function with reduced negativity or reduced
fringe contrast. Also, the influence of the decoherence effects on the non-
classicality of the state depends on the difference between the coherent
states in the superposition $\chi$. Hence, large traditional cat-states
($\ket{\alpha}\pm\ket{-\alpha}$) are more sensitive to decoherence effects.
The decoherence can be quantitatively studied using a noise model that
introduces photon losses due to the interaction with a Gaussian reservoir [93,
149], via the von Neumann entropy $S$, or using the fidelity $F$ (see e.g.
Ref. [150]).

##  VI Discussion

The interaction of intense laser fields with matter leads to a highly non-
linear electronic response, and simultaneously influences the electromagnetic
field by generating radiation in the previously unpopulated field modes, and
leads to a backaction on the fundamental mode which drives the process.
Furthermore, the interaction induces correlations among the modes
participating in the process. With the full quantum electrodynamical
description of intense laser-matter interaction introduced in this tutorial,
we were able to explicitly show how the two most common processes of high
harmonic generation and above-threshold ionization, influence the quantum
state of the electromagnetic field modes. Utilizing conditioning procedures on
the different processes allowed us to perform quantum state engineering
protocols for generating non-classical optical field states with novel
properties, such as high photon number, extending towards extreme wavelength
regime and using ultrashort light pulses.

The fully quantized descriptions of HHG and ATI process allow to obtain
further insights into the dynamics of intense laser-matter interaction, from a
novel perspective by taking into account the quantum nature of the
electromagnetic field. This allows to answer the questions of the non-
classical properties of the electromagnetic field, and to investigate the
backaction of the interaction on the driving laser field in intense laser-
matter interaction. In particular, we uncovered dynamical changes of the
electromagnetic field which could not be addressed by the previous semi-
classical theory.

We introduced the quantum electrodynamical description of the process of HHG
and ATI, considering the natural continuum of electromagnetic field modes.
Considering a conditioning procedure on the emission of harmonic radiation
leads to entangled coherent states, and coherent state superposition with non-
classical signatures.

We obtained the quantum state of the EM field for the conditioning procedure
on ATI, which allowed to obtain the backaction on the fundamental mode. We
emphasize that the backaction on the EM field due to the electron dynamics in
the ATI process can be decomposed in 3 steps. First, the bound electron
wavefunction oscillates in the field, leading to a shift of the EM field modes
itself. In the second step, the electron transitions from the ground state to
the continuum, which influences the EM field state due to the coupling of the
dipole transition matrix element with the electric field operator. In the last
step, after ionization, the electron is driven by the field leading to a shift
of the EM field state due to the charge current associated with the electron’s
motion in the continuum.

With the ability to compute the EM field state conditioned on ATI, we obtain
field observables of the driving laser mode (photon number distribution and
its mean value) to understand the underlying dynamics, and the backaction
associated to the process. Interestingly, the direction of the electron
propagation (determined by the initial electron momentum after ionization)
determines if the shift of the field state amplitude is in or out of phase,
leading to an enhancement or depletion of the fundamental mode, respectively.

This quantum electrodynamical description, together with the conditioning
procedures introduced in theory and experiment, gives rise to the opportunity
of generating non-classical field states. The applicability and use of those
states towards a novel platform for quantum information science will be
outlined in the next section.

##  VII Outlook

This paper provides a tutorial to a new emerging interdisciplinary area of
science, in which ultrafast laser physics and attoscience merge with quantum
optics and quantum information science. In this sense it serves as an
invitation to researchers to join the efforts in this area, and for this
reason we devote special attention to the discussion of the outlook and
possible future investigations. We decided to organize this section using the
different levels: operational (concerning “day-to-day” generalization and
continuations of the research discussed in this tutorial), tactical
(connecting the present state-of-art status to future merge with quantum
technologies, focused on detection, characterization and applications of
massive superposition and entangled states generated with the present
framework), and strategic (relating to long term goals and aims - the
development of the new interdisciplinary platform for quantum technologies
that merges quantum information, attoscience and quantum optics).

###  VII.1 Outlook - Operational level

Typically, the operational plan describes the day to day running of the
operation, in our case the new area of science. The operational plan charts
out “small steps” or milestones to achieve the higher level (tactical) goals
within a realistic time-frame. The present paper sets the QED framework for
attoscience, and focuses on some applications related to conditional
measurements of HHG and ATI in atomic media that lead to the generation of the
massive Schrödinger cat states (massive superposition and entanglement
states). There is an obvious ocean of possible generalizations, and
applications of the present theory that indeed can be achieved on a “day-to-
day” basis. Here we list a few:

  * •

Even in the presently investigated schemes, there is plenty of place for
further studies, such as for instance the role of the CEP, and there is a
whole plethora of possibilities associated with using more complex
polarization states (elliptic, time-dependent elliptic, locally chiral, etc.),
combining polarization with two-colour fields and structuring polarization in
space (cf. [151, 152, 153, 154, 155, 156, 157, 158, 159]). This might lead to
creation of novel types of topological Schrödinger cats, very different from
those discussed so far in the literature (cf. [160]).

  * •

Another important consideration is the spin of the electron, which was
overlooked until recently [161, 162, 163, 164, 165], where it has been shown
that circular strong-field ionization produces spin-polarized electrons. As
such, this process may be exploited to produce entanglement between the spin
and the laser field. The electronic spin provides an additional degree of
freedom useful in heading towards delivering quantum technologies, as the
discrete spin can be more easily measured and manipulated than the continuous
electronic degrees of freedom typically considered.

  * •

The same QED theory, and conditioning methods, can be applied to multielectron
processes in atoms (ionization in general, NSDI by the interplay between
Electron Impact Ionization (EII) and Resonant Excitation with Subsequent
Ionization (RESI)). These can be at least analysed using the SFA approach.

  * •

The present framework can be extended to more complex targets such as
molecules, where e.g. chirality can lead to otherwise symmetry-forbidden
electronic responses, like electronic motion perpendicular to the polarization
plane [166], and electronic orbital angular momentum in ionization with
linearly polarized light [167], thus expanding the spectrum of possible
conditioning setups.

  * •

A completely new avenue is to apply the present framework to solid state and
condensed matter systems. In the first place, this could be done for weakly
correlated systems that can be described by semiconductor Bloch equations [46,
30], or a little more flexible Wannier-Bloch approach [47]. The next step
would be to apply the present framework to imaging of topological, chiral
and/or strongly correlated systems, such as topological insulators [31, 32] or
high $T_{c}$ superconductors (cf. [51]).

###  VII.2 Outlook - Tactical level

In order to seek for the creation of a new feasible platform for quantum
technologies, one needs clearly to relate “what we have” to the basic aspects
of quantum information science: detection, characterization, and finally,
applications of generated quantum correlations and quantum entanglement. At
this stage, experiments generate in a simplified view the massive
superposition of the fundamental laser mode [91], or massively entangled
states of several modes [92, 94].

The next directions of research should then focus on:

  * •

Detection and characterization methods of massive Schrödinger cat states. This
is a formidable challenge [168], which for large cats is extremely demanding
[130].

  * •

These kind of states are in the first place useful for quantum metrology in
general [169, 170, 171], and on photonic platforms [172, 173, 174] in
particular, and perhaps in quantum communications. It will thus be important
to design concrete protocols for their applications.

###  VII.3 Outlook - Strategic level

On the strategic level we aim at developing a novel quantum technology
platform combining attoscience, quantum optics and quantum information. Each
of these areas are already at a step where proof-of-principle experiments have
been developed. The point is to create the new platform that is totally and
intrinsically interdisciplinary. The research presented here, focused in
theory and experiments on conditioning, is only a part of this general vision
aimed at realising a universal and firmly established tool to offer completely
unknown solutions and developments, i.e. a set of stable and reproducible
methods to generate massive entangled states and massive quantum
superposition. Most, if not all of these methods, involve quantum
electrodynamics of intense laser–matter interactions, i.e. the framework
described in the present tutorial. Here is the more complete list of the
future research lines, related to these strategic objectives, aims and goals

  * •

Objective 1: Generation, characterization and concrete applications in QT of
entangled/quantum correlated states using conditioning methods.

  * •

Objective 2: Generation, characterization and concrete applications in QT of
entangled/quantum correlated states in strong–field physics and attosecond
science driven by quantum light; this includes driving by strongly squeezed
light, light with orbital angular momentum (cf. [175]), effects on
photoelectron orbital angular momentum [176], and more.

  * •

Objective 3: Generation, characterization and concrete applications in QT of
entangled/quantum correlated states in Zerfall processes, such as single [177,
178, 179] or double ionization (cf. [35]).

  * •

Objective 4: Studying quantitative and measurable effects of decoherence in
attoscience, by studying reduced density matrices of Zerfall–products in
ionization processes (cf. [180, 181, 182]).

The results of the present tutorial set up the basis for a road map toward
this novel platform of attoscience, quantum optics and quantum information
science for quantum technologies.

## Acknowledgements

We thank Emilio Pisanty for discussion regarding this work. ICFO group
acknowledges support from: ERC AdG NOQIA; Agencia Estatal de Investigación
(R$\&$D project CEX2019-000910-S, funded by MCIN/ AEI/10.13039/501100011033,
Plan National FIDEUA PID2019-106901GB-I00, FPI, QUANTERA MAQS
PCI2019-111828-2, Proyectos de I+D+I “Retos Colaboración” QUSPIN
RTC2019-007196-7); Fundació Cellex; Fundació Mir-Puig; Generalitat de
Catalunya through the European Social Fund FEDER and CERCA program (AGAUR
Grant No. 2017 SGR 134, QuantumCAT U16-011424, co-funded by ERDF Operational
Program of Catalonia 2014-2020); EU Horizon 2020 FET-OPEN OPTOlogic (Grant No
899794); National Science Centre, Poland (Symfonia Grant No.
2016/20/W/ST4/00314); European Union’s Horizon 2020 research and innovation
programme under the Marie-Skłodowska-Curie grant agreement No 101029393
(STREDCH) and No 847648 (“La Caixa” Junior Leaders fellowships ID100010434:
LCF/BQ/PI19/11690013, LCF/BQ/PI20/11760031, LCF/BQ/PR20/11770012,
LCF/BQ/PR21/11840013). P. Tzallas group at FORTH acknowledges LASERLABEUROPE V
(H2020-EU.1.4.1.2 grant no.871124), FORTH Synergy Grant AgiIDA (grand no.
00133), the H2020 framework program for research and innovation under the NEP-
Europe-Pilot project (no. 101007417). ELI-ALPS is supported by the European
Union and co-financed by the European Regional Development Fund (GINOP Grant
No. 2.3.6-15-2015-00001). P.S. acknowledges funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement No 847517. J.R-D. acknowledges support from the
Secretaria d’Universitats i Recerca del Departament d’Empresa i Coneixement de
la Generalitat de Catalunya, as well as the European Social Fund (L’FSE
inverteix en el teu futur)–FEDER. A.S.M. acknowledges funding support from the
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement, SSFI No. 887153. M.F.C. acknowledges
financial support from the Guangdong Province Science and Technology Major
Project (Future functional materials under extreme conditions -
2021B0301030005).

## References

  * Siegman [1986] A. Siegman,  _Lasers_ (University Science Books, California, 1986).
  * Glauber [1963a] R. J. Glauber, Coherent and incoherent states of the radiation field, Phys. Rev. 131, 2766 (1963a).
  * Glauber [1963b] R. J. Glauber, The quantum theory of optical coherence, Phys. Rev. 130, 2529 (1963b).
  * Sudarshan [1963] E. Sudarshan, Equivalence of semiclassical and quantum mechanical descriptions of statistical light beams, Phys. Rev. Lett. 10, 277 (1963).
  * Eberly and Allen [1987] J. H. Eberly and L. Allen,  _Optical Resonance and Two Level Atoms_ (Dover Publication, Unitated States of America, 1987).
  * Piszczatowski  _et al._ [2009] S. Piszczatowski, G. Łach, M. Przybytek, J. Komasa, K. Pachucki, and B. Jeziorski, Theoretical determination of the dissociation energy of molecular hydrogen, Journal of Chemical Theory and Computation 5, 3039 (2009).
  * Pachucki  _et al._ [1996] K. Pachucki, D. Leibfried, M. Weitz, A. Huber, W. König, and T. Hänsch, Theory of the energy levels and precise two-photon spectroscopy of atomic hydrogen and deuterium, Journal of Physics B: Atomic, Molecular and Optical Physics 29, 177 (1996).
  * Drake [2020] G. Drake, Accuracy in atomic and molecular data, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 223001 (2020).
  * Larson and Mavrogordatos [2022] J. Larson and T. Mavrogordatos,  _Jaynes-Cummings Model and Its Descendants: Modern research directions (IOP Series in Quantum Technology)_ (IOP Publishing Ltd, Bristol, 2022).
  * Drake [2022] G. W. Drake,  _Springer Handbook of Atomic, Molecular, and Optical Physics_ (Springer Handbooks, Würzburg, 2022).
  * Amini  _et al._ [2019] K. Amini, J. Biegert, F. Calegari, A. Chacón, M. F. Ciappina, A. Dauphin, D. K. Efimov, C. F. de Morisson Faria, K. Giergiel, P. Gniewek, _et al._ , Symphony on strong field approximation, Rep. Prog. Phys. 82, 116001 (2019).
  * Agostini  _et al._ [1979] P. Agostini, F. Fabre, G. Mainfray, G. Petite, and N. K. Rahman, Free-free transitions following six-photon ionization of xenon atoms, Phys. Rev. Lett. 42, 1127 (1979).
  * Kruit  _et al._ [1983] P. Kruit, J. Kimman, H. G. Muller, and M. J. van der Wiel, Electron spectra from multiphoton ionization of xenon at 1064, 532, and 355 nm, Phys. Rev. A 28, 248 (1983).
  * Becker  _et al._ [2002] W. Becker, F. Grasbon, R. Kopold, D. B. Milošević, G. G. Paulus, and H. Walther, Above-threshold ionization: From classical features to quantum effects, Adv. At. Mol. Opt. Phys. 48, 35 (2002).
  * McPherson  _et al._ [1987] A. McPherson, G. Gibson, H. Jara, U. Johann, T. S. Luk, I. A. McIntyre, K. Boyer, and C. K. Rhodes, Studies of multiphoton production of vacuum-ultraviolet radiation in the rare gases, J. Opt. Soc. Am. B 4, 595 (1987).
  * Ferray  _et al._ [1988] M. Ferray, A. L’Huillier, X. Li, L. Lompre, G. Mainfray, and M. C, Multiple-harmonic conversion of 1064 nm radiation in rare gases, Journal of Physics B: Atomic, Molecular and Optical Physics 21, L31 (1988).
  * Lewenstein and L’Huillier [2009] M. Lewenstein and A. L’Huillier,  _Strong Field Laser Physics_ , edited by T. Brabec (Springer, New York, 2009) Chap. Principles of Single Atom Physics: High-Order Harmonic Generation, Above-Threshold Ionization and Non-Sequential Ionization, pp. 147–183.
  * l’Huillier  _et al._ [1983] A. l’Huillier, L. A. Lompre, G. Mainfray, and C. Manus, Multiply charged ions induced by multiphoton absorption in rare gases at 0.53 $\mu$m, Phys. Rev. A 27, 2503 (1983).
  * Corkum [1993] P. B. Corkum, Plasma perspective on strong field multiphoton ionization, Phys. Rev. Lett. 71, 1994 (1993).
  * Walker  _et al._ [1994] B. Walker, B. Sheehy, L. F. DiMauro, P. Agostini, K. J. Schafer, and K. C. Kulander, Precision measurement of strong field double ionization of helium, Phys. Rev. Lett. 73, 1227 (1994).
  * Feuerstein  _et al._ [2001] B. Feuerstein, R. Moshammer, D. Fischer, A. Dorn, C. D. Schröter, J. Deipenwisch, J. R. Crespo Lopez-Urrutia, C. Höhr, P. Neumayer, J. Ullrich, H. Rottke, C. Trump, M. Wittmann, G. Korn, and W. Sandner, Separation of recollision mechanisms in nonsequential strong field double ionization of Ar: the role of excitation tunneling, Phys. Rev. Lett. 87, 043003 (2001).
  * Becker  _et al._ [2012a] W. Becker, X. Liu, P. Ho, and J. Eberly, Theories of photoelectron correlation in laser-driven multiple atomic ionization, Review of Modern Physics 84, 1011 (2012a).
  * Bandrauk and Kono [2006] A. Bandrauk and H. Kono,  _Advances in Multi-Photon Processes and Spectroscopy: Molecules in intense laser fields: Nonlinea Multiphoton Spectrocopy and Near-femtosecond to Sub-femtosecond (attosecond) Dynamics_ , Vol. 15 (World Scientific, Singapore, 2006) pp. 149–214.
  * Calegari  _et al._ [2016] F. Calegari, G. Sansone, S. Stagira, C. Vozzi, and M. Nisoli, Advances in attosecond science, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 062001 (2016).
  * Nisoli  _et al._ [2017] M. Nisoli, P. Decleva, F. Calegari, A. Palacios, and F. Martín, Attosecond electron dynamics in molecules, Chem. Rev. 117, 10760 (2017).
  * Liu [2015] Y. Liu, ed.,  _Advances of Atoms and Molecules in Strong Laser Fields_ (World Scientific, Singapore, 2015).
  * Lezius  _et al._ [1998] M. Lezius, S. Dobosz, D. Normand, and M. Schmidt, Explosion dynamics of rare gas clusters in strong laser fields, Phys. Rev. Lett. 80, 261 (1998).
  * Jungreuthmayer  _et al._ [2004] C. Jungreuthmayer, M. Geissler, J. Zanghellini, and T. Brabec, Microscopic analysis of large-cluster explosion in intense laser fields, Phys. Rev. Lett. 92, 133401 (2004).
  * Siedschlag  _et al._ [2004] C. Siedschlag, U. Saalmann, and J. Rost,  _Latest Advances in Atomic Cluster Collisions: CLUSTERS IN INTENSE LASER FIELDS_ , Vol. 15 (World Scientific, 2004) pp. 271–278.
  * Vampa and Brabec [2017] G. Vampa and T. Brabec, Merge of high harmonic generation from gases and solids and its implications for attosecond science, J. Phys. B 50, 083001 (2017).
  * Jiménez-Galán  _et al._ [2020] A. Jiménez-Galán, R. E. F. Silva, O. Smirnova, and M. Ivanov, Lightwave control of topological properties in 2d materials for sub-cycle and non-resonant valley manipulation, Nat. Phot. 14, 728 (2020).
  * Chacón  _et al._ [2020] A. Chacón, D. Kim, W. Zhu, S. Kelly, A. Dauphin, E. Pisanty, A. Maxwell, A. Picón, M. Ciappina, D. Kim, C. Ticknor, A. Saxena, and M. Lewenstein, Circular dichroism in higher-order harmonic generation: Heralding topological phases and transitions in chern insulators, Phys. Rev. B 102, 134115 (2020).
  * Eberly  _et al._ [1984] J. Eberly, K. Wódkiewicz, and B. Shore, Noise in strong laser-atom interactions: Phase telegraph noise, Phys. Rev. A 30, 2381 (1984).
  * Wódkiewicz and Eberly [1986] K. Wódkiewicz and J. Eberly, Effects of hidden transient noise processes on multiphoton absorption, Journal of the Optical Society of America B 3, 628 (1986).
  * Maxwell  _et al._ [2022] A. S. Maxwell, L. B. Madsen, and M. Lewenstein, Entanglement of Orbital Angular Momentum in Non-Sequential Double Ionization, Nat. Commun. [in press] (2022), arXiv:2111.10148 .
  * Artemyev  _et al._ [2017a] A. N. Artemyev, L. S. Cederbaum, and P. V. Demekhin, Impact of intense laser pulses on the autoionization dynamics of the $2s2p$ doubly excited state of He, Phys. Rev. A 96, 033410 (2017a).
  * Artemyev  _et al._ [2017b] A. N. Artemyev, L. S. Cederbaum, and P. V. Demekhin, Impact of two-electron dynamics and correlations on high-order-harmonic generation in He, Phys. Rev. A 95, 033402 (2017b).
  * Runge and Gross [1984] E. Runge and E. K. U. Gross, Density-Functional Theory for Time-Dependent Systems, Phys. Rev. Lett. 52, 997 (1984).
  * Castro  _et al._ [2015] A. Castro, A. Rubio, and E. K. U. Gross, Enhancing and controlling single-atom high-harmonic generation spectra: a time-dependent density-functional scheme, The European Physical Journal B 88, 191 (2015).
  * Woźniak  _et al._ [2021] A. P. Woźniak, M. Lesiuk, M. Przybytek, D. K. Efimov, J. S. Prauzner-Bechcicki, M. Mandrysz, M. Ciappina, E. Pisanty, J. Zakrzewski, M. Lewenstein, and R. Moszyński, A systematic construction of Gaussian basis sets for the description of laser field ionization and high-harmonic generation, The Journal of Chemical Physics 154, 094111 (2021).
  * Woźniak  _et al._ [2022] A. P. Woźniak, M. Przybytek, M. Lewenstein, and R. Moszyński, Effects of electronic correlation on the high harmonic generation in helium: a time-dependent configuration interaction singles vs time-dependent full configuration interaction study, arXiv:2201.09339 [physics] (2022).
  * Joachain  _et al._ [2012] C. J. Joachain, N. J. Kylstra, and R. Potvliege,  _Atoms in Intense Laser Fields_ (Cambridge University Press, Cambridge, Cambridge, 2012).
  * Leopold and Percival [1978] J. G. Leopold and I. C. Percival, Microwave Ionization and Excitation of Rydberg Atoms, Phys. Rev. Lett. 41, 944 (1978).
  * Percival [1977] I. C. Percival, Semiclassical theory of Bound States, in  _Advances in Chemical Physics_ (John Wiley & Sons, Ltd, 1977) pp. 1–61.
  * Grochmalicki  _et al._ [1991] J. Grochmalicki, M. Lewenstein, and K. Rz¸ażewski, Stabilization of atoms in superintense laser fields: Is it real?, Phys. Rev. Lett. 66, 1038 (1991).
  * Vampa  _et al._ [2014] G. Vampa, C. R. McDonald, G. Orlando, D. D. Klug, P. B. Corkum, and T. Brabec, Theoretical Analysis of High-Harmonic Generation in Solids, Phys. Rev. Lett. 113, 073901 (2014).
  * Osika  _et al._ [2017] E. N. Osika, A. Chacón, L. Ortmann, N. Suárez, J. A. Pérez-Hernández, B. Szafran, M. F. Ciappina, F. Sols, A. S. Landsman, and M. Lewenstein, Wannier-Bloch Approach to Localization in High-Harmonics Generation in Solids, Phys. Rev. X 7, 021017 (2017).
  * Bauer and Hansen [2018] D. Bauer and K. K. Hansen, High-Harmonic Generation in Solids with and without Topological Edge States, Phys. Rev. Lett. 120, 177401 (2018).
  * Jürß and Bauer [2019] C. Jürß and D. Bauer, High-harmonic generation in Su-Schrieffer-Heeger chains, Phys. Rev. B 99, 195428 (2019).
  * Baldelli  _et al._ [2022] N. Baldelli, U. Bhattacharya, D. González-Cuadra, M. Lewenstein, and T. Graß, Detecting Majorana Zero Modes via Strong Field Dynamics, arXiv:2202.03547 [cond-mat] (2022).
  * Alcalà  _et al._ [2022] J. Alcalà, U. Bhattacharya, J. Biegert, M. Ciappina, U. Elu, T. Graß, P. T. Grochowski, M. Lewenstein, A. Palau, T. P. H. Sidiropoulos, T. Steinle, and I. Tyulnev, High harmonic spectroscopy of quantum phase transitions in a high-${T}_{c}$ superconductor, arXiv:2201.09515 [cond-mat, physics:physics] (2022).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Electronic structure calculations with dynamical mean-field theory, Rev. Mod. Phys. 78, 865 (2006).
  * Vollhardt [2010] D. Vollhardt, Dynamical Mean‐Field Theory of Electronic Correlations in Models and Materials, AIP Conference Proceedings 1297, 339 (2010).
  * Lee [2007] P. A. Lee, From high temperature superconductivity to quantum spin liquid: progress in strong correlation physics, Rep. Prog. Phys. 71, 012501 (2007).
  * Eberly and Javanainen [1988] J. H. Eberly and J. Javanainen, Above-threshold ionisation, Eur J. Phys. 9, 265 (1988).
  * Eberly  _et al._ [1991] J. H. Eberly, J. Javanainen, and K. Rz¸ażewski, Above-threshold ionization, Physics Reports 204, 331 (1991).
  * Suárez  _et al._ [2015] N. Suárez, A. Chacón, M. F. Ciappina, J. Biegert, and M. Lewenstein, Above-threshold ionization and photoelectron spectra in atomic systems driven by strong laser fields, Phys. Rev. A 92, 063421 (2015).
  * Suárez  _et al._ [2016] N. Suárez, A. Chacón, M. F. Ciappina, B. Wolter, J. Biegert, and M. Lewenstein, Above-threshold ionization and laser-induced electron diffraction in diatomic molecules, Phys. Rev. A 94, 043423 (2016).
  * Suárez  _et al._ [2017] N. Suárez, A. Chacón, J. A. Pérez-Hernández, J. Biegert, M. Lewenstein, and M. F. Ciappina, High-order-harmonic generation in atomic and molecular systems, Phys. Rev. A 95, 033415 (2017).
  * Suárez  _et al._ [2018] N. Suárez, A. Chacón, E. Pisanty, L. Ortmann, A. S. Landsman, A. Picón, J. Biegert, M. Lewenstein, and M. F. Ciappina, Above-threshold ionization in multicenter molecules: The role of the initial state, Phys. Rev. A 97, 033415 (2018).
  * Suárez Rojas [2018] N. Suárez Rojas, Strong-field processes in atoms and polyatomic molecules, TDX (Tesis Doctorals en Xarxa) (2018).
  * Schuricke  _et al._ [2011] M. Schuricke, G. Zhu, J. Steinmann, K. Simeonidis, I. Ivanov, A. Kheifets, A. N. Grum-Grzhimailo, K. Bartschat, A. Dorn, and J. Ullrich, Strong-field ionization of lithium, Phys. Rev. A 83, 023413 (2011).
  * Taylor  _et al._ [2003] K. T. Taylor, J. S. Parker, K. J. Meharg, and D. Dundas, Laser-driven helium at 780 nm, The European Physical Journal D - Atomic, Molecular, Optical and Plasma Physics 26, 67 (2003).
  * Parker  _et al._ [2000] J. S. Parker, L. R. Moore, D. Dundas, and K. T. Taylor, Double ionization of helium at 390 nm, Journal of Physics B: Atomic, Molecular and Optical Physics 33, L691 (2000).
  * Moore  _et al._ [2001] L. Moore, J. Parker, D. Dundas, K. Cairns, and K. Taylor, The two electron response in laser driven helium: NATO Advanced Research Workshop on Super-Intense Laser-Atom Physics (2001) pp. 107–116.
  * Muller [2001] H. G. Muller, Non-sequential double ionization of helium and related wave-function dynamics obtained from a five-dimensional grid calculation., Optics Express 8, 417 (2001).
  * Grobe and Eberly [1992] R. Grobe and J. H. Eberly, Photoelectron spectra for a two-electron system in a strong laser field, Phys. Rev. Lett. 68, 10.1103/PhysRevLett.68.2905 (1992).
  * Liu  _et al._ [1999] W.-C. Liu, J. H. Eberly, S. L. Haan, and R. Grobe, Correlation Effects in Two-Electron Model Atoms in Intense Laser Fields, Phys. Rev. Lett. 83, 520 (1999).
  * Becker  _et al._ [2012b] W. Becker, X. Liu, P. J. Ho, and J. H. Eberly, Theories of photoelectron correlation in laser-driven multiple atomic ionization, Reviews of Modern Physics 84, 1011 (2012b).
  * de Morisson Faria and Liu [2011] C. F. de Morisson Faria and X. Liu, Electron–electron correlation in strong laser fields, Journal of Modern Optics 58, 1076 (2011).
  * Sacha and Eckhardt [2001] K. Sacha and B. Eckhardt, Pathways to double ionization of atoms in strong fields, Phys. Rev. A 63, 043414 (2001).
  * Prauzner-Bechcicki  _et al._ [2007] J. S. Prauzner-Bechcicki, K. Sacha, B. Eckhardt, and J. Zakrzewski, Time-Resolved Quantum Dynamics of Double Ionization in Strong Laser Fields, Phys. Rev. Lett. 98, 203002 (2007).
  * Efimov  _et al._ [2021] D. K. Efimov, A. Maksymov, M. Ciappina, J. S. Prauzner-Bechcicki, M. Lewenstein, and J. Zakrzewski, Three-electron correlations in strong laser field ionization, Optics Express 29, 26526 (2021).
  * Maxwell and Figueira de Morisson Faria [2016] A. S. Maxwell and C. Figueira de Morisson Faria, Controlling Below-Threshold Nonsequential Double Ionization via Quantum Interference, Phys. Rev. Lett. 116, 143001 (2016).
  * Guo and Åberg [1988] D. S. Guo and T. Åberg, Quantum electrodynamical approach to multiphoton ionisation in the high-intensity H field, Journal of Physics A: Mathematical and General 21, 4577 (1988).
  * Fu  _et al._ [2001] P. Fu, B. Wang, X. Li, and L. Gao, Interrelation between high-order harmonic generation and above-threshold ionization, Physical Review A - Atomic, Molecular, and Optical Physics 64, 063401 (2001).
  * Wang  _et al._ [2007] B. Wang, L. Gao, X. Li, D. S. Guo, and P. Fu, Frequency-domain theory of high-order above-threshold ionization based on nonperturbative quantum electrodynamics, Physical Review A - Atomic, Molecular, and Optical Physics 75, 063419 (2007).
  * Wang  _et al._ [2012] B. Wang, Y. Guo, J. Chen, Z.-C. Yan, and P. Fu, Frequency-domain theory of nonsequential double ionization in intense laser fields based on nonperturbative QED, Phys. Rev. A 85, 023402 (2012).
  * Gonoskov  _et al._ [2016] I. A. Gonoskov, N. Tsatrafyllis, I. K. Kominis, and P. Tzallas, Quantum optical signatures in strong-field laser physics: Infrared photon counting in high-order-harmonic generation, Scientific Reports 6, 32821 (2016).
  * Tsatrafyllis  _et al._ [2017] N. Tsatrafyllis, I. K. Kominis, I. A. Gonoskov, and P. Tzallas, High-order harmonics measured by the photon statistics of the infrared driving-field exiting the atomic medium, Nature Communications 8, 15170 (2017).
  * Gorlach  _et al._ [2020] A. Gorlach, O. Neufeld, N. Rivera, O. Cohen, and I. Kaminer, The quantum-optical nature of high harmonic generation, Nature Communications 11, 4598 (2020).
  * Rivera and Kaminer [2020] N. Rivera and I. Kaminer, Light–matter interactions with photonic quasiparticles, Nature Reviews Physics 2, 538 (2020).
  * Fuchs  _et al._ [2022] S. Fuchs, J. J. Abel, J. Nathanael, J. Reinhard, F. Wiesner, M. Wünsche, S. Skruszewicz, C. Rödel, D. Born, H. Schmidt, and G. G. Paulus, Photon counting of extreme ultraviolet high harmonics using a superconducting nanowire single-photon detector, Applied Physics B 128, 26 (2022).
  * Spasibko  _et al._ [2017] K. Y. Spasibko, D. A. Kopylov, V. L. Krutyanskiy, T. V. Murzina, G. Leuchs, and M. V. Chekhova, Multiphoton effects enhanced due to ultrafast photon-number fluctuations, Physical Review Letters 119, 223603 (2017).
  * Manceau  _et al._ [2019] M. Manceau, K. Y. Spasibko, G. Leuchs, R. Filip, and M. V. Chekhova, Indefinite-mean pareto photon distribution from amplified quantum noise, Physical review letters 123, 123606 (2019).
  * Gorlach  _et al._ [2022] A. Gorlach, M. E. Tzur, M. Birk, M. Krüger, N. Rivera, O. Cohen, and I. Kaminer, High harmonic generation driven by quantum light, arXiv preprint arXiv:2211.03188 (2022).
  * Varró [2021] S. Varró, Quantum Optical Aspects of High-Harmonic Generation, Photonics 8, 269 (2021).
  * Földi  _et al._ [2021] P. Földi, I. Magashegyi, A. Gombköto, and S. Varró, Describing High-Order Harmonic Generation Using Quantum Optical Models, Photonics 8, 263 (2021).
  * Gombköto  _et al._ [2021] A. Gombköto, P. Földi, and S. Varró, Quantum-optical description of photon statistics and cross correlations in high-order harmonic generation, Phys. Rev. A 104, 033703 (2021).
  * Lewenstein  _et al._ [2021] M. Lewenstein, M. Ciappina, E. Pisanty, J. Rivera-Dean, P. Stammer, Th. Lamprou, and P. Tzallas, Generation of optical schrödinger cat states in intense laser–matter interactions, Nature Physics 17, 1104 (2021).
  * Stammer  _et al._ [2022a] P. Stammer, J. Rivera-Dean, Th. Lamprou, E. Pisanty, M. F. Ciappina, P. Tzallas, and M. Lewenstein, High photon number entangled states and coherent state superposition from the extreme ultraviolet to the far infrared, Phys. Rev. Lett. 128, 123603 (2022a).
  * Rivera-Dean  _et al._ [2022] J. Rivera-Dean, Th. Lamprou, E. Pisanty, P. Stammer, A. Ordóñez, A. Maxwell, M. Ciappina, M. Lewenstein, and P. Tzallas, Strong laser fields and their power to generate controllable high-photon-number coherent-state superpositions, Phys. Rev. A 105, 033714 (2022).
  * Stammer [2022] P. Stammer, Theory of entanglement and measurement in high harmonic generation, arXiv preprint arXiv:2203.04354 (2022).
  * Lamprou  _et al._ [2021] T. Lamprou, R. Lopez-Martens, S. Haessler, I. Liontos, S. Kahaly, J. Rivera-Dean, P. Stammer, E. Pisanty, M. F. Ciappina, M. Lewenstein, and P. Tzallas, Quantum–optical spectrometry in relativistic laser–plasma interactions using the high-harmonic generation process: A proposal, Photonics 8, 192 (2021).
  * Rivera-Dean  _et al._ [2021] J. Rivera-Dean, P. Stammer, E. Pisanty, Th. Lamprou, P. Tzallas, M. Lewenstein, and M. F. Ciappina, New schemes for creating large optical Schrödinger cat states using strong laser fields, Journal of Computational Electronics 20, 2111 (2021).
  * Tannoudji  _et al._ [1992] C. C. Tannoudji, G. Grynberg, and J. Dupont-Roe,  _Atom-photon interactions_ (New York, NY (United States); John Wiley and Sons Inc., 1992).
  * Cohen-Tannoudji  _et al._ [1997] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Photons and Atoms-Introduction to Quantum Electrodynamics_ (Wiley-VCH, Weinheim, 1997).
  * Glauber and Lewenstein [1991] R. J. Glauber and M. Lewenstein, Quantum optics of dielectric media, Phys. Rev. A 43, 467 (1991).
  * Vogel and Welsch [2006] W. Vogel and D.-G. Welsch,  _Quantum optics_ (John Wiley & Sons, Weinheim, 2006).
  * Scully and Zubairy [2001] M. O. Scully and M. S. Zubairy,  _Quantum optics_ (Cambridge University Press, Cambridge, 2001).
  * Haroche and Raimond [2006] S. Haroche and J.-M. Raimond,  _Exploring the quantum: atoms, cavities and photons_ (Oxford Graduate Texts, Oxford, 2006).
  * Hudson [1974] R. L. Hudson, When is the wigner quasi-probability density non-negative?, Reports on Mathematical Physics 6, 249 (1974).
  * Kenfack and Życzkowski [2004] A. Kenfack and K. Życzkowski, Negativity of the Wigner function as an indicator of non-classicality, Journal of Optics B: Quantum and Semiclassical Optics 6, 396 (2004).
  * Asbóth  _et al._ [2004] J. K. Asbóth, P. Adam, M. Koniorczyk, and J. Janszky, Coherent-state qubits: entanglement and decoherence, The European Physical Journal D - Atomic, Molecular, Optical and Plasma Physics 30, 403 (2004).
  * Sanders [1992] B. C. Sanders, Entangled coherent states, Phys. Rev. A 45, 6811 (1992).
  * Gilchrist  _et al._ [2004a] A. Gilchrist, K. Nemoto, W. J. Munro, T. C. Ralph, S. Glancy, S. L. Braunstein, and G. J. Milburn, Schrödinger cats and their power for quantum information processing, Journal of Optics B: Quantum and Semiclassical Optics 6, S828 (2004a).
  * Jouguet  _et al._ [2013] P. Jouguet, S. Kunz-Jacques, A. Leverrier, P. Grangier, and E. Diamanti, Experimental demonstration of long-distance continuous-variable quantum key distribution, Nature Photonics 7, 378 (2013).
  * Lloyd and Braunstein [1999] S. Lloyd and S. L. Braunstein, Quantum Computation over Continuous Variables, Phys. Rev. Lett. 82, 1784 (1999).
  * Ralph  _et al._ [2003] T. C. Ralph, A. Gilchrist, G. J. Milburn, W. J. Munro, and S. Glancy, Quantum computation with optical coherent states, Phys. Rev. A 68, 042319 (2003).
  * Joo  _et al._ [2011] J. Joo, W. J. Munro, and T. P. Spiller, Quantum Metrology with Entangled Coherent States, Phys. Rev. Lett. 107, 083601 (2011).
  * Aichelburg and Grosse [1977] P. C. Aichelburg and H. Grosse, Exactly soluble system of relativistic two-body interaction, Physical Review D 16, 1900 (1977).
  * Aichelburg [1977] P. C. Aichelburg, Model for Relativistic Two-body Interaction with Radiation Reaction, Phys. Rev. Lett. 38, 451 (1977).
  * Rz̧ażewski and Żakowicz [1976] K. Rz̧ażewski and W. Żakowicz, Initial value problem and causality of radiating oscillator, Journal of Physics A: Mathematical and General 9, 1159 (1976).
  * Rz̧ażewski and Żakowicz [1980] K. Rz̧ażewski and W. Żakowicz, Initial value problem for two oscillators interacting with electromagnetic field, Journal of Mathematical Physics 21, 378 (1980).
  * Tong and Chu [1997] X.-M. Tong and S.-I. Chu, Density-functional theory with optimized effective potential and self-interaction correction for ground states and autoionizing resonances, Phys. Rev. A 55, 3406 (1997).
  * Tong and Lin [2005] X. M. Tong and C. D. Lin, Empirical formula for static field ionization rates of atoms and molecules by lasers in the barrier-suppression regime, Journal of Physics B: Atomic, Molecular and Optical Physics 38, 2593 (2005).
  * Lewenstein  _et al._ [1994] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, and P. B. Corkum, Theory of high-harmonic generation by low-frequency laser fields, Phys. Rev. A 49, 2117 (1994).
  * Sundaram and Milonni [1990] B. Sundaram and P. W. Milonni, High-order harmonic generation: simplified model and relevance of single-atom theories to experiment, Phys. Rev. A 41, 6571 (1990).
  * Bauer and Koval [2006] D. Bauer and P. Koval, Qprop: A Schrödinger-solver for intense laser–atom interaction, Computer Physics Communications 174, 396 (2006).
  * Nielsen and Chuang [2000] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, Cambridge, 2000).
  * Krause  _et al._ [1992] J. L. Krause, K. J. Schafer, and K. C. Kulander, High-order harmonic generation from atoms and ions in the high intensity regime, Phys. Rev. Lett. 68, 3535 (1992).
  * van Enk and Hirota [2001] S. J. van Enk and O. Hirota, Entangled coherent states: Teleportation and decoherence, Physical Review A 64, 022313 (2001).
  * Othman and Yevick [2018] A. Othman and D. Yevick, Quantum properties of the superposition of two nearly identical coherent states, International Journal of Theoretical Physics 57, 2293 (2018).
  * Madsen [2021] L. B. Madsen, Strong-field approximation for high-order harmonic generation in infrared laser pulses in the accelerated Kramers-Henneberger frame, Physical Review A 104, 033117 (2021).
  * Paulus  _et al._ [2003] G. G. Paulus, F. Lindner, H. Walther, A. Baltuška, E. Goulielmakis, M. Lezius, and F. Krausz, Measurement of the Phase of Few-Cycle Laser Pulses, Phys. Rev. Lett. 91, 253004 (2003).
  * Milošević  _et al._ [2006] D. B. Milošević, G. G. Paulus, D. Bauer, and W. Becker, Above-threshold ionization by few-cycle pulses, Journal of Physics B: Atomic, Molecular and Optical Physics 39, R203 (2006).
  * Lvovsky and Raymer [2009] A. I. Lvovsky and M. G. Raymer, Continuous-variable optical quantum-state tomography, Rev. Mod. Phys. 81, 299 (2009).
  * Breitenbach  _et al._ [1997] G. Breitenbach, S. Schiller, and J. Mlynek, Measurement of the quantum states of squeezed light, Nature 387, 471 (1997).
  * Skotiniotis  _et al._ [2017] M. Skotiniotis, W. Dür, and P. Sekatski, Macroscopic superpositions require tremendous measurement devices, Quantum 1, 34 (2017).
  * Chatziathanasiou  _et al._ [2006] S. Chatziathanasiou, S. Kahaly, E. Skantzakis, G. Sansone, R. Lopez-Martens, S. Haessler, K. Varju, G. D. Tsakiris, D. Charalambidis, and P. Tzallas, Generation of attosecond light pulses from gas and solid state media, Photonics 4, 26 (2006).
  * Paulus  _et al._ [2001] G. G. Paulus, F. Grasbon, H. Walther, P. Villoresi, M. Nisoli, S. Stagira, E. Priori, and S. D. Silvestri, Absolute-phase phenomena in photoionization with few–cycle laser pulses, Nature 414, 182 (2001).
  * Tsatrafyllis  _et al._ [2019] N. Tsatrafyllis, S. Kühn, M. Dumergue, P. Foldi, S. Kahaly, E. Cormier, I. A. Gonoskov, B. Kiss, K. Varju, S. Varro, and P. Tzallas, Quantum optical signatures in a strong laser pulse after interaction with semiconductors, Phys. Rev. Lett. 122, 193602 (2019).
  * Kim  _et al._ [2005] I. J. Kim, C. M. Kim, H. T. Kim, G. H. Lee, Y. S. Lee, J. Y. Park, D. J. Cho, and C. H. Nam, Highly efficient high–harmonic generation in an orthogonally polarized two–color laser field, Phys. Rev. Lett. 94, 243901 (2005).
  * Mauritsson  _et al._ [2006] J. Mauritsson, P. Johnsson, E. Gustafsson, A. L’Huillier, K. J. Schafer, and M. B. Gaarde, Attosecond pulse trains generated using two color laser fields, Phys. Rev. Lett. 97, 013001 (2006).
  * Fleischer  _et al._ [2014] A. Fleischer, O. Kfir, T. Diskin, P. Sidorenko, , and O. Cohen, Spin angular momentum and tunable polarization in highharmonic generation, Nat. Photonics 8, 543 (2014).
  * Gilchrist  _et al._ [2004b] A. Gilchrist, K. Nemoto, W. J. Munro, T. C. Ralph, S. Glancy, S. L. Braunstein, and G. J. Milburn, Schrödinger cats and their power for quantum information processing, J. Opt B: Quantum Semiclass. 6, 5828 (2004b).
  * Laghaout  _et al._ [2013] A. Laghaout, J. S. Neergaard-Nielsen, I. Rigas, C. Kragh, A. Tipsmark, and U. L. Andersen, Amplification of realistic Schr\”odinger-cat-state-like states by homodyne heralding, Phys. Rev. A 87, 043826 (2013).
  * Sychev  _et al._ [2017] D. V. Sychev, A. E. Ulanov, A. A. Pushkina, M. W. Richards, I. A. Fedorov, and A. I. Lvovsky, Enlargement of optical Schrödinger’s cat states, Nat. Phot. 11, 379 (2017).
  * Bachor and Ralph [2019] H. Bachor and T. Ralph,  _A Guide to Experiments in Quantum Optics_ (Wiley‐VCH Verlag, Weinheim, Germany, 2019).
  * Yuen and Chan [1983] H. P. Yuen and V. W. S. Chan, Noise in homodyne and heterodyne detection, Opt. Lett. 8, 177 (1983).
  * Lvovsky [2004] A. I. Lvovsky, Iterative maximum-likelihood reconstruction in quantum homodyne tomography, Journal of Optics B: Quantum and Semiclassical Optics 6, S556 (2004).
  * Schleich [2001] W. P. Schleich,  _Quantum Optics in Phase Space_ (Wiley-VHC Verlag, Weinheim, Germany, 2001).
  * Gerry and Knight [2005] C. Gerry and P. Knight,  _Introductory Quantum Optics_ (Cambridge University Press, Cambridge, UK, 2005).
  * Leonhardt [a] U. Leonhardt, Measuring the quantum state of light, (Eds., P. Knight and A. Miller, Cambridge University Press, UK, 1997) (a).
  * Leonhardt [b] U. Leonhardt, Essential quantum optics, (Cambridge University Press, UK, 2010) (b).
  * Herman [1980] G. T. Herman,  _Image Reconstruction from projections: The fundamentals of Computerized tomography_ (Academic Press, New York, 1980).
  * Banaszek [1999] K. Banaszek, Quantum homodyne tomography with a priori constraints, Phys. Rev. A 59, 4797 (1999).
  * M.Leonhardt [1993] M.Leonhardt, Quantum statistics of a lossless beam splitter: Su(2) symmetry in phase space, Phys. Rev. A 48, 3265 (1993).
  * Zhang  _et al._ [2021] M. Zhang, H. Kang, M. Wang, F. Xu, X. Su, and K. Peng, Quantifying quantum coherence of optical cat states, Photonics Research 9, 887 (2021).
  * Pisanty  _et al._ [2019a] E. Pisanty, G. J. Machado, V. Vicuña-Hernández, A. Picón, A. Celi, J. P. Torres, and M. Lewenstein, Knotting fractional-order knots with the polarization state of light, Nature Photonics 13, 569 (2019a).
  * Pisanty  _et al._ [2019b] E. Pisanty, L. Rego, J. San Román, A. Picón, K. M. Dorney, H. C. Kapteyn, M. M. Murnane, L. Plaja, M. Lewenstein, and C. Hernández-García, Conservation of torus-knot angular momentum in high-order harmonic generation, Phys. Rev. Lett. 122, 203201 (2019b).
  * Rego  _et al._ [2019] L. Rego, K. M. Dorney, N. J. Brooks, Q. L. Nguyen, C.-T. Liao, J. San Román, D. E. Couch, A. Liu, E. Pisanty, M. Lewenstein, _et al._ , Generation of extreme-ultraviolet beams with time-varying orbital angular momentum, Science 364, eaaw9486 (2019).
  * Neufeld  _et al._ [2019] O. Neufeld, D. Podolsky, and O. Cohen, Floquet group theory and its application to selection rules in harmonic generation, Nature Communications 10, 405 (2019).
  * Ayuso  _et al._ [2019] D. Ayuso, O. Neufeld, A. F. Ordonez, P. Decleva, G. Lerner, O. Cohen, M. Ivanov, and O. Smirnova, Synthetic chiral light for efficient control of chiral light–matter interaction, Nature Photonics 13, 866 (2019).
  * Stammer  _et al._ [2020] P. Stammer, S. Patchkovskii, and F. Morales, Evidence of ac-stark-shifted resonances in intense two-color circularly polarized laser fields, Physical Review A 101, 033405 (2020).
  * Ayuso  _et al._ [2021] D. Ayuso, A. F. Ordonez, P. Decleva, M. Ivanov, and O. Smirnova, Enantio-sensitive unidirectional light bending, Nature Communications 12, 3951 (2021).
  * Katsoulis  _et al._ [2021] G. P. Katsoulis, Z. Dube, P. Corkum, A. Staudte, and A. Emmanouilidou, Momentum scalar triple product as a measure of chirality in electron ionization dynamics of strongly-driven atoms, arXiv:2112.02670 [physics] (2021), arXiv: 2112.02670.
  * Mayer  _et al._ [2021] N. Mayer, M. Ivanov, and O. Smirnova, Imprinting chirality on atoms using synthetic chiral light, arXiv:2112.02658 [physics] (2021), arXiv: 2112.02658.
  * Dziarmaga  _et al._ [2012] J. Dziarmaga, W. H. Zurek, and M. Zwolak, Non-local quantum superpositions of topological defects, Nature Physics 8, 49 (2012).
  * Barth and Smirnova [2013] I. Barth and O. Smirnova, Spin-polarized electrons produced by strong-field ionization, Physical Review A 88, 013401 (2013).
  * Hartung  _et al._ [2016] A. Hartung, F. Morales, M. Kunitski, K. Henrichs, A. Laucke, M. Richter, T. Jahnke, A. Kalinin, M. Schöffler, L. P. H. Schmidt, M. Ivanov, O. Smirnova, and R. Dörner, Electron spin polarization in strong-field ionization of xenon atoms, Nature Photonics 10, 526 (2016).
  * Eckart  _et al._ [2018] S. Eckart, M. Kunitski, M. Richter, A. Hartung, J. Rist, F. Trinter, K. Fehre, N. Schlott, K. Henrichs, L. P. H. Schmidt, T. Jahnke, M. Schöffler, K. Liu, I. Barth, J. Kaushal, F. Morales, M. Ivanov, O. Smirnova, and R. Dörner, Ultrafast preparation and detection of ring currents in single atoms, Nature Physics 14, 701 (2018).
  * Kaushal and Smirnova [2018] J. Kaushal and O. Smirnova, Looking inside the tunnelling barrier III: Spin polarisation in strong field ionisation from orbitals with high angular momentum, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 174003 (2018).
  * Trabert  _et al._ [2018] D. Trabert, A. Hartung, S. Eckart, F. Trinter, A. Kalinin, M. Schöffler, L. P. H. Schmidt, T. Jahnke, M. Kunitski, and R. Dörner, Spin and Angular Momentum in Strong-Field Ionization, Physical Review Letters 120, 43202 (2018), arXiv:1711.03935 .
  * Ayuso  _et al._ [2022] D. Ayuso, A. F. Ordonez, and O. Smirnova, Ultrafast chirality: the road to efficient chiral measurements, arXiv:2203.00580 [physics] (2022), arXiv: 2203.00580.
  * Planas  _et al._ [2022] X. B. Planas, A. Ordóñez, M. Lewenstein, and A. S. Maxwell, Strong-field chiral imaging with twisted photoelectrons, arXiv preprint arXiv:2202.07289 (2022).
  * Cirac  _et al._ [1998] J. I. Cirac, M. Lewenstein, K. Mølmer, and P. Zoller, Quantum superposition states of bose-einstein condensates, Phys. Rev. A 57, 1208 (1998).
  * Giovannetti  _et al._ [2004] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum-enhanced measurements: beating the standard quantum limit, Science 306, 1330 (2004).
  * Tóth and Apellaniz [2014] G. Tóth and I. Apellaniz, Quantum metrology from a quantum information science perspective, Journal of Physics A: Mathematical and Theoretical 47, 424006 (2014).
  * Pezze  _et al._ [2018] L. Pezze, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, Reviews of Modern Physics 90, 035005 (2018).
  * Polino  _et al._ [2020] E. Polino, M. Valeri, N. Spagnolo, and F. Sciarrino, Photonic quantum metrology, AVS Quantum Science 2, 024703 (2020).
  * Göbel and Siegner [2015] E. O. Göbel and U. Siegner,  _Quantum Metrology: Foundation of Units and Measurements_ (John Wiley & Sons, 2015).
  * Barbieri [2022] M. Barbieri, Optical quantum metrology, PRX Quantum 3, 010202 (2022).
  * Bhattacharya  _et al._ [2022] U. Bhattacharya, S. Chaudhary, T. Grass, A. S. Johnson, S. Wall, and M. Lewenstein, Fermionic chern insulator from twisted light with linear polarization, Physical Review B 105, L081406 (2022).
  * Maxwell  _et al._ [2021] A. Maxwell, G. Armstrong, M. Ciappina, E. Pisanty, Y. Kang, A. Brown, M. Lewenstein, and C. F. de Morisson Faria, Manipulating twisted electrons in strong-field ionization, Faraday Discussions 228, 394 (2021).
  * Vrakking [2021] M. J. Vrakking, Control of attosecond entanglement and coherence, Phys. Rev. Lett. 126, 113203 (2021).
  * Koll  _et al._ [2022] L.-M. Koll, L. Maikowski, L. Drescher, T. Witting, and M. J. Vrakking, Experimental control of quantum-mechanical entanglement in an attosecond pump-probe experiment, Phys. Rev. Lett. 128, 043201 (2022).
  * Vrakking [2022] M. J. Vrakking, Ion-photoelectron entanglement in photoionization with chirped laser pulses, Journal of Physics B: Atomic, Molecular and Optical Physics (2022).
  * Bourassin-Bouchet  _et al._ [2020] C. Bourassin-Bouchet, L. Barreau, V. Gruson, J.-F. Hergott, F. Quéré, P. Salières, and T. Ruchon, Quantifying decoherence in attosecond metrology, Phys. Rev. X 10, 031048 (2020).
  * Laurell  _et al._ [2022] H. Laurell, D. Finkelstein-Shapiro, C. Dittel, C. Guo, R. Demjaha, M. Ammitzböll, R. Weissenbilder, L. Neoričić, S. Luo, M. Gisselbrecht, _et al._ , Continuous variable quantum state tomography of photoelectrons, arXiv preprint arXiv:2202.06798 (2022).
  * Busto  _et al._ [2021] D. Busto, H. Laurell, D. F. Shapiro, C. Alexandridi, M. Isinger, S. Nandi, R. Squibb, M. Turconi, S. Zhong, C. Arnold, _et al._ , Probing electronic decoherence with high-resolution attosecond photoelectron interferometry, arXiv preprint arXiv:2111.12037 (2021).
  * Stammer  _et al._ [2022b] P. Stammer, J. Rivera-Dean, A. Maxwell, T. Lamprou, A. Ordóñez, M. F. Ciappina, P. Tzallas, and M. Lewenstein, Quantum electrodynamics of ultra-intense laser-matter interactions, Zenodo (2022b).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nat. Meth. 17, 261 (2020).

##  Appendix A Field commutation relations

The commutation relation between the canonical field variables
$\mathbf{A}(\mathbf{r})$ and $\mathbf{\Pi}(\mathbf{r})$ defined respectively
in Eqs. (3) and (5) is given by

$\displaystyle\left[A_{m}(\mathbf{r}),\Pi_{n}(\mathbf{r}^{\prime})\right]=i\hbar\delta^{\perp}_{mn}(\mathbf{r}-\mathbf{r}^{\prime}),$
(97)

with the transverse $\delta$-function

$\displaystyle\delta_{mn}^{\perp}(\mathbf{r}-\mathbf{r}^{\prime})=\frac{1}{(2\pi)^{3}}\int\differential^{3}k\left[\delta_{mn}-\frac{k_{m}k_{n}}{k^{2}}\right]e^{i\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}^{\prime})}.$
(98)

These canonical conjugate field variables are related to the field operators,
$\mathbf{B}=\nabla\times\mathbf{A}$, and
$\mathbf{E}^{\perp}=-\mathbf{\Pi}/\epsilon_{0}$, of the magnetic field and the
transverse part of electric field, respectively. The free-field Hamiltonian
$H_{f}$ can then be expressed in terms of the physical operators

$\displaystyle H_{f}=\frac{1}{2}\int
d^{3}r\left[\epsilon_{0}\mathbf{E}^{2}(\mathbf{r})+\frac{1}{\mu_{0}}\mathbf{B}^{2}(\mathbf{r})\right],$
(99)

which has the same form as the energy of the classical electromagnetic field,
and it is required that the observables obey the commutation relation

$\displaystyle\left[E_{m}(\mathbf{r}),B_{n}(\mathbf{r}^{\prime})\right]=-\epsilon_{mnl}\frac{i\hbar}{\epsilon_{0}}\partialderivative{x_{l}}\delta(\mathbf{r}-\mathbf{r}^{\prime}).$
(100)

##  Appendix B Interaction Hamiltonian

The minimal-coupling Hamiltonian given in Eq. (17) can be written in terms of
the transformed dynamical variables from (23), (24), to obtain the multipolar-
coupling Hamiltonian

$\displaystyle H^{\prime}$
$\displaystyle=\frac{1}{2}\int\differential^{3}r\left\\{\frac{1}{\epsilon_{0}}\left[\mathbf{\Pi}^{\prime}(\mathbf{r})+\mathbf{P}_{\perp}(\mathbf{r})\right]^{2}+\frac{1}{\mu_{0}}\left[\nabla\times\mathbf{A}(\mathbf{r})\right]^{2}\right\\}$
(101)
$\displaystyle\quad+\sum_{\alpha}\bigg{[}\frac{1}{2m}\Big{\\{}\mathbf{p}^{\prime}_{\alpha}-e\int_{0}^{1}\differential
s\
s(\mathbf{r}_{\alpha}-\mathbf{R}_{\alpha})\times\left[\nabla\times\mathbf{A}(\mathbf{R}_{\alpha}+s(\mathbf{r}_{\alpha}-\mathbf{R}_{\alpha}))\right]\Big{\\}}^{2}+V_{\rm
at}(\mathbf{r}_{\alpha},\mathbf{R}_{\alpha}))\bigg{]},$

where, as mentioned in Section III, $V_{\rm
at}(\mathbf{r}_{\alpha},\mathbf{R}_{\alpha})$ is the effective potential felt
by the single active electron, and $\mathbf{\Pi}^{\prime}$ and
$\mathbf{p}^{\prime}$ are given by Eqs. (23) and (24) respectively.

Using the mode expansion of the dynamical field variables of the transverse
part of the electromagnetic field, the total Hamiltonian can be decomposed in
the following way

$\displaystyle H^{\prime}=H_{f}^{\prime}+H_{\rm
at}^{\prime}+H_{int}^{\prime},$ (102)

where the free-field Hamiltonian $H_{f}^{\prime}$ has the same form as (2) but
in terms of the new mode operators

$\displaystyle H_{f}^{\prime}$
$\displaystyle=\frac{1}{2}\int\differential^{3}r\left\\{\frac{1}{\epsilon_{0}}{\mathbf{\Pi}^{\prime}}^{2}(\mathbf{r})+\frac{1}{\mu_{0}}\left[\nabla\times\mathbf{A}(\mathbf{r})\right]^{2}\right\\}=\sum_{\mu}\int\differential^{3}k\,\hbar\omega_{k}\left({a^{\prime}}^{\dagger}_{\bf{k}\mu}a^{\prime}_{\bf{k}\mu}+\frac{1}{2}\right).$
(103)

We have separated the purely electronic (atomic or molecular) Hamiltonian

$\displaystyle H_{\rm at}^{\prime}$
$\displaystyle=\sum_{\alpha}\left[\frac{{\mathbf{p}^{\prime}}^{2}_{\alpha}}{2m}+V_{\rm
at}(\mathbf{r}_{\alpha},\mathbf{R}_{\alpha})+\frac{1}{2\epsilon_{0}}\int\differential^{3}r\mathbf{P}_{\perp}^{2}(\mathbf{r})\right]=\sum_{\alpha}\left[\frac{{\mathbf{p}^{\prime}}^{2}_{\alpha}}{2m_{\alpha}}+\tilde{V}_{\rm
at}(\mathbf{r}_{\alpha},\mathbf{R}_{\alpha})\right],$ (104)

where $\tilde{V}_{\rm at}(\mathbf{r}_{\alpha},\mathbf{R}_{\alpha})$ is the
final effective potential for single active electrons that includes effects of
polarization. Finally, the interaction Hamiltonian reads

$\displaystyle H_{int}^{\prime}=$
$\displaystyle-\frac{e}{\epsilon_{0}}\sum_{\alpha}\int_{0}^{1}\differential
s(\mathbf{r}_{\alpha}-\mathbf{R}_{\alpha})\cdot\mathbf{\Pi}^{\prime}(\mathbf{R}_{\alpha}+s(\mathbf{r}_{\alpha}-\mathbf{R}_{\alpha}))$
(105) $\displaystyle+\frac{e}{2m}\sum_{\alpha}\int_{0}^{1}\differential s\
s\left\\{[(\mathbf{r}_{\alpha}-\mathbf{R}_{\alpha})\times\mathbf{p}^{\prime}_{\alpha}][\nabla\times\mathbf{A}(\mathbf{R}_{\alpha}+s(\mathbf{r}_{\alpha}-\mathbf{R}_{\alpha}))]+\operatorname{h.c.}\right\\}$
$\displaystyle+\frac{e^{2}}{2m}\sum_{\alpha}\left\\{\int_{0}^{1}\differential
s\
s(\mathbf{r}_{\alpha}-\mathbf{R}_{\alpha})\times[\nabla\times\mathbf{A}(\mathbf{R}_{\alpha}+s(\mathbf{r}_{\alpha}-\mathbf{R}_{\alpha}))]\right\\}^{2}.$

The first term is the most relevant in our analysis since it couples the
transverse displacement field
$\mathbf{\Pi}^{\prime}(\mathbf{r})=-[\epsilon_{0}\mathbf{E}_{\perp}(\mathbf{r})+\mathbf{P}_{\perp}(\mathbf{r})]$
with the atomic polarization $\mathbf{P}(\mathbf{r})$. The second and third
term describe the interaction of the paramagnetic magnetization with the
magnetic induction field
$\mathbf{B}(\mathbf{r})=\nabla\times\mathbf{A}(\mathbf{r})$, and the
diamagnetic energy of the system which is quadratic in the magnetic induction
field, respectively.

##  Appendix C HHG spectrum

The coherent contribution to the HHG spectrum is proportional to

$\displaystyle S(\omega_{\mathbf{k},\mu})$
$\displaystyle\propto\lim_{t\to\infty}\absolutevalue{\chi_{\mathbf{k},\mu}(t)}^{2}$
(106)
$\displaystyle=\tilde{g}(k)^{2}N\absolutevalue{\mathbf{\epsilon}_{\mathbf{k},\mu}\cdot\expectationvalue{\mathbf{d}(\omega_{\mathbf{k}})}}^{2}+\tilde{g}(k)^{2}\sum_{\alpha\neq\alpha^{\prime}}e^{-i\mathbf{k}(\mathbf{R}_{\alpha}-\mathbf{R}_{\alpha^{\prime}})}\mathbf{\epsilon}_{\mathbf{k},\mu}\cdot\expectationvalue{\mathbf{d}_{\alpha}(\omega_{\mathbf{k}})}\mathbf{\epsilon}_{\mathbf{k},\mu}\cdot\expectationvalue{\mathbf{d}_{\alpha^{\prime}}(\omega_{\mathbf{k}})}^{*},$

where we have assumed that each of the $N$ atoms contributes equally to the
HHG spectrum such that

$\displaystyle\sum_{\alpha=\alpha^{\prime}}\absolutevalue{\mathbf{\epsilon}_{\mathbf{k},\mu}\cdot\expectationvalue{\mathbf{d}_{\alpha}(\omega_{\mathbf{k}})}}^{2}=N\absolutevalue{\mathbf{\epsilon}_{\mathbf{k},\mu}\cdot\expectationvalue{\mathbf{d}(\omega_{\mathbf{k}})}}^{2}.$
(107)

Assuming that the $N$ atoms contribute equally to the HHG spectrum of course
implies the additional assumption that all atoms see the same radiation field
of the driving source. In realistic experiments this is almost never the case,
and focal averaging over the spatial beam profile needs to be taken into
account. However, for the concepts and methods introduced in this manuscript
this only plays a minor role since the post-selection scheme does not depend
on the particular structure of the HHG spectra, and holds true for generic
harmonic amplitudes $\chi_{\mathbf{k},\mu}$.

Following the same assumptions that lead to the dipole approximation, i.e.
$\mathbf{k}\cdot\mathbf{R}_{\alpha}\ll 1$, and again assuming that all $N$
atoms equally contribute to the HHG process, the second term in (106) is given
by

$\displaystyle\tilde{g}(k)^{2}$
$\displaystyle\sum_{\alpha\neq\alpha^{\prime}}e^{-i\mathbf{k}(\mathbf{R}_{\alpha}-\mathbf{R}_{\alpha^{\prime}})}\mathbf{\epsilon}_{\mathbf{k},\mu}\cdot\expectationvalue{\mathbf{d}_{\alpha}(\omega_{\mathbf{k}})}\mathbf{\epsilon}_{\mathbf{k},\mu}\cdot\expectationvalue{\mathbf{d}_{\alpha^{\prime}}(\omega_{\mathbf{k}})}^{*}=\tilde{g}(k)^{2}(N^{2}-N)\absolutevalue{\mathbf{\epsilon}_{\mathbf{k},\mu}\cdot\expectationvalue{\mathbf{d}(\omega_{\mathbf{k}})}}^{2},$
(108)

such that the spectrum is given by

$\displaystyle
S(\omega_{\mathbf{k},\mu})\propto\lim_{t\to\infty}\absolutevalue{\chi_{\mathbf{k},\mu}(t)}^{2}=\tilde{g}(k)^{2}N^{2}\absolutevalue{\mathbf{\epsilon}_{\mathbf{k},\mu}\cdot\expectationvalue{\mathbf{d}(\omega_{\mathbf{k}})}}^{2}.$
(109)

##  Appendix D Derivation of electronic transition matrix elements

In the main part of the manuscript we need the time-dependent dipole
transition matrix elements between the ground state $\ket{g}$ and a continuum
state $\ket{\mathbf{v}}$, or between two different continuum states, see for
instance (69) or (57), respectively. To obtain the transition matrix elements
we use that $\mathbf{r}(t)=U_{sc}^{\dagger}(t)\mathbf{r}U_{sc}(t)$, and we
will further use that the action of the propagator on the electronic state is
known from the semi-classical theory [118]

$\displaystyle\ket{\psi_{g/\mathbf{v}}(t)}$
$\displaystyle=U_{sc}(t)\ket{g/\mathbf{v}}=e^{iI_{p}t/\hbar}\left(a_{g/\mathbf{v}}(t)\ket{g}+\int\differential^{3}v^{\prime}b_{g/\mathbf{v}}(\mathbf{v}^{\prime},t)\ket{\mathbf{v}^{\prime}}\right),$
(110)

which are the solutions of the time-dependent Schrödinger equation with the
ground state $\ket{g}$ or the continuum state $\ket{\mathbf{v}}$ as the
initial condition, respectively.

###  D.1 Electron ground and continuum state population amplitudes

The amplitudes (110) of the ground and continuum states obey the following
coupled differential equations

$\derivative{t}a_{g/\mathbf{v}}(t)=\frac{i}{\hbar}\int\differential^{3}v^{\prime}\mathbf{E}_{cl}(t)\cdot\mathbf{d}^{*}(\mathbf{v}^{\prime})b_{g/\mathbf{v}}(\mathbf{v}^{\prime},t)$
(111) $\displaystyle\derivative{t}b_{g/\mathbf{v}}(\mathbf{v}^{\prime},t)=$
$\displaystyle-\frac{i}{\hbar}\left(\frac{\mathbf{v^{\prime}}^{2}}{2m}+I_{p}\right)b_{g/\mathbf{v}}(\mathbf{v}^{\prime},t)+\frac{i}{\hbar}\mathbf{E}_{cl}(t)\cdot\mathbf{d}(\mathbf{v}^{\prime})a_{g/\mathbf{v}}(t)$
(112)
$\displaystyle+e\mathbf{E}_{cl}(t)\cdot\nabla_{\mathbf{v}^{\prime}}b_{g/\mathbf{v}}(\mathbf{v}^{\prime},t)+i\mathbf{E}_{cl}(t)\cdot\int\differential^{3}v^{\prime\prime}b_{g/\mathbf{v}}(\mathbf{v}^{\prime\prime},t)\mathbf{g}(\mathbf{v}^{\prime},\mathbf{v}^{\prime\prime}),$

where
$\mathbf{d}(\mathbf{v}^{\prime})=-e\bra{\mathbf{v}^{\prime}}\mathbf{r}\ket{g}$.
However, to solve this coupled differential equations we shall perform the
following approximation by neglecting re-scattering events and only consider
the contribution from direct electrons, i.e. we neglect the continuum-
continuum transitions in (112) (term with
$\mathbf{g}(\mathbf{v}^{\prime},\mathbf{v}^{\prime\prime})$). We can therefore
solve the differential equation for the continuum amplitude with the initial
condition, $b_{g}(\mathbf{v}^{\prime},t_{0})=0$ and $a_{g}(t_{0})=1$, i.e.
initial ground state population

$\displaystyle\derivative{t}b_{g}(\mathbf{v}^{\prime},t)$
$\displaystyle=-\frac{i}{\hbar}\left(\frac{\mathbf{v^{\prime}}^{2}}{2m}+I_{p}\right)b_{g}(\mathbf{v}^{\prime},t)+\frac{i}{\hbar}\mathbf{E}_{cl}(t)\cdot\mathbf{d}(\mathbf{v}^{\prime})a_{g}(t)+e\mathbf{E}_{cl}(t)\cdot\nabla_{\mathbf{v}^{\prime}}b_{g}(\mathbf{v}^{\prime},t),$
(113)

with the solution given by

$\displaystyle b_{g}(\mathbf{p}^{\prime},t)$
$\displaystyle=\frac{i}{\hbar}\int_{t_{0}}^{t}\differential
t^{\prime}\mathbf{E}_{cl}(t^{\prime})\cdot\mathbf{d}\left(\mathbf{p}^{\prime}+\frac{e}{c}\mathbf{A}_{cl}(t^{\prime})\right)\exp{-\frac{i}{\hbar}S(\mathbf{p}^{\prime},t,t^{\prime})}a_{g}(t^{\prime}),$
(114)

where we have substituted the canonical momentum
$\mathbf{p}^{\prime}=m\mathbf{v}^{\prime}-\frac{e}{c}\mathbf{A}(t)$, such that
the gauge invariant kinetic momentum
$m\dot{\mathbf{r}}(t)=\mathbf{p}+\frac{e}{c}\mathbf{A}(t)$ appears in the
expression. With this, we denote by $S(\mathbf{p},t,t^{\prime})$ the semi-
classical action given by

$S(\mathbf{p},t,t^{\prime})=\int^{t}_{t^{\prime}}\differential
t^{\prime\prime}\bigg{[}\dfrac{1}{2m}\Big{(}\mathbf{p}^{\prime}+\dfrac{e}{c}\mathbf{A}_{cl}(t^{\prime\prime})\Big{)}^{2}+I_{p}\bigg{]}.$
(115)

This can now be inserted into the differential equation for the ground state
amplitude (111) and we obtain

$\displaystyle\derivative{t}a_{g}(t)=-\int_{t_{0}}^{t}\differential
t^{\prime}\gamma(t,t^{\prime})a_{g}(t^{\prime}),$ (116)

where

$\displaystyle\gamma(t,t^{\prime})$
$\displaystyle=\frac{1}{\hbar^{2}}\int\differential^{3}p\mathbf{E}_{cl}(t)\cdot\mathbf{d}^{*}(\mathbf{p}+\frac{e}{c}\mathbf{A}_{cl}(t))\mathbf{E}_{cl}(t^{\prime})\cdot\mathbf{d}\left(\mathbf{p}+\frac{e}{c}\mathbf{A}_{cl}(t^{\prime})\right)\exp{-\frac{i}{\hbar}S(\mathbf{p},t,t^{\prime})}.$
(117)

We now assume that the rate of change of $a_{g}(t)$ is slow compared to the
field, such that we can replace $a_{g}(t^{\prime})$ by $a_{g}(t)$ in the
integrand, and the solution is given by

$\displaystyle a_{g}(t)$ $\displaystyle=\exp{-\int_{t_{0}}^{t}\differential
t^{\prime}W(t^{\prime})}a_{g}(t_{0})=\exp{-\int_{t_{0}}^{t}\differential
t^{\prime}W(t^{\prime})},$ (118)

with the ionization rate given by

$\displaystyle W(t)=\int_{t_{0}}^{t}\differential
t^{\prime}\gamma(t,t^{\prime}).$ (119)

We further need the amplitudes with the continuum state as the initial
condition
$b_{\mathbf{v}}(\mathbf{v}^{\prime},t_{0})=\delta(\mathbf{v}-\mathbf{v^{\prime}})$
and $a_{\mathbf{v}}(t_{0})=0$. The corresponding solution of (113) for
$b_{\mathbf{v}}(\mathbf{v}^{\prime},t)$ is given by

$\displaystyle b_{\mathbf{v}}(\mathbf{v}^{\prime},t)$
$\displaystyle=\frac{i}{\hbar}\int_{t_{0}}^{t}\differential
t^{\prime}\mathbf{E}_{cl}(t^{\prime})\cdot\mathbf{d}\left(\mathbf{p}^{\prime}+\frac{e}{c}\mathbf{A}_{cl}(t^{\prime})\right)\exp{-\frac{i}{\hbar}S(\mathbf{p}^{\prime},t,t^{\prime})}a_{\mathbf{v}}(t^{\prime})+\exp{-\frac{i}{\hbar}S(\mathbf{p}^{\prime},t,t_{0})}\delta(\mathbf{v}-\mathbf{v^{\prime}}),$
(120)

where the last term takes into account the initial condition
$b_{\mathbf{v}}(\mathbf{v}^{\prime},t_{0})=\delta(\mathbf{v}-\mathbf{v^{\prime}})$.
Under the same assumptions for the solution of the differential equation of
$a_{\mathbf{v}}(t)$ in (111), we find that for $a_{\mathbf{v}}(t_{0})=0$

$\displaystyle a_{\mathbf{v}}(t)$
$\displaystyle=\frac{i}{\hbar}\int\differential^{3}v^{\prime}\int_{t_{0}}^{t}\differential
t^{\prime}\mathbf{E}_{cl}(t^{\prime})\cdot\mathbf{d}^{*}(\mathbf{v}^{\prime})\exp{-\frac{i}{\hbar}S(\mathbf{p}^{\prime},t^{\prime},t_{0})}\exp{-\int_{t^{\prime}}^{t}\differential\tau
W(\tau)}\delta(\mathbf{v}-\mathbf{v^{\prime}})$ (121)
$\displaystyle=\frac{i}{\hbar}\int_{t_{0}}^{t}\differential
t^{\prime}\mathbf{E}_{cl}(t^{\prime})\cdot\mathbf{d}^{*}(\mathbf{v})\exp{-\frac{i}{\hbar}S(\mathbf{p},t^{\prime},t_{0})}\exp{-\int_{t^{\prime}}^{t}\differential\tau
W(\tau)}.$

We shall assume that the probability of recombination to the ground state is
small when the electron is initially in the continuum. This means we
approximate $a_{\mathbf{v}}(t)\approx 0$, which simplifies the amplitude of
the continuum population

$\displaystyle b_{\mathbf{p}}(\mathbf{p}^{\prime},t)$
$\displaystyle=\exp{-\frac{i}{\hbar}S(\mathbf{p}^{\prime},t,t_{0})}\delta(\mathbf{p}-\mathbf{p}^{\prime}).$
(122)

We further assume that the depletion of the ground state can be neglected such
that $a_{g}(t)\simeq 1$ and accordingly the continuum population (114) reads

$\displaystyle b_{g}(\mathbf{p}^{\prime},t)$
$\displaystyle=\frac{i}{\hbar}\int_{t_{0}}^{t}\differential
t^{\prime}\mathbf{E}_{cl}(t^{\prime})\cdot\mathbf{d}\left(\mathbf{p}^{\prime}(t)+\frac{e}{c}\mathbf{A}_{cl}(t^{\prime})\right)\exp{-\frac{i}{\hbar}S(\mathbf{p},t,t^{\prime})}.$
(123)

We have expressed both amplitudes in terms of the gauge invariant kinetic
momentum $\mathbf{p}=\mathbf{v}-\frac{e}{c}\mathbf{A}(t)$.

###  D.2 Bound-continuum transition matrix element

In (69) we need the time-dependent transition matrix element between the
ground state $\ket{g}$ and the continuum state $\mathbf{v_{\alpha}}$ for the
process of ATI. We can use (110), and have

$\displaystyle\bra{\mathbf{v}}\mathbf{r}(t)\ket{g}=\bra{\psi_{\mathbf{v}}(t)}\mathbf{r}\ket{\psi_{g}(t)}\simeq\int\differential^{3}v^{\prime}b_{\mathbf{v}}^{*}(\mathbf{v}^{\prime},t)\bra{\mathbf{v}^{\prime}}\mathbf{r}\ket{g},$
(124)

where we have used that $a_{\mathbf{v}}(t)\approx 0$ and $a_{g}(t)\simeq 1$.
And we have further neglected C-C transition matrix elements. Inserting (122)
we have

$\displaystyle\bra{\mathbf{v}}\mathbf{r}_{\alpha}(t)\ket{g}\simeq\exp{\frac{i}{\hbar}S(\mathbf{p},t,t_{0})}\bra{\mathbf{p}+e/c\mathbf{A}(t)}\mathbf{r}\ket{g}.$
(125)

###  D.3 Continuum-Continuum transition matrix element

For the time-dependent transition matrix element between the two continuum
states in (57) we have

$\displaystyle\bra{\mathbf{v}}\mathbf{r}(t)\ket{\mathbf{v}^{\prime}}$
$\displaystyle=\bra{\psi_{\mathbf{v}}(t)}\mathbf{r}\ket{\psi_{\mathbf{v}^{\prime}}(t)}\simeq\int\differential^{3}v^{\prime\prime}\int\differential^{3}v^{\prime\prime\prime}b^{*}_{\mathbf{v}}(\mathbf{v}^{\prime\prime},t)b_{\mathbf{v}^{\prime}}(\mathbf{v}^{\prime\prime\prime},t)\matrixelement{\mathbf{v}^{\prime\prime}}{\mathbf{r}}{\mathbf{v}^{\prime\prime\prime}},$
(126)

where we have again used that $a_{\mathbf{v}}(t)\approx 0$ and $a_{g}(t)\simeq
1$.

##  Appendix E Field observables for direct ATI

The pure state of the EM field conditioned on ATI (68) is given by

$\displaystyle\ket{\phi_{ATI}(t,\mathbf{p})}$
$\displaystyle=\frac{-ie}{\hbar}\int_{t_{0}}^{t}\differential
t^{\prime}\prod_{\mathbf{k},\mu}e^{i\varphi_{\mathbf{k}\mu}(t,t^{\prime},\mathbf{v}_{\alpha})}D[\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})]\mathbf{E}_{Q}(t^{\prime},\mathbf{R}_{\alpha})\cdot\bra{\mathbf{v}_{\alpha}}\mathbf{r}_{\alpha}(t^{\prime})\ket{g}$
(127) $\displaystyle\hskip
170.71652pt\bigotimes_{\mathbf{k},\mu\simeq\mathbf{k}_{L,\mu}}\ket{\chi_{\mathbf{k},\mu}(t^{\prime})}\bigotimes_{\mathbf{k},\mu\gg\mathbf{k}_{L,\mu}}\ket{\chi_{\mathbf{k},\mu}(t^{\prime})}.$

We can now use that

$\displaystyle\prod_{\mathbf{k},\mu}e^{i\varphi_{\mathbf{k}\mu}(t,t^{\prime},\mathbf{v}_{\alpha})}D[\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})]\mathbf{E}_{Q}(t^{\prime},\mathbf{R}_{\alpha})=\left[\mathbf{E}_{Q}(t^{\prime},\mathbf{R}_{\alpha})-\mathbf{E}_{\delta}(t^{\prime},\mathbf{R}_{\alpha})\right]\prod_{\mathbf{k},\mu}e^{i\varphi_{\mathbf{k}\mu}(t,t^{\prime},\mathbf{v}_{\alpha})}D[\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})],$
(128)

where

$\displaystyle\mathbf{E}_{\delta}(t^{\prime},\mathbf{R}_{\alpha})$
$\displaystyle=i\sum_{\mu}\int\differential^{3}k\tilde{g}(k)\epsilon_{k\mu}\Big{[}\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})e^{-i\omega_{k}t+i\mathbf{k}\cdot\mathbf{R}_{\alpha}}-\text{c.c.}\Big{]}$
(129)

such that the state of the field now reads

$\displaystyle\ket{\phi_{ATI}(t,\mathbf{p})}=$
$\displaystyle\frac{-ie}{\hbar}\int_{t_{0}}^{t}\differential
t^{\prime}\left[\mathbf{E}_{Q}(t^{\prime},\mathbf{R}_{\alpha})-\mathbf{E}_{\delta}(t^{\prime},\mathbf{R}_{\alpha})\right]\cdot\bra{\mathbf{v}_{\alpha}}\mathbf{r}_{\alpha}(t^{\prime})\ket{g}$
(130)
$\displaystyle\bigg{[}\bigotimes_{\mathbf{k},\mu\simeq\mathbf{k}_{L,\mu}}e^{i[\varphi_{\mathbf{k}\mu}(t,t^{\prime},\mathbf{v}_{\alpha})+\varphi_{\mathbf{k}\mu}^{\chi,\delta}(t,t^{\prime})]}\ket{\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})}\bigg{]}$
$\displaystyle\bigg{[}\bigotimes_{\mathbf{k},\mu\gg\mathbf{k}_{L,\mu}}e^{i[\varphi_{\mathbf{k}\mu}(t,t^{\prime},\mathbf{v}_{\alpha})+\varphi_{\mathbf{k}\mu}^{\chi,\delta}(t,t^{\prime})]}\ket{\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})}\bigg{]},$

where

$\displaystyle\varphi_{\mathbf{k}\mu}^{\chi,\delta}(t,t^{\prime})=\operatorname{Im}\left[\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\chi_{\mathbf{k},\mu}^{*}(t^{\prime})\right].$
(131)

We can now do the transformation back to the original laboratory frame

$\displaystyle\ket{\Phi_{ATI}(t,\mathbf{p})}=e^{-iH_{f}t}D(\alpha_{\mathbf{k}_{L}})\ket{\phi_{ATI}(t,\mathbf{p})},$
(132)

and we find that

$\displaystyle\ket{\Phi_{ATI}(t,\mathbf{p})}$
$\displaystyle=\frac{-ie}{\hbar}\int_{t_{0}}^{t}\differential
t^{\prime}\left[\mathbf{E}_{Q}(t^{\prime}-t,\mathbf{R}_{\alpha})-\mathbf{E}_{\delta}(t^{\prime},\mathbf{R}_{\alpha})-\mathbf{E}_{cl}(t^{\prime},\mathbf{R}_{\alpha})\right]\cdot\bra{\mathbf{v}_{\alpha}}\mathbf{r}_{\alpha}(t^{\prime})\ket{g}$
(133)
$\displaystyle\bigotimes_{\mathbf{k},\mu\simeq\mathbf{k}_{L,\mu}}e^{i[\varphi_{\mathbf{k}\mu}(t,t^{\prime},\mathbf{v}_{\alpha})+\varphi_{\mathbf{k}\mu}^{\chi,\delta}(t,t^{\prime})+\varphi_{\mathbf{k}\mu}^{\chi,\delta,\alpha}(t,t^{\prime})]}\ket{\left[\alpha_{\mathbf{k}_{L}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t}}$
$\displaystyle\bigotimes_{\mathbf{k},\mu\gg\mathbf{k}_{L,\mu}}e^{i[\varphi_{\mathbf{k}\mu}(t,t^{\prime},\mathbf{v}_{\alpha})+\varphi_{\mathbf{k}\mu}^{\chi,\delta}(t,t^{\prime})]}\ket{\left[\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t}},$

where

$\displaystyle\varphi_{\mathbf{k}\mu}^{\chi,\delta,\alpha}(t,t^{\prime})=\operatorname{Im}\left\\{\alpha_{\mathbf{k}_{L}}\left[\delta_{\alpha}^{*}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})+\chi_{\mathbf{k},\mu}^{*}(t^{\prime})\right]\right\\}.$
(134)

We have further used that

$\displaystyle
e^{-iH_{f}t}\mathbf{E}_{Q}(t^{\prime},\mathbf{R}_{\alpha})e^{iH_{f}t}=\mathbf{E}_{Q}(t^{\prime}-t,\mathbf{R}_{\alpha})=i\sum_{\mu}\int\differential^{3}k\tilde{g}(k)\epsilon_{k\mu}\left[a_{k\mu}e^{-i\omega_{k}(t^{\prime}-t)+i\mathbf{k}\cdot\mathbf{R}_{\alpha}}-\text{h.c.}\right].$
(135)

This state (133) is the field state conditioned on ATI with a single electron
momentum in the laboratory frame. We shall now project this state on the
harmonic vacuum $\ket{\\{0\\}_{HH}}$ in order to obtain the state of the
driving field mode in the situation where no harmonic radiation is generated.
We thus find that

$\displaystyle\ket{\Phi_{ATI}^{\mathbf{k}_{L}}(t,\mathbf{p})}=$
$\displaystyle\bra{\\{0\\}_{HH}}\ket{\Phi_{ATI}(t,\mathbf{p})}$ (136)
$\displaystyle=$ $\displaystyle\frac{-ie}{\hbar}\int_{t_{0}}^{t}\differential
t^{\prime}\left[\mathbf{E}_{Q,\mathbf{k}_{L}}(t^{\prime}-t,\mathbf{R}_{\alpha})+i\sum_{\mu}\int\differential^{3}k\Theta(\mathbf{k}_{HH})\tilde{g}(k)\epsilon_{\mathbf{k},\mu}e^{-i\omega_{k}t^{\prime}+i\mathbf{k}\cdot\mathbf{R}_{\alpha}}\left[\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]\right]$
$\displaystyle\left.-\mathbf{E}_{\delta}(t^{\prime},\mathbf{R}_{\alpha})-\mathbf{E}_{cl}(t^{\prime},\mathbf{R}_{\alpha})\right]\cdot\bra{\mathbf{v}_{\alpha}}\mathbf{r}_{\alpha}(t^{\prime})\ket{g}\prod_{\mathbf{k},\mu\gg\mathbf{k}_{L,\mu}}\Big{[}e^{i[\varphi_{\mathbf{k}\mu}(t,t^{\prime},\mathbf{v}_{\alpha})+\varphi_{\mathbf{k}\mu}^{\chi,\delta}(t,t^{\prime})]}e^{-\frac{1}{2}\absolutevalue{\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})}^{2}}\Big{]}$
$\displaystyle\bigotimes_{\mathbf{k},\mu\simeq\mathbf{k}_{L,\mu}}e^{i[\varphi_{\mathbf{k}\mu}(t,t^{\prime},\mathbf{v}_{\alpha})+\varphi_{\mathbf{k}\mu}^{\chi,\delta}(t,t^{\prime})+\varphi_{\mathbf{k}\mu}^{\chi,\delta,\alpha}(t,t^{\prime})]}\ket{\left[\alpha_{\mathbf{k}_{L}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t}},$

where $\Theta(\mathbf{k}_{HH})$ is a Heaviside function taking into account
only the contributions with field momentum $\mathbf{k}_{HH}\gg\mathbf{k}_{L}$.
Having in mind that the action of the electric field operator over a generic
coherent state $\ket{\beta_{\mathbf{k}}}$ can be written as

$\displaystyle\mathbf{E}_{Q,\mathbf{k}}(t,\mathbf{R}_{\alpha})\ket{\beta_{\mathbf{k}}}$
$\displaystyle=\mathbf{E}_{\beta.\mathbf{k}}(t,\mathbf{R}_{\alpha})\ket{\beta_{\mathbf{k}}}-i\tilde{g}(k)\epsilon_{\mathbf{k},\mu}D(\beta_{\mathbf{k}})\ket{1_{\mathbf{k}}},$
(137)

we can then write

$\displaystyle\mathbf{E}_{Q,\mathbf{k}_{L}}(t^{\prime}-t,\mathbf{R}_{\alpha})\ket{\left[\alpha_{\mathbf{k}_{L}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t}}$
(138) $\displaystyle\hskip
56.9055pt=\big{[}\mathbf{E}_{cl}(t^{\prime},\mathbf{R}_{\alpha})+\mathbf{E}_{\delta}(t^{\prime},\mathbf{R}_{\alpha})+\mathbf{E}_{\chi}(t^{\prime},\mathbf{R}_{\alpha})\big{]}\ket{\left[\alpha_{\mathbf{k}_{L}}+\xi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t}}$
$\displaystyle\hskip
56.9055pt\quad-i\tilde{g}(\mathbf{k}_{L})\epsilon_{\mathbf{k}_{L},\mu}D\big{(}\left[\alpha_{\mathbf{k}_{L}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t}\big{)}\ket{1}$

such that the conditioned to ATI state reads

$\displaystyle\ket{\Phi_{ATI}^{\mathbf{k}_{L}}(t,\mathbf{p})}$
$\displaystyle=\bra{\\{0\\}_{HH}}\ket{\Phi_{ATI}(t,\mathbf{p})}$ (139)
$\displaystyle=\dfrac{-ie}{\hbar}\int^{t}_{t_{0}}\differential
t^{\prime}\bigg{[}\mathbf{E}_{\chi}(t^{\prime},\mathbf{R}_{\alpha})+\int\differential^{3}k\Theta(\mathbf{k}_{HH})\tilde{g}(k)\epsilon_{\mathbf{k},\mu}e^{-i\omega_{k}t^{\prime}+i\mathbf{k}\cdot\mathbf{R}_{\alpha}}\left[\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]\bigg{]}$
$\displaystyle\quad\cdot\bra{\mathbf{v}_{\alpha}}\mathbf{r}_{\alpha}(t^{\prime})\ket{g}\mathcal{C}_{\text{HH}}(p,t,t^{\prime})\bigotimes_{\mathbf{k},\mu\simeq\mathbf{k}_{L,\mu}}e^{i\theta_{\mathbf{k},\mu}(t,t^{\prime},\mathbf{p},\alpha)}\ket{\left[\alpha_{\mathbf{k}_{L}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t}}$
$\displaystyle\quad-\dfrac{e}{\hbar}\int\differential^{3}k\Pi(\mathbf{k}_{L})\tilde{g}(\mathbf{k})\int^{t}_{t_{0}}\differential
t^{\prime}\bm{\epsilon}_{\mathbf{k_{L}},\mu}\cdot\bra{\mathbf{v}_{\alpha}}\mathbf{r}_{\alpha}(t^{\prime})\ket{g}\mathcal{C}_{\text{HH}}(p,t,t^{\prime})\quad$
$\displaystyle\quad\quad\times\prod_{\mathbf{k,\mu}\simeq\mathbf{k}_{L,\mu}}D\big{(}\left[\alpha_{\mathbf{k}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t}\big{)}\bigotimes_{\mathbf{k}^{\prime},\mu\simeq\mathbf{k}_{L,\mu}}e^{i\theta_{\mathbf{k}^{\prime},\mu}(t,t^{\prime},\mathbf{p},\alpha)}\ket{1_{\mathbf{k}},\\{0\\}_{\mathbf{k^{\prime}}\neq\mathbf{k}}},$

where $\Theta(\mathbf{k}_{\text{HH}})$ is a Heaviside function that takes into
account contributions for which $\mathbf{k}_{\text{HH}}\gg\mathbf{k}_{L}$,
$\Pi(\mathbf{k}_{L})$ is a rectangular function that is one whenever
$\mathbf{k}\simeq\mathbf{k}_{L}$, and zero otherwise. We have further defined

$\displaystyle\theta_{\mathbf{k},\mu}(t,t^{\prime},\mathbf{p},\alpha)$
$\displaystyle=\phi_{\mathbf{k},\mu}(t,t^{\prime},\mathbf{v}_{\alpha})+\operatorname{Im}\left[\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\chi_{\mathbf{k},\mu}^{*}(t^{\prime})\right]+\operatorname{Im}\left\\{\alpha_{\mathbf{k}_{L}}\left[\delta_{\alpha}^{*}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})+\chi_{\mathbf{k},\mu}^{*}(t^{\prime})\right]\right\\},$
(140)

which is the BCH phase that appears together with other phases arising from
applying the commutation relation between two displacement operators. Finally,
we have defined

$\displaystyle C_{\text{HH}}(p,t,t^{\prime})$
$\displaystyle=\prod_{\mathbf{k},\mu\gg\mathbf{k}_{L,\mu}}\Big{[}e^{i\big{\\{}\varphi_{\mathbf{k}\mu}(t,t^{\prime},\mathbf{v}_{\alpha})+\operatorname{Im}\left[\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\chi_{\mathbf{k},\mu}^{*}(t^{\prime})\right]\big{\\}}}e^{-\frac{1}{2}\absolutevalue{\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})}^{2}}\Big{]},$
(141)

which is a weight function that arises as a consequence of the conditioning
onto the vacuum state of the harmonics modes. We have also defined

$\displaystyle\mathbf{E}_{\chi}(t^{\prime},\mathbf{R}_{\alpha})$
$\displaystyle=i\sum_{\mu}\int\differential^{3}k\tilde{g}(k)\epsilon_{k\mu}\big{[}\chi_{\mathbf{k},\mu}(t^{\prime})e^{-i\omega_{k}t+i\mathbf{k}\cdot\mathbf{R}_{\alpha}}-\chi_{\mathbf{k},\mu}^{*}(t^{\prime})e^{i\omega_{k}t-i\mathbf{k}\mathbf{R}_{\alpha}}\big{]}.$
(142)

Finally, we project the previous expression with respect to the a Fock state
basis $\ket{n}$ to get the photon number probability amplitude for the direct
term. Note that, in order to get expressions that can be analyzed
computationally, we consider a single fundamental mode, corresponding to the
central frequency of the employed electric field

$\displaystyle\ket{\Phi_{ATI}^{\mathbf{k}_{L}}(t,\mathbf{p})}$
$\displaystyle=\bra{\\{0\\}_{HH}}\ket{\Phi_{ATI}(t,\mathbf{p})}$ (143)
$\displaystyle=\dfrac{-ie}{\hbar}\int^{t}_{t_{0}}\differential
t^{\prime}\bigg{[}\mathbf{E}_{\chi}(t^{\prime},\mathbf{R}_{\alpha})+\int
d^{3}k\Theta(\mathbf{k}_{HH})\tilde{g}(k)\epsilon_{\mathbf{k},\mu}e^{-i\omega_{k}t^{\prime}+i\mathbf{k}\cdot\mathbf{R}_{\alpha}}\left[\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]\bigg{]}$
$\displaystyle\quad\cdot\bra{\mathbf{v}_{\alpha}}\mathbf{r}_{\alpha}(t^{\prime})\ket{g}\mathcal{C}_{\text{HH}}(p,t,t^{\prime})e^{i\theta_{\mathbf{k},\mu}(t,t^{\prime},\mathbf{p},\alpha)}e^{-\tfrac{1}{2}\lvert\alpha_{\mathbf{k}_{L}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\rvert^{2}}$
$\displaystyle\quad\times\dfrac{(\left[\alpha_{\mathbf{k}_{L}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t})^{n}}{\sqrt{n!}}$
$\displaystyle\quad-\dfrac{e}{\hbar}\int\differential^{3}k\Pi(\mathbf{k}_{L})\tilde{g}(\mathbf{k})\int^{t}_{t_{0}}\differential
t^{\prime}\bm{\epsilon}_{\mathbf{k_{L}},\mu}\cdot\bra{\mathbf{v}_{\alpha}}\mathbf{r}_{\alpha}(t^{\prime})\ket{g}\mathcal{C}_{\text{HH}}(p,t,t^{\prime})\quad$
$\displaystyle\quad\quad\times
e^{i\theta_{\mathbf{k}^{\prime},\mu}(t,t^{\prime},\mathbf{p},\alpha)}e^{-\tfrac{1}{2}\lvert\alpha_{\mathbf{k}_{L}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\rvert^{2}}$
$\displaystyle\quad\Big{[}\sqrt{n}\dfrac{(\left[\alpha_{\mathbf{k}_{L}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t})^{n-1}}{\sqrt{n-1!}}-\dfrac{(\left[\alpha_{\mathbf{k}_{L}}+\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(t,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\right]e^{-i\omega_{\mathbf{k}}t})^{n}}{\sqrt{n!}}\Big{]},$

such that, up to a normalization factor, the photon number probability
distribution for the $\mathbf{k}_{L,\mu}$ mode is given by

$P_{n_{\mathbf{k}_{L}}}(t,\mathbf{p})=\Big{\lvert}\innerproduct{n}{\Phi_{ATI}^{\mathbf{k}_{L}}(t,\mathbf{p})}\Big{\rvert}^{2}.$
(144)

##  Appendix F Perturbation theory for rescattering in ATI

From the perturbation theory analysis considered in the main part of the
manuscript (79) for the introduction of the rescattering part of the state, we
get the following set of equations for each of the perturbation orders

$\displaystyle
i\hbar\derivative{t}\ket{\phi^{(0)}(\mathbf{v},t)}=e\mathbf{E}_{Q}(t,\mathbf{R})\cdot\matrixelement{\mathbf{v}}{\mathbf{r}(t)}{g}\ket{\phi(t)}+e\mathbf{E}_{Q}(t,\mathbf{R})\cdot\Delta\mathbf{r}(t,\mathbf{v}_{\alpha})\ket{\phi^{(0)}(\mathbf{v},t)}$
(145) $\displaystyle
i\hbar\derivative{t}\ket{\phi^{(1)}(\mathbf{v},t)}=e\mathbf{E}_{Q}(t,\mathbf{R})\Delta\mathbf{r}(t,\mathbf{v})\ket{\phi^{(1)}(\mathbf{v},t)}-\dfrac{\hbar}{e}\int\differential^{3}v^{\prime}\exp[\frac{i}{\hbar}S(\mathbf{v}^{\prime},t,t_{0})]\mathbf{E}_{Q}(t,\mathbf{R})\cdot\mathbf{g}(\mathbf{v},\mathbf{v}^{\prime})$
$\displaystyle\hskip
85.35826pt\times\exp[-\frac{i}{\hbar}S(\mathbf{v}^{\prime},t^{\prime},t_{0})]\ket{\phi^{(0)}(\mathbf{v}^{\prime},t)}.$

Here, we define $\Delta\mathbf{r}(t,\mathbf{v})$ as in Eq. (62), so that the
first equation describes the dynamics of the zeroth order perturbation term
which characterizes direct ionization phenomena, and which has been previously
discussed. On the other hand, the second term describes the rescattering
dynamics, i.e., electrons that, once ionized, and can then re-scatter at the
potential of the parent ion. Similarly to what we had in the direct ionization
term, the solution to this differential equation is given as the superposition
of an homogeneous term similar to that in Eq. (63), plus an inhomogeneous term
of the form

$\displaystyle\ket{\phi^{(1)}(\mathbf{v},t)}=$
$\displaystyle-i\int^{t}_{t_{0}}\differential
t_{2}\int\differential^{3}v^{\prime}\prod_{\mathbf{k},\mu}e^{i\varphi_{\mathbf{k}\mu}(t,t_{2},\mathbf{v})}D\big{(}\delta(t,t_{2},\omega_{k},\mathbf{v})\big{)}\exp[\frac{i}{\hbar}S(\mathbf{v},t,t_{0})]$
(146)
$\displaystyle\times\mathbf{E}_{Q}(t_{2},\mathbf{R})\cdot\mathbf{g}(\mathbf{v},\mathbf{v}^{\prime})\exp[-\frac{i}{\hbar}S(t_{2},t_{0},\mathbf{v}^{\prime})]\ket{\phi^{(0)}(\mathbf{v}^{\prime},t_{2})}.$

However, since the electron is initially in the ground state so that
$\ket{\phi^{(1)}(\mathbf{v},t_{0})}=0$ which implies that the homogeneous
solution does not participate in the final solution. Then, the rescattering
part of the state is entirely described by Eq. (146).

##  Appendix G Description of the numerical implementation

In the following we describe the numerical procedure used to obtain the field
observables in the main part of the manuscript (see Section IV.3.2). The
analysis is implemented in Python, and the codes can be found in [183].

Our numerical implementation consists of the following steps:

  1. 1.

Once the vector potential $\mathbf{A}(t)$ of the pulse under consideration is
defined, we compute the expectation value of the time-dependent dipole moment
$\langle\mathbf{d}(t)\rangle$ under the SFA framework.

  2. 2.

Given a time interval $[t_{0},T]$, where $t_{0}$ and $T$ corresponds to the
beginning and end of the pulse respectively, we compute $\chi(t)$ and
$C_{\text{HH}}(p,T,t)$ given by Eq. (141). We use a step size of $\Delta
t=1.0$. For the displacement via $\delta(t,t_{0},\omega_{k},p)$ we compute the
analytic exact expression for a $\sin^{2}$ envelope pulse.

  3. 3.

Finally, we use an interpolating scheme such that we can have a numerical
function for $\chi(t)$ and $C_{\text{HH}}(p,T,t)$, which can adapt to the time
integrals that need to be done for computing the photon-number probability
distribution.

Furthermore, along this section we will work with atomic units (a.u.), such
that $\hbar=1,e^{2}=1,m=1$ and $\epsilon_{0}=4\pi$.

###  G.1 SFA expressions for the time-dependent dipole

The expectation value of the time-dependent dipole moment is defined as

$\langle\mathbf{d}(t)\rangle=\matrixelement{\psi(t)}{\mathbf{d}}{\psi(t)},$
(147)

where $\ket{\psi(t)}$ is the wavefunction of a single electron in the intense
laser field. To obtain the wavefunction we introduce the SFA approximations
[118], namely (i) the contribution of all bound states except the ground state
are neglected; (ii) the depletion of the ground state can be neglected; and
(iii) in the continuum, the electron can be treated as a free particle, that
is, we can neglect the interaction with the binding potential. Thus, we can
write the semi-classical quantum state describing the electron as

$\ket{\psi(t)}=e^{iI_{p}t}\ket{g}+\int\differential^{3}p\
b(\mathbf{p},t)\ket{\mathbf{p}+\mathbf{A}(t)},$ (148)

where the probability amplitude of finding the electron in the continuum state
$\ket{\mathbf{p}+\mathbf{A}(t)}$ is given by

$\displaystyle b(\mathbf{p},t)$ $\displaystyle=i\int^{t}_{t_{0}}\differential
t^{\prime}\mathbf{E}(t^{\prime})\cdot\mathbf{d}(\mathbf{p}+\mathbf{A}(t^{\prime}))\exp{-i\int^{t}_{t^{\prime}}\differential\tau\Big{[}\big{(}\mathbf{p}+\mathbf{A}(\tau)\big{)}^{2}+I_{p}\Big{]}},$
(149)

where
$\mathbf{d}(\mathbf{p}+\mathbf{A}(t))=\bra{\mathbf{p}+\mathbf{A}(t)}\mathbf{d}\ket{g}$
is the transition dipole moment between the ground and the continuum state.
Introducing these expressions in (147), and neglecting the contributions from
the continuum-continuum transitions, we get the SFA expression for the dipole
moment expectation value

$\displaystyle\langle\mathbf{d}(t)\rangle$
$\displaystyle=i\int\differential^{3}p\int^{t}_{t_{0}}\differential
t^{\prime}\mathbf{d}^{*}(\mathbf{p}+\mathbf{A}(t))\mathbf{E}(t^{\prime})\cdot\mathbf{d}(\mathbf{p}+\mathbf{A}(t^{\prime}))\exp{-i\int^{t}_{t^{\prime}}\differential\tau\Big{[}\big{(}\mathbf{p}+\mathbf{A}(\tau)\big{)}^{2}+I_{p}\Big{]}}+\text{c.c.}$
(150)

For the numerical analysis, we considered a linearly polarized field with
vector potential given by

$\mathbf{A}(t)=\epsilon_{z}\dfrac{E_{0}}{\omega}\sin[2](\dfrac{\omega}{2n_{\text{cyc}}})\sin(\omega
t),$ (151)

where $n_{\text{cyc}}$ is the total number of cycles of the pulse with a
$\sin^{2}$ envelope, $\epsilon_{z}$ is a unitary vector along the $z$ axis
describing the polarization of the state. In all numerical analysis, we use
$\omega=0.057$ a.u., and $E_{0}=0.053$ a.u. On the other hand, we consider the
ground state wavefunction of the electron to be described by

$\innerproduct{x}{g}=\sqrt{\dfrac{\lambda^{3}}{\pi}}e^{-\lambda r},$ (152)

such that the dipole transition matrix element connecting the continuum state
with the ground state is given by

$\displaystyle\mathbf{d}(\mathbf{p}+\mathbf{A}(t))$
$\displaystyle=\bra{\mathbf{p}+\mathbf{A}(t))}\mathbf{d}\ket{g}=-i\sqrt{\dfrac{\lambda^{3}}{\pi}}\dfrac{\big{(}\mathbf{p}+\mathbf{A}(t)\big{)}}{(2\pi)^{3/2}}\dfrac{32\pi\lambda}{(\lambda^{2}+(p+A(t))^{2})^{3}}.$
(153)

Figure 18: (a) Dipole moment expectation value obtained from the numerical
calculations under the SFA approach. (b) Vector potential used to obtain the
dipole moment shown in (a). The pulse shown has $n_{\text{cyc}}=5$ cycles.

More particularly, for a 1s state of an hydrogen atom we set $I_{p}=0.5$ a.u.,
and $\lambda=1$. In Fig. 18 we show in (a) the dipole moment expectation value
(150) obtained, and in (b) the corresponding vector potential with
$n_{\text{cyc}}=5$. In order to perform the two integrals, we have first
discretized the one for the momentum (which we perform in 1D), and used the
quad function from the SciPy package [184].

###  G.2 Calculation and interpolation of the quantum optical displacements

With the dipole moment obtained in the previous subsection within the SFA
approach, we can now compute the photonic displacement given by $\xi(t)$,
which is defined as

$\chi(t)\propto\int^{t}_{t_{0}}\ \differential t^{\prime}e^{i\omega
t^{\prime}}\langle\mathbf{d}(t^{\prime})\rangle.$ (154)

We further compute

$\displaystyle C_{\text{HH}}(p,T,t^{\prime})$
$\displaystyle=\prod_{\mathbf{k},\mu\gg\mathbf{k}_{L,\mu}}\Big{[}e^{i\big{\\{}\varphi_{\mathbf{k}\mu}(T,t^{\prime},\mathbf{v}_{\alpha})+\operatorname{Im}\left[\delta_{\alpha}(T,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})\chi_{\mathbf{k},\mu}^{*}(t^{\prime})\right]\big{\\}}}e^{-\frac{1}{2}\absolutevalue{\chi_{\mathbf{k},\mu}(t^{\prime})+\delta_{\alpha}(T,t^{\prime},\omega_{k},\mathbf{v}_{\alpha})}^{2}}\Big{]},$
(155)

where $\delta(p,T,t^{\prime},\omega_{k})$ is obtained analytically. In
particular, the computation of $C_{\text{HH}}(p,T,t^{\prime})$ involves the
calculation of the BCH phase
$\varphi_{\mathbf{k}\mu}(T,t^{\prime},\mathbf{v}_{\alpha})$ from (131). For
this, we consider a discrete set of harmonic modes consisting of the 2nd up to
the 21st harmonic, and perform an interpolation for each of the BCH phase
factors. In principle we could consider more harmonic orders, but their
amplitude $\absolutevalue{\delta(t,t^{\prime},\omega_{k},\mathbf{v})}$ rapidly
decrease for increasing harmonic orders, and therefore their contributions are
negligible. Since the effect of $\chi_{\mathbf{k},\mu}$ can be neglected in
comparison to that introduced by $\delta(t,t^{\prime},\omega_{k},\mathbf{v})$,
we neglect their effect in the previous conditioning operation. Finally, once
the $C_{\text{HH}}(p,T,t^{\prime})$ has been obtained for each value of the
ionization time $t^{\prime}$, we proceed to use an interpolation scheme to
compute the photon number probability distribution.

In both cases, similarly to what we did before, we perform this integral with
the quad function, and is followed by performing an interpolation using the
interp1d function provided by SciPy. The reason behind this is that, after the
interpolation, we get two functions that exactly behave as $\chi(t)$ and
$C_{\text{HH}}(p,T,t^{\prime})$ within the $[t_{0},T]$ interval. Thus, by
doing this we avoid the nested integrals that appear within the definition of
the photon-number probability distribution. This provides a huge numerical
speed up. In particular, for doing the interpolation, we used the cubic option
which considers polynomials of zeroth to third order, and that better fits the
oscillating behaviour of both functions. The results of the interpolation are
shown with the dashed orange curves in Fig. 19.

Figure 19: Behaviour of the (a) $\chi_{\mathbf{k}_{L}}(t)$ and (b)
$C_{\text{HH}}(p=2U_{p},T,t)$ with a time interval of $\Delta t=1.0$ (solid
black curve), and the function obtained from the interpolation (dashed orange
curve). For the interpolation plot, we considered a time step $\Delta t=0.01$.

For computing the dipole moment expectation value, we defined the time step to
be equal to $\Delta t=1.0$. In Fig. 19, we show the behaviour of $\chi(t)$ and
$\delta(p,T,t)$ (with $p=0.1U_{p}$) defined by the integral (solid black
line), and after the interpolation (cyan dashed line). For the interpolation
plot, we considered a time step $\Delta t=0.01$.

###  G.3 Further considerations

As written in the main part of the manuscript, the photon number probability
distribution is computed in the original reference frame, i.e. after undoing
the displacement operation $D(\alpha)$. The usual regime of intense laser
field physics has coherent states amplitudes of the driving laser in the order
of $|\alpha|\sim 10^{6}$. Such high photon numbers can only be hardly
implemented in the numerical analysis due to the numerical accuracy. Thus, in
order to avoid this problem, we artificially reduce the value of $\alpha$ and
multiply the $g(k)$ by an extra factor $\lambda<1$ to account for this
reduction in the computed values of $\delta(p,t,t^{\prime})$ and $\xi(t)$.

On the other hand, the specific value of $\alpha$ that we used in our
numerical analysis is $\alpha=7i$. This value is obtained from the definition
of the vector potential, which is proportional to

$\hat{A}(t)\propto f(t)\big{[}\hat{a}^{\dagger}e^{i\omega
t}+\hat{a}e^{-i\omega t}\big{]}.$ (156)

Note that in the last expression we have introduced the envelope function via
$f(t)$, which appears as a consequence of the multimode nature of our quantum
optical description. If we now compute the average value of the vector
potential operator with respect to the selected coherent state, we get

$\langle A(t)\rangle=\matrixelement{7i}{\hat{A}(t)}{7i}\propto
14f(t)\sin(\omega t),$ (157)

which has the form of the vector potential we used in the previous part of our
numerical analysis. Note that one can similarly get the same phase by looking
at the real and imaginary part of $\alpha_{\mathbf{k},\mu}$ in Fig. 3 of the
main text.