metadata

# The imaginary part of the high-harmonic cutoff

Emilio Pisanty<EMAIL_ADDRESS>ICFO – Institut de Ciencies Fotoniques,
The Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona)  Marcelo F. Ciappina  ICFO – Institut de Ciencies Fotoniques, The
Barcelona Institute of Science and Technology, 08860 Castelldefels (Barcelona)
Maciej Lewenstein  ICFO – Institut de Ciencies Fotoniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona)  ICREA,
Passeig de Lluís Companys, 23, 08010 Barcelona, Spain

(22 July 2020)

###### Abstract

High-harmonic generation – the emission of high-frequency radiation by the
ionization and subsequent recombination of an atomic electron driven by a
strong laser field – is widely understood using a quasiclassical trajectory
formalism, derived from a saddle-point approximation, where each saddle
corresponds to a complex-valued trajectory whose recombination contributes to
the harmonic emission. However, the classification of these saddle points into
individual quantum orbits remains a high-friction part of the formalism. Here
we present a scheme to classify these trajectories, based on a natural
identification of the (complex) time that corresponds to the harmonic cutoff.
This identification also provides a natural complex value for the cutoff
energy, whose imaginary part controls the strength of quantum-path
interference between the quantum orbits that meet at the cutoff. Our
construction gives an efficient method to evaluate the location and brightness
of the cutoff for a wide class of driver waveforms by solving a single saddle-
point equation. It also allows us to explore the intricate topologies of the
Riemann surfaces formed by the quantum orbits induced by nontrivial waveforms.

Accepted Manuscript for J. Phys. Photonics 2, 034013 (2020), available as
arXiv:2003.00277.

†† © The authors, 2020. Licensed under CC BY 4.0.

High-harmonic generation (HHG) is an extremely nonlinear optical process in
which a strong laser field drives the emission of a train of short bursts of
high-frequency radiation Krausz2009 ; Corkum2007 , which can cover hundreds of
harmonic orders of the driving field, over a broad plateau that terminates at
a cutoff. This emission comes from a three-step process in which the laser
ionizes the target atom via tunnel ionization, and then propels the released
electron back to the parent ion, where it recombines with the hole it left
behind, releasing its kinetic energy as a photon Corkum1993 ; Kulander1993 .

This emission can be modelled using a wide range of approaches, from classical
heuristics Corkum1993 ; Kulander1993  to intensive numerical computations
Scrinzi2014 , but the quantitative models that most closely follow the overall
intuition are quasiclassical methods Lewenstein1994 ; Amini2019 , known as the
Strong-Field Approximation (SFA), where the emission amplitude is given by a
path-integral sum over discrete emission events. These are known as quantum
orbits Salieres2001 ; Paulus2000 ; Kopold2002 , i.e., quasiclassical
trajectories whose start and end times are complex Ivanov2014 ; Nayak2019 .

The quantum-orbit formalism arises naturally under the approximation that the
electron’s motion in the continuum is completely controlled by the driving
laser, which allows an exact solution in terms of highly oscillatory
integrals. These are then reduced, using a steepest-descent method known as
the saddle-point approximation (SPA), to discrete contributions coming from
the saddle points of the integrand BruijnAsymptotics ; BleisteinIntegrals ;
GerlachSPAonline . These saddle points represent the quasiclassical
trajectories, and they typically come in pairs – most notably the ‘short’ and
‘long’ trajectories Lewenstein1995 ; Zair2008 . Each pair of trajectories
approaches each other over the harmonic plateau and then performs a Stokes
transition at the harmonic cutoff Figueira2002 ; Milosevic2002 , giving way to
an exponential drop where only one of the saddles is used. In practice,
however, classifying the saddle points into these pairs of trajectories is one
of the highest-friction points when applying this method Chipperfield2007 ;
Hoffmann2011 ; Das2017 , particularly since the saddles tend to move quickly,
and approach each other very closely, at the harmonic cutoff.

---

(a)

(b)

(c)

Figure 1:  Solving the HHG saddle-point equations returns a discrete set of
saddle points as a function of the harmonic photon energy $\Omega$, shown in
(b) for the monochromatic field in (a). Our main result is the trajectory
classification in (c), which consists of the harmonic-cutoff points
$t_{\mathrm{hc}}$ (triangles) and the separatrices through them, which allow
the cloud of saddle-point solutions in (b) to be organized into individual
trajectories; the different colors correspond to different quantum orbits. We
model neon in a field of wavelength $800\text{\,}\mathrm{n}\mathrm{m}$ and
intensity
$2\text{\times}{10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$.

In this work we construct a quantum-orbit classification scheme based on a
natural notion of complex-valued harmonic-cutoff times $t_{\mathrm{hc}}$.
These are points – given by zeros of the second derivative of the action –
which sit between the two saddle points as they approach each other, and which
provide an organic separation between the two. Thus, an unordered set of
saddle-point solutions like the one shown in Fig. 1b can be cleanly organized
programmatically into families of trajectories, as shown in Fig. 1c, in a
flexible and general fashion which is robust to changes in the quantum orbits
over broad families of optical fields.

Once these second-order saddle points have been identified, they naturally
fill in the role of the quantum orbits corresponding to the high-harmonic
cutoff, and they give the photon energy at which it occurs – a quantity which
can often be hard to pin down precisely – as the real part of the time
derivative of the action at the harmonic-cutoff time. We benchmark this
understanding of the harmonic cutoff against the standard ‘cutoff law’, which
describes the scaling as $\Omega_{\mathrm{max}}\sim I_{p}+3.17U_{p}$ with the
ionization potential $I_{p}$ and the ponderomotive energy $U_{p}$ of the field
(previously derived both classically Corkum1993 ; Kulander1993  and via a
systematic expansion in $I_{p}/U_{p}$ Lewenstein1994 ). However, our method
extends trivially to drivers with higher harmonic content as well as to
tailored polarizations.

Moreover, the information at $t_{\mathrm{hc}}$, which can be obtained by
solving a single second-order saddle point equation, is sufficient to
calculate an accurate evaluation of the harmonic yield at the cutoff, as well
as a good qualitative estimate (which we term the Harmonic-Cutoff
Approximation) for the shape of the spectrum at the upper plateau. The
efficiency of this approach makes it a good tool when optimizing both the
position and the brightness of the cutoff over the high-dimensional parameter
spaces available on current optical synthesizers Manzoni2015 .

This understanding of the high-harmonic cutoff $\Omega_{\mathrm{hc}}$ also
assigns a natural value to its imaginary part, whose direct impact is to
control the closeness of the approach between the pair of saddle points at the
cutoff. Since this closeness regulates, in turn, the distinguishability
between the two trajectories, the imaginary part of the harmonic cutoff
ultimately controls the strength of quantum-path interference (QPI) Zair2008
between the pair of orbits.

In particular, the zeros of the imaginary part of the harmonic cutoff energy
$\imaginary(\Omega_{\mathrm{hc}})$ pinpoint the configurations where the
saddle points for the two quantum orbits have an exact coalescence at the
cutoff, in which case they cannot be distinguished from each other. For
tailored polarizations, as well as other polychromatic fields with one or more
nontrivial shape parameters, $\imaginary(\Omega_{\mathrm{hc}})$ will
generically oscillate, indicating that the quantum orbits have reconnection
events – where, in essence, the evanescent post-cutoff saddle point will
transfer from the ‘short’ quantum orbit to the ‘long’ one.

We showcase this behaviour in bichromatic counter-rotating circularly-
polarized ‘bicircular’ fields Fleischer2014 ; Kfir2015 ; Eichmann1995 ;
Milosevic1996 ; Milosevic2000bicircular ; Baykusheva2016 ; Dorney2017 ;
JimenezGalan2018 ; Pisanty2014 ; Pisanty2017 ; Milovsevic2019  as the relative
intensity between the two components changes: the various quantum orbits then
recombine with each other, making for a complicated landscape within a unified
topology. This illustrates the fact that the various saddle points are simply
different branches of a single, unified Riemann surface, with our harmonic-
cutoff times $t_{\mathrm{hc}}$ taking on the role of the branch points of this
Riemann surface. More practically, the reconnections make the quantum-orbit
classification challenging, but we show that our algorithm can seamlessly
handle these changes in the trajectories.

In the more structural terms of catastrophe theory Poston1978 , the harmonic
cutoff is a spectral caustic, in the ‘fold’ class of diffraction catastrophes
Raz2012 . In this sense, the complex harmonic-cutoff energy
$\Omega_{\mathrm{hc}}$ is the bifurcation set for the catastrophe – suitably
generalized to the complex variables involved – and it marks the location of
the caustic. The formal study of caustics has seen increasing interest within
attoscience Raz2012 ; Goulielmakis2012 ; Austin2014 ; Facciala2016 ;
Facciala2018 ; Hamilton2017 ; Birulia2019 ; Uzan2020 ; Kelvich2017 , and our
approach provides a useful tool for exploring and describing the higher-
dimensional bifurcation sets at the heart of the higher-order catastrophes
that become available as our control over short pulses of light becomes more
finely tuned.

This work is structured as follows. In Section I we construct the harmonic-
cutoff times and examine their structure, first summarizing the standard SFA
in Section I.1 and constructing a model with a cubic action in Section I.2,
which we explore in depth in Sections I.3 and I.4, where we construct the
classification algorithm; in Section I.5 we explore the quantum-orbit Riemann
surface uncovered by this perspective. In Section II we construct the
Harmonic-Cutoff Approximation, by explicitly integrating the Airy integral of
the cubic action of our model, and we benchmark it against the standard
approximations. In Section III we examine how the complex cutoff energy
$\Omega_{\mathrm{hc}}$ scales with the field parameters, and show that it
agrees with the known cutoff law and that it extends it to the complex domain.
Finally, in Section IV, we explore the branch-cut classification for
bicircular fields, showing that our algorithm can handle the quantum-orbit
reconnections, as well as how these reconnections give rise to a nontrivial
topology for the quantum orbits.

The functionality we describe here has been integrated into the RBSFA package
for Mathematica RBSFA , and our specific implementation is available from Ref.
FigureMaker, . Interactive versions of 3D Fig. 1 and 3D Fig. 2, as well as
3D-printable models of those surfaces, are available as Supplementary Material
SupplementaryMaterial .

##  I The complex harmonic-cutoff times

###  I.1 The Strong-Field Approximation

In the SFA, high-harmonic emission is calculated as the expectation value of
the dipole operator, under the assumption that there is a single active
electron that is either in the atomic potential’s ground state $\ket{g}$ or in
a laser-driven continuum where the effect of the atomic potential is
negligible. The calculation Ivanov2014 ; Nayak2019 ; Amini2019 ;
Lewenstein1994  then gives the harmonic dipole in the form

$\displaystyle\mathbf{D}(\Omega)$
$\displaystyle=\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!\\!\textrm{d}t\int_{-\infty}^{t}\\!\\!\\!\\!\\!\\!\textrm{d}t^{\prime}\\!\int\\!\textrm{d}\mathbf{p}\,\mathbf{d}(\mathbf{p}+\mathbf{A}(t))\Upsilon(\mathbf{p}+\mathbf{A}(t^{\prime}))$
$\displaystyle\qquad\qquad\qquad\qquad\times
e^{-iS_{V}\\!(\mathbf{p},t,t^{\prime})+i\Omega t}+\mathrm{c.c.},$ (1)

where the three integrals over the times of ionization and recollision,
$t^{\prime}$ and $t$, and the intermediate momentum, $\mathbf{p}$, form a
reduced Feynman path integral summing over a restricted set of relevant orbits
Salieres2001 . This Feynman sum is modulated by ionization and recollision
amplitudes given by the transition dipole matrix element
$\mathbf{d}(\mathbf{k})=\matrixelement{\mathbf{k}}{\hat{r}}{g}$ and the scalar
function
$\Upsilon(\mathbf{k})=(I_{p}+\tfrac{1}{2}\mathbf{k}^{2})\innerproduct{\mathbf{k}}{g}$,
evaluated at the ionization and recollision velocities: $\mathbf{A}(t)$ is the
vector potential of the laser and $\mathbf{p}$ is the canonical momentum, so
$\mathbf{v}(t)=\mathbf{p}+\mathbf{A}(t)$ is the recollision velocity and
$m_{e}\mathbf{v}(t)=\mathbf{v}(t)$ is the kinematic momentum. (We use atomic
units with $m_{e}=\hbar=1$ throughout, and we consider neon in a linearly-
polarized field $\mathbf{F}(t)=F_{0}\hat{\mathbf{e}}_{z}\cos(\omega t)$ of
wavelength $800\text{\,}\mathrm{n}\mathrm{m}$ and intensity
$2\text{\times}{10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$
unless otherwise stated.)

More importantly, the contribution from each orbit in the Feynman sum in (1)
is modulated by a phase given by the total action accumulated during
propagation,

$\displaystyle S(\mathbf{p},t,t^{\prime})$
$\displaystyle=S_{V}(\mathbf{p},t,t^{\prime})-\Omega t$ (2a)
$\displaystyle\text{for}\ S_{V}(\mathbf{p},t,t^{\prime})$
$\displaystyle=\frac{1}{2}\int_{t^{\prime}}^{t}\left(\mathbf{p}+\mathbf{A}(\tau)\right)^{2}\textrm{d}\tau+I_{p}(t-t^{\prime}),$
(2b)

with $I_{p}$ the ionization potential of the atom, which is normally dominated
by the Volkov component $S_{V}$. For a field with amplitude $F_{0}$ and
frequency $\omega$, the action scales with the ratio $z=U_{p}/\omega$ of the
field’s ponderomotive energy $U_{p}=F_{0}^{2}/4\omega^{2}$ to its frequency.
This is typically large, so the exponential term $e^{-iS}$ changes much faster
than the amplitudes in the prefactor, making the integral in (1) highly
oscillatory.

The highly oscillatory nature of this amplitude then allows us to deal with
this integral using the method of steepest descents BruijnAsymptotics ;
BleisteinIntegrals ; GerlachSPAonline , also known as the saddle-point
approximation (SPA). In this method, we deform the integration contours in (1)
into the complex plane, so that they will pass through the saddle points of
the exponent (2). There the exponent is locally quadratic, so the integral can
be reformulated as a gaussian and integrated exactly, under the assumption
that the prefactor is slow. These points can be found through the saddle-point
equations

$\displaystyle 0$ $\displaystyle=\frac{\partial S}{\partial
t}=\frac{1}{2}\left(\mathbf{p}+\mathbf{A}(t)\right)^{2}+I_{p}-\Omega,$ (3a)
$\displaystyle 0$ $\displaystyle=\frac{\partial S}{\partial
t^{\prime}}=\frac{1}{2}\left(\mathbf{p}+\mathbf{A}(t^{\prime})\right)^{2}+I_{p},$
(3b) $\displaystyle 0$ $\displaystyle=\frac{\partial
S}{\partial\mathbf{p}}=\int_{t^{\prime}}^{t}\left(\mathbf{p}+\mathbf{A}(\tau)\right)\textrm{d}\tau,$
(3c)

which can be interpreted on physical grounds as encoding the requirements of
energy conservation at recollision and ionization ((3a) and (3b), resp.) as
well as the return of the quasiclassical laser-driven trajectory
$\boldsymbol{\alpha}(t)=\int_{t^{\prime}}^{t}\mathbf{v}(\tau)\textrm{d}\tau$
to the ion at the time of recollision.111We use the term ‘quasiclassical’ here
to distinguish this formalism from more general semiclassical approaches,
which use Feynman sums over allowed classical trajectories and add quantum
corrections coming from the gaussian (or other similar) spread of the integral
around them.  One key feature of these conditions is that the ionization
equation (3b), in particular, cannot be satisfied with real variables, forcing
all of the variables involved to take complex values.

Normally, the return equation (3c) is solved separately, since its linearity
in $\mathbf{p}$ guarantees a unique solution,

$\mathbf{p}_{s}(t,t^{\prime})=-\frac{1}{t-t^{\prime}}\int_{t^{\prime}}^{t}\mathbf{A}(\tau)\textrm{d}\tau,$
(4)

for any arbitrary pair $(t,t^{\prime})$ of ionization and recollision times.
Once the momentum integral has been performed in this way, the expression for
the harmonic dipole takes the form

$\displaystyle\mathbf{D}(\Omega)$
$\displaystyle=\int_{-\infty}^{\infty}\\!\\!\\!\\!\\!\textrm{d}t\int_{-\infty}^{t}\\!\\!\\!\\!\\!\\!\textrm{d}t^{\prime}\>\mathbf{d}(\mathbf{p}_{s}(t,t^{\prime}){+}\mathbf{A}(t))\Upsilon(\mathbf{p}_{s}(t,t^{\prime}){+}\mathbf{A}(t^{\prime}))$
$\displaystyle\qquad\times\left(\frac{2\pi}{i(t-t^{\prime})}\right)^{3/2}e^{-iS_{V}\\!(t,t^{\prime})+i\Omega
t},$ (5)

where the fractional power of the excursion time $\tau=t-t^{\prime}$
represents the dispersion of the released wavepacket in position space.222If
the time integrals are performed numerically, this factor needs to be
regularized at $\tau\to 0$, to account for a breakdown of the saddle-point
approximation in that limit Pisanty2016 . If the time integrals are also
evaluated via saddle-point methods, however, this is not required, as that
limit is not used. We notate
$S(t,t^{\prime})=S(\mathbf{p}_{s}(t,t^{\prime}),t,t^{\prime})$ where it does
not lead to confusion, and we drop the added complex conjugate for simplicity.

The resulting two-dimensional integral, (5), is now in its minimal form; the
saddle-point equations for the amended action read

$\displaystyle 0$ $\displaystyle=\frac{\partial S}{\partial
t}=\frac{1}{2}\left(\mathbf{p}_{s}(t,t^{\prime})+\mathbf{A}(t)\right)^{2}+I_{p}-\Omega,$
(6a) $\displaystyle 0$ $\displaystyle=\frac{\partial S}{\partial
t^{\prime}}=\frac{1}{2}\left(\mathbf{p}_{s}(t,t^{\prime})+\mathbf{A}(t^{\prime})\right)^{2}+I_{p},$
(6b)

and they can only be solved numerically.333This is often done by gradient
descent on the modulus of the right-hand side Ivanov2014 ; Nayak2019 , but it
is also possible to use Newton’s method directly, as it readily extends to
multiple complex variables Chipperfield2007 ; RBSFA .  A typical set of
solutions for these saddle-point equations is shown in Fig. 1b: the solutions
form a discrete set of points, which shift when the harmonic frequency
$\Omega$ changes. These discrete points thus trace out individual curves on
the complex $t$ and $t^{\prime}$ planes, which form the individual quantum
orbits.

---

(a)

(b)

(c)

Figure 2:  (a,b) Detail of the first pair of quantum orbits from Fig. 1,
labelled by harmonic order, $\Omega/\omega$. At the avoided crossing of the
two saddles, they go through a Stokes and then an anti-Stokes transition
(solid and hollow arrow, resp.), at the points where $\real(S)$ and
$\imaginary(S)$ of the two are equal, resp. (c) Saddle-point approximation to
the core elements of the harmonic dipole, $\tau^{-3/2}e^{-iS}$, for the first
pair of quantum orbits (green and blue curves). At the anti-Stokes transition
(dot-dashed line), one of the two starts growing exponentially, but it must be
discarded at the Stokes transition (solid line), when the required
integration-contour topology changes. This produces a discontinuous change in
the total SPA dipole (dashed line), which can be fixed by using the Uniform
Approximation (solid line). The oscillations over the plateau are quantum path
interference between the two trajectories.

Once the saddle points have been found, the SPA gives the harmonic dipole as

$\displaystyle\mathbf{D}(\Omega)$
$\displaystyle=\sum_{s}\frac{2\pi}{\sqrt{-\det[S^{\prime\prime}(t_{s},t^{\prime}_{s})]}}\left(\frac{2\pi}{i(t_{s}-t^{\prime}_{s})}\right)^{3/2}$
$\displaystyle\qquad\qquad\times\mathbf{d}(\mathbf{p}_{s}{+}\mathbf{A}(t_{s}))\Upsilon(\mathbf{p}_{s}{+}\mathbf{A}(t^{\prime}_{s}))\,$
(7) $\displaystyle\qquad\qquad\times
e^{-iS_{V}\\!(t_{s},t^{\prime}_{s})+i\Omega t_{s}},$

with
$\det[S^{\prime\prime}]=\partial_{tt}^{2}S\,\partial_{t^{\prime}t^{\prime}}^{2}S-(\partial_{tt^{\prime}}^{2}S)^{2}$
the Hessian of the action, thus giving the harmonic yield as a coherent sum of
discrete contributions coming from each of the quantum orbits. Fig. 2c shows a
typical example of how the individual contributions (blue and green curves)
get combined into a total harmonic yield, with quantum-path interference
beatings in the plateau as the two contributions go in and out of phase
Lewenstein1995 ; Zair2008 .

Within the SPA expression (7), the key component is the summation: this runs
over all the saddle points of the action which are compatible with a suitable
steepest-descent integration contour, a property that can be nontrivial to
determine. Fig. 2b shows a typical configuration, with the short- and long-
trajectory saddle points approaching each other over the harmonic plateau,
close to the real axis, performing an ‘avoided crossing’ at the harmonic
cutoff, and then advancing into imaginary time.

At this avoided crossing, the saddle-point pair experiences two important
transitions, known as the Stokes and anti-Stokes lines Figueira2002 ;
Milosevic2002 ; Chipperfield2007 . At the anti-Stokes transition, which is
defined as the point where $\imaginary(S)$ for the two orbits is equal, the
short-trajectory contribution begins to grow exponentially, and it must be
discarded from the summation in order to keep reasonably physical results.

This elimination is enforced by the Stokes transition, which typically happens
earlier, defined as the point where $\real(S)$ for the two orbits is equal and
then changes order. This is an important change, as the steepest-descent
method requires integration contours to follow the contour lines of
$\real(S)$: thus, the change in the ordering of
$\real(S(t_{s},t_{s}^{\prime}))$ for the long- and short-trajectory saddle
points means that, after the Stokes transition, the short-trajectory saddle
point (in this example) can no longer form part of a suitable integration
contour, and it needs to be discarded from the summation.

This means, however, that the SPA harmonic yield at the cutoff is a
discontinuous function of $\Omega$, coming from the discrete jump at the
points where one of the trajectories is eliminated, and this discontinuity is
clearly incompatible with the initial, obviously continuous, expressions for
$\mathbf{D}(\Omega)$ in (1) and (5). This apparent paradox is resolved by
noting that the SPA is not valid when saddles are close together, as the
quadratic approximation to the exponent fails. Instead, one must use a cubic
form for the exponent, which can be integrated in terms of Airy or fractional-
order Bessel functions. This gives a continuous spectrum from the saddle-point
pair, known as the Uniform Approximation Figueira2002 ; Milosevic2002  (UA),
which is shown in Fig. 2c.

From a more general viewpoint, it is important to stress that applying these
considerations is only possible once the various saddle points have been
classified into continuously-connected quantum orbits, as without that step it
is impossible to even define the objects that will be included in the
summation or discarded from it.444Similarly, keeping good track of the saddle
points in continuous quantum orbits is also essential to ensure that the
Hessian square root in (7),
$\sqrt{-\det[S^{\prime\prime}(t_{s},t^{\prime}_{s})]}$, does not cross any
branch cuts.  In practice, this is a high-friction point in the calculation,
requiring expensively tight grid spacings in $\Omega$ to accurately resolve
the avoided crossings, or the additional design of an adaptive grid to
increase the energy resolution there Chipperfield2007 ; Hoffmann2011 ; Das2017
. Moreover, this problem is compounded in more complex polychromatic fields,
where the saddle-point structure changes depending on the details of the laser
pulse. It is the goal of this manuscript to provide a simple and effective
method to classify these quantum orbits, separating the points in Fig. 2b
along the diagonal into the two curves shown.

###  I.2 A model for the saddle points at the cutoff

In order to build this method, we first consider a simple model for the saddle
points at and near the harmonic cutoff, in order to isolate the key features
of the problem. In essence, Fig. 2b shows us two saddle points approaching
each other and then receding, so we consider the simplest model action with
only two saddles, of the polynomial form

$\frac{\textrm{d}S}{\textrm{d}t}=A(t-t_{s,1})(t-t_{s,2}).$ (8)

(For simplicity, we restrict our attention for now to actions with only one
variable, $t$, which corresponds to solving (6b) first for
$t_{s}^{\prime}=t_{s}^{\prime}(t)$, giving $S(t)=S(t,t_{s}^{\prime}(t))$ and
then examining (6a). We shall lift this restriction later.) Moreover, we use
the clear symmetry evident in Fig. 2b, with both saddles (approximately)
symmetrically placed about the center of the plot, and enforce the condition
$t_{s,1}=-t_{s,2}$ in our model, so that its action obeys

$\displaystyle\frac{\textrm{d}S}{\textrm{d}t}$
$\displaystyle=A(t-t_{s})(t+t_{s})$ $\displaystyle=A(t^{2}-t_{s}^{2}),$ (9)

which can be integrated to find

$S(t)=\frac{1}{3}A\,t^{3}-At_{s}^{2}\,t+C$ (10)

as a cubic polynomial, where we set $C=0$ for simplicity.

Turning now to the oscillatory integral for this action, we can define it in
the form

$\displaystyle\int e^{-iS(t)}\textrm{d}t=\int
e^{-\frac{i}{3}At^{3}+iAt_{s}^{2}\,t}\,\textrm{d}t$ (11)

and then compare the linear term, $e^{+iAt_{s}^{2}\,t}$, with the Fourier
kernel $e^{+i\Omega t}$ from (1), which has the same structure. Thus, in order
to turn (11) into a form that is clearly analogous to (1), we separate
$At_{s}^{2}=\Omega-(\Omega_{c}+i\eta)$ into a variable and a constant part,
and thus define

$\displaystyle D(\Omega)=\int
e^{-\frac{i}{3}At^{3}+i(\Omega-\Omega_{c}-i\eta)t}\,\textrm{d}t.$ (12)

Here the constant part $\Omega_{\mathrm{hc}}=\Omega_{c}+i\eta$ has a nonzero
imaginary part $+i\eta$, as we do not have any guarantees that $At_{s}^{2}$ is
real, but its real part $\Omega_{c}$ can be set to zero if desired, since it
simply acts as an offset for the variable $\Omega$. Here the functional
dependence

$At_{s}^{2}=\Omega-(\Omega_{c}+i\eta)=\Omega-\Omega_{\mathrm{hc}}$ (13)

on $\Omega$ is an additional postulate, justified only by analogy with the
Fourier kernel of the full integral which our model attempts to mimic, but, as
we shall see shortly, the ‘quantum orbits’ $t_{s}(\Omega)$ that result from
this identification form a good model for the avoided-crossing behaviour in
Fig. 2b.

In addition, if we now separate $S(t)=S_{V}(t)-\Omega t$ as we did in (2)
above, the saddle-point equation (9) can now be rephrased as the requirement
that

$\displaystyle\frac{\textrm{d}S_{V}}{\textrm{d}t}=At^{2}+(\Omega_{c}+i\eta)=\Omega,$
(14a) with $\Omega$ running over the real axis, in direct analogy to (3a) and
(6a). This is our first key insight: the curves traced out by the solutions of
the saddle-point equation (9) can also be described as the contour line
$\imaginary\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}\right]\mathclose{}=0$
(14b)

of the derivative of the model Volkov action $S_{V}$ (and also of the full
action $S$, since they differ by a real number).

---
Figure 3:  Contour map of
$\imaginary\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}\right]\mathclose{}$
for the model action in (10), with blue (yellow) indicating negative
(positive) values, and with the
$\imaginary\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}\right]\mathclose{}=0$
contour highlighted in black. The gray dot-dashed lines are the contours of
$\real\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}\right]\mathclose{}$
passing through the central saddle point. We set $\eta=-1$ and
$A=e^{$10\text{\,}\mathrm{\SIUnitSymbolDegree}$i}$.

This insight then lights the way further: our principal task is to look for
objects between the two curves in Fig. 2b that will help us separate them, and
the contour map of
$\imaginary\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}\right]\mathclose{}$
is clearly the place to look. We show this contour map in Fig. 3, with the
zero contour of (14b) highlighted as a thick black curve, clearly showing the
avoided-crossing structure of Fig. 2b that we want to model. More importantly,
this plot shows us our second key insight: there is indeed a nontrivial object
separating the two quantum orbits, in the form of a saddle point in this
contour map. This is the central object we are after: the harmonic-cutoff time
$t_{\mathrm{hc}}$ for this model. Like all saddle points, this can be found as
a zero of the derivative of the contour map’s original function, or in other
words, via

$\frac{\textrm{d}^{2}S_{V}}{\textrm{d}t^{2}}(t_{\mathrm{hc}})=0.$ (15)

In the particular case of our model action (14a), for which
$\frac{\textrm{d}^{2}S_{V}}{\textrm{d}t^{2}}=2A\,t$, this saddle point lies at
the origin, due to the explicit choice of center of symmetry made by setting
$t_{s,1}=-t_{s,2}$ above. In the general case, (15) will find the centerpoint
between
$\imaginary\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}\right]\mathclose{}$
contours whenever they approach each other.

The appearance of a saddle point at the midpoint between contour lines in
close proximity is a generic feature of contour landscapes. Similar structures
have been explored previously Pisanty2016slalom ; Keil2016 , in the context of
tunnel ionization.

###  I.3 The landscape at the harmonic-cutoff point

Having found the saddle point, we can now use it directly, since we can fully
reconstruct the model action $S_{V}(t)$ using only its behaviour at the
center, by taking derivatives:

$\displaystyle\frac{\textrm{d}S_{V}}{\textrm{d}t}(t_{\mathrm{hc}})$
$\displaystyle=\Omega_{\mathrm{hc}}=\Omega_{c}+i\eta,$ (16a)
$\displaystyle\frac{\textrm{d}^{3}S_{V}}{\textrm{d}t^{3}}(t_{\mathrm{hc}})$
$\displaystyle=2A,$ (16b)

with the second derivative vanishing by construction.

---
Figure 4:  Contour map of $\imaginary\mathopen{}\left[\frac{\partial
S_{V}}{\partial t}(t,t_{s}^{\prime}(t))\right]$ over the complex recollision-
time plane for the full Volkov action from (2b). The saddle-point trajectories
of Fig. 1 appear clearly as the zero contour (thick lines), i.e. as the
complex times where $\frac{\partial S_{V}}{\partial t}(t,t_{s}^{\prime}(t))$
is real-valued. The harmonic-cutoff points of Fig. 1c (triangles) are,
correspondingly, the saddle points of this derivative; the coloured lines show
the contours at these saddle points.

We begin by focusing on the first derivative, which gives the linear term in
the action, whose coefficient is directly connected to the harmonic frequency
via $At_{s}^{2}=\Omega-(\Omega_{c}+i\eta)$ as set above. In particular, this
term directly encodes the saddle-point solutions, and we can recover these
explicitly by inverting that relationship:

$t_{s}=\pm\sqrt{\frac{\Omega-(\Omega_{c}+i\eta)}{A}}.$ (17)

Here the square in $At_{s}^{2}$ implies that $t_{s}$ appears as an explicit
square root, with the sign ambiguity giving us the two saddle-point solutions
of the problem.

This is also a key structural insight, which has largely remained unobserved
in the literature: the various solutions of saddle-point equations are simply
different branches of a unified Riemann surface, examined at the real axis of
the target space. This observation is made explicit in the branch choice given
by the $\pm$ sign in (17), but the same is true in the full problem: if we
approach the coupled equations in (6) by solving (6b) first for
$t_{s}^{\prime}=t_{s}^{\prime}(t)$, then (6a) can be expressed in the form

$\frac{\partial S_{V}}{\partial
t}(t,t_{s}^{\prime}(t))=\frac{1}{2}\left(\mathbf{p}_{s}(t,t_{s}^{\prime}(t))+\mathbf{A}(t)\right)^{2}+I_{p}=\Omega,$
(18)

where it explicitly involves the inverse of the many-to-one, complex-valued
function $\frac{\partial S_{V}}{\partial t}(t,t_{s}^{\prime}(t))$, and this is
precisely the type of problem encoded by Riemann surfaces. (We explore this
perspective in detail in Section I.5 below.)

Moreover, as far as the full HHG problem is concerned, our insight in (14b)
above that the saddle-point trajectories are simply contour lines of a
derivative of the action remains unchanged, with the obvious alterations to

$\imaginary\mathopen{}\left[\frac{\partial S_{V}}{\partial
t}(t,t_{s}^{\prime}(t))\right]\mathclose{}=0.$ (19)

We exemplify this directly in Fig. 4, showing the contour map of
$\imaginary\mathopen{}\left[\frac{\partial S_{V}}{\partial
t}(t,t_{s}^{\prime}(t))\right]\mathclose{}$ for the full HHG problem as in
Fig. 1, and there the zero contour lines (highlighted in black) are indeed the
curves followed by the typical quantum orbits shown in Fig. 1c. Similarly, the
harmonic-cutoff times $t_{\mathrm{hc}}$ used there are found as the saddle
points of the $\frac{\partial S_{V}}{\partial t}(t,t_{s}^{\prime}(t)){}$
landscape in Fig. 4.

That said, the quantum-orbit curves of Fig. 4 are not simply unstructured
curves, as they have an explicit parametrization in terms of the harmonic
frequency $\Omega$, but the same is true for the contours of
$\imaginary\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}\right]\mathclose{}$,
which can also be seen as parametrized by the real part of the function. This
brings physical content to the first part of the derivatives we found in (16),
where the real part of $\frac{\textrm{d}S_{V}}{\textrm{d}t}(t_{\mathrm{hc}})$
is simply the frequency offset $\Omega$. This offset is carried from the
saddle point $t_{\mathrm{hc}}$ to the quantum-orbit lines by means of the
contour lines of
$\real\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}\right]\mathclose{}$.
These are shown dot-dashed in Fig. 4, and they clearly intersect the quantum-
orbit tracks at the point where they are closest to each other, and thus at
the transition itself. This is then what allows us to identify

$\Omega_{c}=\real\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}(t_{\mathrm{hc}})\right]\mathclose{}$
(20a) as the frequency of the harmonic cutoff.

Having made this leap, we must now confront the fact that
$\frac{\textrm{d}S_{V}}{\textrm{d}t}(t_{\mathrm{hc}})$ as given in (16a) also
has a nonzero imaginary part,

$\eta=\imaginary\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}(t_{\mathrm{hc}})\right]\mathclose{},$
(20b)

which similarly needs to be addressed. This imaginary part has already played
a role, and it is explicitly shown in Fig. 3 as the height of the
$\frac{\textrm{d}S_{V}}{\textrm{d}t}$ saddle with respect to the zero contour
where the quantum orbits lie. As such, the imaginary part $\eta$ of the
harmonic cutoff, defined by (16a), controls how closely the two quantum orbits
approach each other at the cutoff, and thus how distinguishable they are. As
we shall see in Section II.1 below, this distinguishability further controls
the strength of quantum-path interference between them.

Moreover, the sign of $\eta$ controls the direction in which the quantum
orbits turn when they approach each other, and thus which of the two post-
cutoff evanescent solutions, at positive and negative imaginary time,
corresponds to the ‘short’ and ‘long’ trajectory to the left and right of
$t_{\mathrm{hc}}$. Generally, this imaginary part of
$\frac{\textrm{d}S_{V}}{\textrm{d}t}$ is fixed by the problem and it can be
either positive or negative, but, as we shall see in Section IV, sign changes
in $\eta$ are generic behaviour when the driving pulse shape changes depending
on one or more parameters, which implies reconnections and changes in topology
for the curves traced out by the quantum orbits. This further emphasizes the
fact that the different saddle points are essentially one and the same object,
corresponding only to different branches of one unified Riemann surface.

However, this point can be observed more cleanly, by simply allowing the
harmonic frequency $\Omega$ to acquire complex values, i.e., by considering
the analytical continuation of $D(\Omega)$. Doing this amounts to directly
affecting the value of $\eta$ in (14a), and if $\Omega$ is taken at the
complex harmonic cutoff itself, $\Omega_{c}+i\eta$, then the constant term in
the saddle-point equation vanishes, leaving a double zero of the form

$At^{2}=0$ (21)

at $t=t_{\mathrm{hc}}$. In other words, the harmonic-cutoff times
$t_{\mathrm{hc}}$ are the locations where the saddle points for the two
quantum orbits fully coalesce into a single point.

###  I.4 The orientation of the harmonic-cutoff saddle and quantum-orbit
classification

Having examined the first derivative of $S_{V}$ at the harmonic-cutoff time,
we now turn to the role of the third derivative, given in (16b). This term
gives the coefficient of the quadratic term in
$\frac{\textrm{d}S_{V}}{\textrm{d}t}$ and, as such, it controls the
orientation of the saddle in $\frac{\textrm{d}S_{V}}{\textrm{d}t}$ at
$t_{\mathrm{hc}}$, as shown in Fig. 3. (Indeed, we chose a nonzero phase for
$A=e^{$10\text{\,}\mathrm{\SIUnitSymbolDegree}$i}$ for that configuration to
emphasize that this saddle need not be neatly oriented, as observed in Fig.
4.) Retrieving this orientation is essential in order to use our harmonic-
cutoff times $t_{\mathrm{hc}}$ to classify the saddle points into quantum
orbits: in order to turn the clouds of points in Fig. 1b into the ordered
curves of Fig. 1c, we also require the direction of the separatrix that goes
through the missed approach.

To obtain this direction, we look for the vector that goes from
$t_{\mathrm{hc}}$ to the closest point on the
$\imaginary\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}\right]\mathclose{}=0$
contour that the quantum orbits follow. As we noted above, this point occurs
at $\Omega=\Omega_{c}$, which means that we simply need to solve the saddle-
point equation (14a) at that frequency, i.e.,

$\displaystyle At_{\mathrm{sep}}^{2}+i\eta=0\ \implies\
t_{\mathrm{sep}}=\sqrt{-i\eta/A},$ (22)

or, in terms of the explicit derivatives of the action and with an explicit
relation to the saddle center $t_{\mathrm{hc}}$,

$\displaystyle\delta
t_{\mathrm{sep}}=t_{\mathrm{sep}}-t_{\mathrm{hc}}=\sqrt{-2i\frac{\imaginary\mathopen{}\left[\displaystyle\frac{\textrm{d}S_{V}}{\textrm{d}t}(t_{\mathrm{hc}})\right]\mathclose{}}{\displaystyle\frac{\textrm{d}^{3}S_{V}}{\textrm{d}t^{3}}(t_{\mathrm{hc}})}}.$
(23)

We can then write the explicit condition for the separatrix by treating
$\delta t_{\mathrm{sep}}$ and $t_{s}-t_{\mathrm{hc}}$ as vectors and asking
for the sign their inner product. Thus, the criterion that separates the two
quantum orbits on either side of a given harmonic-cutoff time
$t_{\mathrm{hc}}$ is the sign comparison

$\real\mathopen{}\Big{[}(t_{s}-t_{\mathrm{hc}})^{*}\,\delta
t_{\mathrm{sep}}\Big{]}\mathclose{}\lessgtr 0,$ (24)

with $t_{\mathrm{hc}}$ defined as in (15) and $\delta t_{\mathrm{sep}}$
defined as in (23).

In a practical calculation, there will typically be multiple harmonic-cutoff
times, alternating between near-threshold harmonic frequencies, where the
quantum orbits first appear, and high-frequency cutoffs. To complete our
classification scheme – shown as the background coloured zones in Fig. 1c,
which we repeat in some additional detail in Fig. 5 – we separate the complex
time plane into strips using the mid-points between the real parts of the
successive harmonic-cutoff times, and then divide each strip in two using the
criterion in (24).

Here it is important to note that the definition in (23) contains a sign
ambiguity coming from the choice of sign for the square root. In practice, the
principal branch of the square root tends to work well for most cases, but the
branch choice there does occasionally require dedicated attention,
particularly for the first low-energy harmonic cutoff at the start of the
series.

The results of our classification procedure are shown in Fig. 5: the
separatrices obtained from the sign-comparison criterion (24) break up the
complex plane into trapezoids that contain one and only one quantum orbit
each, and they can thus be used to classify the saddle-point solutions into
well-defined families in a uniform and robust fashion.

---
Figure 5:  Results of our saddle-point classification procedure, as shown
previously in Fig. 1c. We find harmonic-cutoff times $t_{\mathrm{hc}}$
(triangles) as solutions of (15), with the separatrix through each
$t_{\mathrm{hc}}$ given by the sign comparison in (24) with $\delta
t_{\mathrm{sep}}$ defined in (23); we apply this criterion on vertical strips
obtained by taking the midpoints between the real parts of the
$t_{\mathrm{hc}}$. This then defines clear regions occupied by each of the
quantum orbits, which can be unambiguously labelled.
---
List of tdfigures 1  Topology of the Riemann surface of the quantum orbits.
This is the manifold
$\mathcal{S}=\\{(t,\Omega)\in\mathbb{C}^{2}:\Omega=\frac{\textrm{d}S_{V}}{\textrm{d}t}(t)\\}$,
which has complex dimension 1 and is embedded in the space $\mathbb{C}^{2}$,
with real dimension 4, so we can only plot it by projecting out one component:
we show its projections to $\real(\omega t)$ in (a, b) and to
$\imaginary(\omega t)$ in (c-f). The projection to $\real(\omega t)$ shows the
topology as a single sheet which wraps around itself, covering the
$\Omega\in\mathbb{C}$ plane multiple times when the $\real(\omega t)$
coordinate is projected out, so that $\frac{\textrm{d}S_{V}}{\textrm{d}t}$ has
a multi-valued inverse. To separate this multi-valued inverse into valid
single-valued functions, the sheet needs to be cut into separate branches:
this is the role of the separation into coloured regions from the
classification scheme in Fig. 5, which carry over to the multiple coloured
sections in (a). The sheets connect to each other via square-root-type branch
cuts, and we show a detail of the first such connection in (b). The projection
onto $\imaginary(\omega t)$ forms a self-intersecting surface which loops
around itself multiple times, so we plot each of the branches separately (with
half of each of its neighbours shown half-transparently) in (c-f). Interactive
versions and 3D-printable models of these plots are available at imaginary-
harmonic-cutoff.github.io SupplementaryMaterial .

###  I.5 The quantum-orbit Riemann surface

We now return to an observation we made earlier: the saddle-point equations
for HHG represent the inverse of a many-to-one analytical function,
$\frac{\partial S_{V}}{\partial t}$ as given in (18), evaluated on the real
axis in $\Omega$, and this is precisely the definition of a Riemann surface.
As such, it is pertinent to study this Riemann surface as a whole, since its
topology and geometry are the key factors that govern the quantum-orbit
classification.

The Riemann surface here is the set

$\mathcal{S}=\left\\{(t,\Omega)\in\mathbb{C}^{2}:\Omega=\frac{\textrm{d}S_{V}}{\textrm{d}t}(t)\right\\},$
(25)

and it forms a manifold of complex dimension 1 (i.e., a manifold where each
point has a neighbourhood homeomorphic with the complex unit disk) embedded in
a space of complex dimension 2 (and therefore of real dimension 4). This space
is too big to be visualized directly, so we approach it by projecting down to
the real and imaginary parts of $t$, i.e., by projecting the surface to the
spaces $(\real(\Omega),\imaginary(\Omega),\real(\omega t))$ and
$(\real(\Omega),\imaginary(\Omega),\imaginary(\omega t))$, which we show in 3D
Fig. 1.

The topology is most clearly displayed in the $\real(\omega t)$ projection,
shown in 3D Fig. 11: the Riemann surface consists of a series of connected
sheets which connect sequentially to each other, one by one, as $\real(\omega
t)$ increases. In essence, the topology and the overall geometry here are
those of the Riemann surface of $\sin(z)$, with the sheets connected pairwise
by square-root branch points (locally homeomorphic to the quadratic model
discussed in Section I.2) as shown in detail in 3D Fig. 11.

This Riemann surface carries an image of the $\omega t$ complex plane (so, in
particular, the mesh on the surface represents a square grid on that plane),
but it forms a multiple cover of the image plane $\Omega\in\mathbb{C}$. To get
the inverse, then, the surface must be split into separate branches, and this
is precisely the role of the coloured regions shown in Fig. 5 as produced by
our saddle-point classification algorithm: this colouring is retained in the
Riemann surface as displayed in 3D Fig. 1, and each individual region, when
projected down to the $\Omega$ plane, forms a single full cover of the complex
plane.555This single-cover property is approximate, as there is still a small
amount of double-cover overlap in regions next to the separatrix when away
from the branch points at the $t_{\mathrm{hc}}$. This can be fixed if
necessary, but it is not central to our argument here.  In other words, each
of the coloured regions in Fig. 5 forms the image of a single-valued branch of
the inverse of $\frac{\textrm{d}S_{V}}{\textrm{d}t}(t)$.

The imaginary-part projection onto $\imaginary(\omega t)$ is slightly harder
to represent and visualize, as the corresponding surface folds back on itself,
with multiple self-intersections. However, the separation of this surface into
individual branches also solves this problem: when plotting one branch at a
time, as we show in 3D Fig. 11-1, the self-intersections disappear – as they
must if the surface represents a single-valued function – and they are only
present as intersections between different branches.

Physically, the most important part of this Riemann surface is its
intersection with the $\Omega\in\mathbb{R}^{+}$ positive real axis in the
$\Omega$ plane, which, as discussed in Section I.3 above, forms the quantum
orbits themselves. These are shown highlighted as curves in 3D Fig. 1, and
their role here is most clearly apparent in 3D Fig. 11: this intersection
forms the $(\real(\omega t),\real(\Omega))$ energy-time mapping for the
quantum orbits, a well-known (and deeply physical) part of the theory
Ivanov2014 ; Nayak2019 . The quantum-orbit Riemann surface, together with its
topology, is nothing more than the analytical continuation of this standard
energy-time mapping.

###  I.6 Derivatives of the action in the two-variable case

Before moving on, it is important to define in more detail the relationship
between our simple model and the full-blown HHG integral, whose action has two
variables instead of one. This is possible, as mentioned above, by using (6b)
to define $t_{s}^{\prime}=t_{s}^{\prime}(t)$ and with that the single-variable
$S_{V}(t)=S_{V}(t,t_{s}^{\prime}(t))$, but that only works explicitly for the
values of the action, and its derivatives need to be considered more
carefully.

The first derivative is not affected, since

$\displaystyle\frac{\textrm{d}S_{V}}{\textrm{d}t}(t)$
$\displaystyle=\frac{\textrm{d}}{\textrm{d}t}S_{V}(t,t_{s}^{\prime}(t))$
$\displaystyle=\frac{\partial S_{V}}{\partial
t}(t,t_{s}^{\prime}(t))+\frac{\textrm{d}t_{s}^{\prime}}{\textrm{d}t}(t)\frac{\partial
S_{V}}{\partial t^{\prime}}(t,t_{s}^{\prime}(t)),$ (26)

and here the partial derivative in the second term, $\frac{\partial
S_{V}}{\partial t^{\prime}}(t,t_{s}^{\prime}(t))$, vanishes by the definition
of $t_{s}^{\prime}(t)$. However, if we now turn to the second derivative, the
procedure no longer works, and the equivalent calculation,

$\displaystyle\frac{\textrm{d}^{2}S_{V}}{\textrm{d}t^{2}}(t)$
$\displaystyle=\frac{\textrm{d}}{\textrm{d}t}\frac{\partial S_{V}}{\partial
t^{\prime}}(t,t_{s}^{\prime}(t))$
$\displaystyle=\frac{\partial^{2}S_{V}}{\partial
t^{2}}(t,t_{s}^{\prime}(t))+\frac{\textrm{d}t_{s}^{\prime}}{\textrm{d}t}(t)\frac{\partial^{2}S_{V}}{\partial
t\,\partial t^{\prime}}(t,t_{s}^{\prime}(t)),$ (27)

returns a term that includes
$\frac{\textrm{d}t_{s}^{\prime}}{\textrm{d}t}(t)$. This cannot be evaluated
explicitly within this system, as $t_{s}^{\prime}(t)$ itself is only defined
implicitly.

To resolve this, we use the implicit definition of $t_{s}^{\prime}(t)$,

$\frac{\partial S_{V}}{\partial t^{\prime}}(t,t_{s}^{\prime}(t))\equiv 0,$
(28)

and then differentiate with respect to $t$, to obtain

$\displaystyle\frac{\partial^{2}S_{V}}{\partial t\,\partial
t^{\prime}}(t,t_{s}^{\prime}(t))+\frac{\textrm{d}t_{s}^{\prime}}{\textrm{d}t}(t)\frac{\partial^{2}S_{V}}{\partial
t^{\prime 2}}(t,t_{s}^{\prime}(t))\equiv 0,$ (29)

which can then be solved for the implicit derivative as

$\displaystyle\frac{\textrm{d}t_{s}^{\prime}}{\textrm{d}t}(t)=-\frac{\partial^{2}S_{V}}{\partial
t\,\partial
t^{\prime}}(t,t_{s}^{\prime}(t))\bigg{/}\frac{\partial^{2}S_{V}}{\partial
t^{\prime 2}}(t,t_{s}^{\prime}(t)).$ (30)

Finally, this can be substituted into (27) to give the form

$\displaystyle\frac{\textrm{d}^{2}S_{V}}{\textrm{d}t^{2}}(t)$
$\displaystyle=\frac{\displaystyle\frac{\partial^{2}S_{V}}{\partial
t^{2}}\frac{\partial^{2}S_{V}}{\partial t^{\prime
2}}-\left(\frac{\partial^{2}S_{V}}{\partial t\,\partial
t^{\prime}}\right)^{2}}{\displaystyle\frac{\partial^{2}S_{V}}{\partial
t^{\prime 2}}}.$ (31)

Here we have dropped the right-hand side evaluations at
$(t,t_{s}^{\prime}(t))$ for simplicity, but it is also possible to understand
the right-hand side as a function of $(t,t^{\prime})$ as independent
variables, whose zeros are to be found in conjunction with those of (6b), and
this is an easier approach in practice.

Finally, to obtain the higher-order derivatives (in particular, the third
derivative, which gives the cubic coefficient $A$) we use symbolic computation
RBSFA ; FigureMaker , using the substitution rule (30) at each step to retain
explicit formulations.

###  I.7 The threshold harmonics for linear fields

For a general driving field, the harmonic cutoffs given by our definition will
appear at both ends of the harmonic plateau, oscillating between low
frequencies near the ionization threshold and high energies at the classical
cutoff. In general, moreover, all of these harmonic cutoffs will have nonzero
imaginary parts, corresponding to missed approaches at a finite distance
between the quantum orbits involved. However, as can be seen in Fig. 4 and
Fig. 1c, this is not the case for the linearly-polarized field we use above,
since half of the harmonic-cutoff saddles $t_{\mathrm{hc}}$ lie directly on
the quantum-orbit line, i.e., the quantum orbits have a full crossing there.

This behaviour is generic to all linearly-polarized fields, whether
monochromatic or not, and it implies that the low-energy cutoffs always lie at
exactly $\Omega_{c}+i\eta=I_{p}$, and the threshold harmonics at
$\Omega=I_{p}$ always involve an exact saddle coalescence. In other words, for
threshold harmonics, the saddles of the action landscape are not gaussian:
instead, they are exact monkey saddles, as shown in Fig. 7a below; this was
noticed as early as Ref. Lewenstein1994, , but has not drawn much attention
since then.666This is largely because the approximations that produced the SFA
integral are known to be inaccurate for the lower plateau and threshold
harmonics.

To understand this behaviour, we return to the saddle-point equations (6),
which for linearly-polarized fields can be rephrased in the simpler, scalar
form

$\displaystyle\frac{\partial S_{V}}{\partial
t}=\frac{1}{2}\left(p_{s}(t,t^{\prime})+A(t)\right)^{2}+I_{p}$
$\displaystyle=\Omega,$ (32a) $\displaystyle-\frac{\partial S_{V}}{\partial
t^{\prime}}=\frac{1}{2}\left(p_{s}(t,t^{\prime})+A(t^{\prime})\right)^{2}+I_{p}$
$\displaystyle=0.$ (32b)

Here the first equation is the crucial one, since at $\Omega=I_{p}$ it tells
us that

$\frac{1}{2}v_{r}(t,t^{\prime})^{2}=\frac{1}{2}\left(p_{s}(t,t^{\prime})+A(t)\right)^{2}=0,$
(33)

i.e., the recollision velocity $v_{r}(t,t^{\prime})=p_{s}(t,t^{\prime})+A(t)$
vanishes exactly, with a double zero which still vanishes after taking one
derivative.

This becomes particularly important when we consider the second derivative of
the action, (31), whose zeros determine the harmonic-cutoff times. The partial
derivatives involved in (31) can be evaluated by differentiating (32a) with
respect to $t$ and $t^{\prime}$, which yields

$\displaystyle\frac{\partial^{2}S_{V}}{\partial t^{2}}$
$\displaystyle=\left(p_{s}+A(t)\right)\left[\frac{\partial p_{s}}{\partial
t}+\dot{A}(t)\right],$ (34a) $\displaystyle\frac{\partial S_{V}}{\partial
t\,\partial t^{\prime}}$ $\displaystyle=\left(p_{s}+A(t)\right)\frac{\partial
p_{s}}{\partial t^{\prime}}.$ (34b)

Both of these derivatives share a common factor of
$v_{r}(t,t^{\prime})=p_{s}(t,t^{\prime})+A(t)$, and this then means that the
second derivative $\frac{\textrm{d}^{2}S_{V}}{\textrm{d}t^{2}}(t)$ will always
vanish when evaluated at solutions of the first-order saddle-point equations
(32) at $\Omega=I_{p}$, completing the proof.

It is important to note that, when a saddle-point coalescence like this
occurs, the classification scheme embodied by our test in (24) breaks down
since, at the double zero,
$\frac{\textrm{d}S_{V}}{\textrm{d}t}(t_{\mathrm{hc}})$ also vanishes, so the
direction vector $\delta t_{\mathrm{sep}}$ from (23) also vanishes. In
practice, this can cause numerical instability, with
$\imaginary\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}(t_{\mathrm{hc}})\right]\mathclose{}$
evaluating to machine precision with a fluctuating sign, and this needs to be
handled explicitly – which could be as simple as assigning an artificial
nonzero
$\imaginary\mathopen{}\left[\frac{\textrm{d}S_{V}}{\textrm{d}t}(t_{\mathrm{hc}})\right]\mathclose{}$
to be used there, as either direction will work well for the separatrix. (For
the plots in Fig. 4 and Fig. 1c, we used a small ellipticity of
$\epsilon=0.01\%$ to break the degeneracy and avoid this instability.)

At an exact coalescence, it is impossible to classify the saddles in a unique
way: two saddles come in and two saddles come out, but there is no unique way
to match either of the outgoing saddles with either of the incoming ones. This
happens for all linearly-polarized fields at the $\Omega=I_{p}$ threshold, but
it can also happen at the harmonic cutoff in isolated cases when the field
depends on a variable parameter, as we will exemplify in Section IV below.

##  II The Harmonic-Cutoff Approximation

Having examined the structure of the action at the harmonic-cutoff times
$t_{\mathrm{hc}}$, in this section we will explore the behaviour of the
oscillatory integral around it, which we will employ to build the Harmonic-
Cutoff Approximation (HCA), an efficient method for estimating the harmonic
dipole around the cutoff using only information at $t_{\mathrm{hc}}$ itself.

###  II.1 Airy-function representation for the model case

To do this, we return to the model integral from (12), where we set as the
object of interest the integral the exponential of our model action (10),

$\displaystyle
D(\Omega)=\int_{C}e^{-\frac{i}{3}At^{3}+i(\Omega-\Omega_{\mathrm{hc}})t}\,\textrm{d}t,$
(35)

over an integration contour $C$ which should start at the ionization time
$t^{\prime}$ (or some suitably compatible valley of $\imaginary(S)$ to the
left of $t_{\mathrm{hc}}$) and end at $t\to\infty$. This integral is
essentially in Airy form, and it is almost structurally identical to the Airy
function’s integral representation (NIST_handbook, , Eq. (9.5.E4)),

$\operatorname{Ai}(z)=\frac{1}{2\pi i}\int_{\infty e^{-i\pi/3}}^{\infty
e^{i\pi/3}}\exp(\tfrac{1}{3}t^{3}-zt)\textrm{d}t.$ (36)

This is expected, since the harmonic cutoff is a ‘fold’ catastrophe, whose
associated diffraction-catastrophe integral is the Airy function Berry1980 .

To bring our representation (35) into the canonical form (36), the core
transformation is to eliminate the coefficient in front of the cubic term,

$-\frac{i}{3}At^{3}\longmapsto\frac{1}{3}\tilde{t}\,^{3}.$ (37)

In a sense, this is relatively simple, as it boils down to a change in
integration variable to $\tilde{t}=iA^{1/3}\,t,$ but there is an added
complication in that the radical $A^{1/3}$ of the cubic coefficient admits
three separate branches,

$\tilde{t}=ie^{2\pi ik/3}A^{1/3}\,t,$ (38)

for $k\in\\{0,1,2\\}$, which requires dedicated attention.

The variable change itself is essentially trivial, and (35) transforms under
(38) to

$\displaystyle D(\Omega)$ $\displaystyle=\frac{e^{-2\pi
ik/3}}{iA^{1/3}}\\!\\!\int_{\tilde{C}}\\!\exp(\frac{1}{3}\tilde{t}\,^{3}-\frac{\Omega_{\mathrm{hc}}-\Omega}{e^{2\pi
ik/3}A^{1/3}}\tilde{t})\,\textrm{d}\tilde{t},$ (39a) with
$\displaystyle\tilde{C}=ie^{2\pi ik/3}A^{1/3}\,C.$ (39b)

This is essentially in explicit Airy form, so long as the contour is correct,
and thus we can substitute in the Airy function as

$\displaystyle D(\Omega)$ $\displaystyle=\frac{2\pi}{e^{2\pi
ik/3}A^{1/3}}\operatorname{Ai}\mathopen{}\left(\frac{\Omega_{\mathrm{hc}}-\Omega}{e^{2\pi
ik/3}A^{1/3}}\right)\mathclose{},$ (40)

where $k$ needs to be chosen such that the altered contour $\tilde{C}=ie^{2\pi
ik/3}A^{1/3}\,C$ is compatible with the standard contour in (36), which we
depict in Fig. 6.

|  |
---|---|---

(a)

(b)

(c)

Figure 6:  Contour map of the exponent $\real(\tfrac{1}{3}t^{3}-zt)$ of the
integral representation of the Airy function, (36), with the standard contour
shown as the black arrows, for (a) coalescent saddles, (b) the oscillatory
regime, and (c) the evanescent part. The details of the contour can be altered
as required (such as e.g. to pass through the two saddle points in (b) along
steepest-descent contours), under the constraints that the contour start at
infinity at $-\tfrac{\pi}{2}<\arg(t)<-\tfrac{\pi}{6}$ and end at
$\tfrac{\pi}{6}<\arg(t)<\tfrac{\pi}{2}$, i.e., that it go down the valleys
that surround $e^{-i\pi/3}$ and $e^{i\pi/3}$.

---

(a)

(b)

(c)

Figure 7:  Contour map of $\imaginary[S(t,t_{s}^{\prime}(t))-\Omega t]$,
presented as in Fig. 4, taking a complex $\Omega$ equal to $\Omega_{c}+i\eta$
at the various harmonic-cutoff times $t_{\mathrm{hc}}$, causing the nearby
saddles to fully coalesce, and marking a complete breakdown of the SPA. To
apply the Harmonic-Cutoff Approximation, each monkey saddle (marked by the
converging coloured contours) must be rotated by $\tilde{C}=ie^{2\pi
ik/3}A^{1/3}\,C$, choosing the branch index $k$ so that the integration
contour matches the canonical Airy contours of Fig. 6. The phase $\arg(A)$ is
inset in (b) and (c), and it equals
$-87.97\text{\,}\mathrm{\SIUnitSymbolDegree}$,
$-$0.011\text{\,}\mathrm{\SIUnitSymbolDegree}$$ and
$-$0.155\text{\,}\mathrm{\SIUnitSymbolDegree}$$, resp., for the three monkey
saddles shown in (a), consistently with their different orientations.

In principle, a rigorous approach to the contour-choice problem requires a
careful examination of the constrains on the contour of the original integral,
as shown in Fig. 7, to determine the correct integration contour $C$ that is
compatible with the original integration limits; this is then rotated by
$iA^{1/3}$ and any additional factors of $e^{2\pi ik/3}$ necessary for the
contour to be in canonical form. (Moreover, care must be taken when taking the
radical $A^{1/3}$ over a parameter scan, since, as mentioned in Fig. 7, $A$
can lie very close to commonly-used choices for the branch cut of the
radical.)

In practice, if there are known constraints on the behaviour of the integral
(in particular, exponential decay at $\Omega>\Omega_{c}$), there will
typically only be one choice of $k$ compatible with the constraints, which can
then be selected on physical grounds. As a general rule, exponential decay at
$\Omega>\Omega_{c}$ requires $e^{2\pi ik/3}A^{1/3}$ to have a negative real
part.

Our result for the model integral, (40), now allows us to have a closer look
at the role of the imaginary part of the complex harmonic-cutoff energy
$\Omega_{\mathrm{hc}}$, which we introduced in (20b). As we mentioned then,
this imaginary part controls the strength of quantum-path interference between
the two quantum orbits that meet at the relevant cutoff, and we show this in
Fig. 8. The HCA approximates the integral (35) as an explicit Airy function in
$\Omega$ (red line), which transitions from oscillatory to exponential-decay
behaviour at $\Omega_{c}=\real(\Omega_{\mathrm{hc}})$.

However, in order for the interference fringes in the oscillatory regime to
have a full contrast and pass through the zeros of the Airy function, the
argument must be real, and this is generally impossible if
$\eta=\imaginary(\Omega_{\mathrm{hc}})$ is nonzero. To increase the
interference features, thus, we can artificially reduce (blue line) or
eliminate (green line) this imaginary part; conversely, increasing $\eta$
(purple and magenta lines) further damps the interference, and it introduces
exponential behaviour into the plateau. That said, to get full contrast in the
interference it is also necessary for the denominator to be real-valued, which
we show as the gray line by artificially adjusting it.

The nonzero imaginary parts of the constants in the argument of the Airy
function also makes this approximation distinct from previous approaches that
regularize the cutoff in terms of single Airy functions Ivanov2014 ;
Frolov2009 ; Frolov2010 , making it better able to capture the QPI contrast in
the neighbourhood of the cutoff. On the other hand, our approximation sits one
level below the Uniform Approximation Figueira2002 ; Milosevic2002 ; Wong2001
(as well as the equivalent constructions in Refs. Frolov2012, ;
Sarantseva2013, ; Okajima2012, ), which expands the prefactor to subleading
order and is thus able to capture variations in the QPI contrast coming from
the prefactor’s effect on the saddles. That said, the HCA is, at least in
principle, strictly local to the cutoff, and this makes it especially suited
for efficient evaluation and clear analysis of the harmonic strength (as well
as phase properties Khokhlova2016 ) at the cutoff.

###  II.2 The full HHG integral

We now return to the full two-dimensional integral for HHG, (5), and transform
it into the model form (35) so that the HCA can be used to estimate it. For
simplicity, we write the integral in explicit prefactor-action form,

$\displaystyle\mathbf{D}(\Omega)$
$\displaystyle=\int_{-\infty}^{\infty}\\!\\!\\!\textrm{d}t\int_{-\infty}^{t}\\!\\!\\!\textrm{d}t^{\prime}\
\mathbf{f}(t,t^{\prime})\,e^{-iS_{V}\\!(t,t^{\prime})+i\Omega t}$ (41a)
$\displaystyle\mathbf{f}(t,t^{\prime})$
$\displaystyle=\frac{(2\pi/i)^{3/2}}{(t-t^{\prime})^{3/2}}\mathbf{d}(\mathbf{p}_{s}(t,t^{\prime}){+}\mathbf{A}(t))\Upsilon(\mathbf{p}_{s}(t,t^{\prime}){+}\mathbf{A}(t^{\prime})).$
(41b)

To transform this into the model form (35), we need to reduce the integral to
a single dimension, and then translate the origin to the relevant harmonic-
cutoff time $t_{\mathrm{hc}}$.777However, it is important to point out that
the process is rather more general than this. As was pointed out in the
construction of the Uniform Approximation Figueira2002  (as well as its
earlier analogues in a semiclassical context Schomerus1997 ), the coordinate
separation done here is in essence an application of the splitting lemma of
catastrophe theory Poston1978 , which allows the ‘fold’ catastrophe encoded by
the Airy function to be isolated into a single coordinate axis. As such, the
simple model of Section I.2 is not a ‘toy’ model in any sense: instead, it is
a universal model, which is fully capable of capturing the (local) behaviour
of the integral.

---
Figure 8:  Control of the imaginary part of the harmonic cutoff on quantum-
path interference. The red line shows the Airy-function dependence of the
Harmonic-Cutoff Approximation, as in (40), for $\Omega_{\mathrm{hc}}$ and $A$
corresponding to the first-return cutoff in Fig. 1. The coloured lines have
the imaginary part $\eta=\imaginary(\Omega_{\mathrm{hc}})$ artificially
reduced or amplified by a variable factor $r$, which respectively amplifies or
damps the QPI features; this corresponds to choosing a lower or higher
$\imaginary(\frac{\partial S_{V}}{\partial t})$ contour to follow in Fig. 3.
To obtain full interference, $r$ must be set to $0$, but the denominator must
also be real and negative, i.e., we replace $e^{2\pi ik/3}A^{1/3}$ with
$-|A|^{1/3}$. The black vertical line is at $\Omega=\Omega_{c}$.

The first is achieved, as mentioned earlier, by doing a saddle-point
approximation on the $t^{\prime}$ integral only, and this returns the harmonic
dipole in the form

$\displaystyle\mathbf{D}(\Omega)$
$\displaystyle=\sum_{s}\int_{-\infty}^{\infty}\textrm{d}t\,\sqrt{\frac{2\pi}{i\frac{\partial^{2}S_{V}}{\partial
t^{\prime 2}}(t,t_{s}^{\prime}(t))}}\mathbf{f}(t,t_{s}^{\prime}(t))$ (42)
$\displaystyle\qquad\qquad\qquad\qquad\times
e^{-iS_{V}\\!(t,t_{s}^{\prime}(t))+i\Omega t},$

with $t^{\prime}_{s}(t)$ defined implicitly via (6b). After this, we perform a
Taylor expansion of the Volkov action $S_{V}(t,t_{s}^{\prime}(t))$ at the
harmonic-cutoff time $t_{\mathrm{hc}}$ up to third order,

$\displaystyle\mathbf{D}(\Omega)$
$\displaystyle=\sum_{s,\mathrm{hc}}\int_{C}\\!\\!\textrm{d}t\sqrt{\frac{2\pi/i}{\frac{\partial^{2}S_{V}}{\partial
t^{\prime
2}}(t,t_{s}^{\prime}(t))}}\mathbf{f}(t,t_{s}^{\prime}(t))e^{-iS_{V}\\!(t_{\mathrm{hc}},t_{s}^{\prime}(_{\mathrm{hc}}))}$
(43)
$\displaystyle\qquad\times\exp(-\frac{i}{3}A(t-t_{\mathrm{hc}})^{3}-i\Omega_{\mathrm{hc}}(t-t_{\mathrm{hc}})+i\Omega
t),$

where $A$ and $\Omega_{\mathrm{hc}}$ are the third and first derivatives of
$S_{V}(t,t_{s}^{\prime}(t))$, as set in (16) and evaluated as constructed in
Section I.6; here we have also allowed the integration contour $C$ to vary so
that it can pass through $t_{\mathrm{hc}}$ as appropriate. Finally, we move
the integration origin to $t_{\mathrm{hc}}$,

$\displaystyle\mathbf{D}(\Omega)$
$\displaystyle=\sum_{s,\mathrm{hc}}\int_{\bar{C}}\\!\\!\textrm{d}\bar{t}\sqrt{\frac{2\pi/i}{\frac{\partial^{2}S_{V}}{\partial
t^{\prime
2}}(t,t_{s}^{\prime}(t))}}\mathbf{f}(t,t_{s}^{\prime}(t))e^{-iS_{V}\\!(t_{\mathrm{hc}},t_{s}^{\prime}(_{\mathrm{hc}}))}$
(44)
$\displaystyle\qquad\times\exp(-\frac{i}{3}A\bar{t}\,^{3}-i\Omega_{\mathrm{hc}}\bar{t}+i\Omega(\bar{t}+t_{\mathrm{hc}})),$

keeping the explicit $t=\bar{t}+t_{\mathrm{hc}}$ for notational simplicity.

To apply the Harmonic-Cutoff Approximation, we now assume that the prefactor
varies slowly enough at $t_{\mathrm{hc}}$ that it can be pulled out of the
integral (i.e., that it changes slower than the decay into the valleys shown
in Fig. 7),888On a more rigorous footing, this amounts to a zeroth-order
Taylor expansion of the prefactor at $t_{\mathrm{hc}}$, as is done in the SPA
GerlachSPAonline . Derivatives of the prefactor can be added as required, and
the corresponding terms will change the Airy function to its derivatives (as
can be seen by differentiating (36) with respect to $z$, which brings down a
factor of $t$ into the prefactor). For the cases we plot here, these terms are
about two orders of magnitude weaker than the leading-order contribution.
which gives

$\displaystyle\mathbf{D}(\Omega)$
$\displaystyle=\sum_{\mathrm{hc}}\sqrt{\frac{2\pi}{i\frac{\partial^{2}S_{V}}{\partial
t^{\prime
2}}(t_{\mathrm{hc}}^{\phantom{a}},t^{\prime}_{\mathrm{hc}})}}\mathbf{f}(t_{\mathrm{hc}}^{\phantom{a}},t^{\prime}_{\mathrm{hc}})e^{-iS_{V}\\!(t_{\mathrm{hc}}^{\phantom{a}},t_{\mathrm{hc}}^{\prime})}$
(45) $\displaystyle\qquad\times e^{i\Omega
t_{\mathrm{hc}}}\int_{\bar{C}}\exp(-\frac{i}{3}A\bar{t}\,^{3}+i(\Omega-\Omega_{\mathrm{hc}})\bar{t})\textrm{d}\bar{t}.$

The integral is now in the form of (35), with an added functional dependence
on $\Omega$ coming from the term $e^{i\Omega t_{\mathrm{hc}}}$. This term
results from the translation of the time origin, and it can have a sizeable
effect on the harmonic yield if $\imaginary(t_{\mathrm{hc}})$ is nonzero. That
said, we can now use the result (40) for the model form directly.

We obtain thus our final result for the HCA,

$\displaystyle\mathbf{D}(\Omega)$
$\displaystyle=\sum_{\mathrm{hc}}\sqrt{\frac{2\pi}{i\frac{\partial^{2}S_{V}}{\partial
t^{\prime
2}}(t_{\mathrm{hc}}^{\phantom{a}},t^{\prime}_{\mathrm{hc}})}}\frac{2\pi}{e^{2\pi
ik/3}A_{\mathrm{hc}}^{1/3}}$ (46)
$\displaystyle\quad\times\mathbf{f}(t_{\mathrm{hc}}^{\phantom{a}},t^{\prime}_{\mathrm{hc}})e^{-iS_{V}\\!(t_{\mathrm{hc}}^{\phantom{a}},t_{\mathrm{hc}}^{\prime})+i\Omega
t_{\mathrm{hc}}}\operatorname{Ai}\mathopen{}\left(\frac{\Omega_{\mathrm{hc}}-\Omega}{e^{2\pi
ik/3}A_{\mathrm{hc}}^{1/3}}\right)\mathclose{},$

where the harmonic-cutoff time pairs
$(t_{\mathrm{hc}}^{\phantom{a}},t^{\prime}_{\mathrm{hc}})$ are found by
solving the simultaneous equations

$\displaystyle\frac{\partial^{2}S_{V}}{\partial
t^{2}}\frac{\partial^{2}S_{V}}{\partial t^{\prime
2}}-\left(\frac{\partial^{2}S_{V}}{\partial t\,\partial
t^{\prime}}\right)^{2}$ $\displaystyle=0$ (47a) $\displaystyle\frac{\partial
S_{V}}{\partial t^{\prime}}$ $\displaystyle=0$ (47b)

and where the coefficients $\Omega_{\mathrm{hc}}$ and $A_{\mathrm{hc}}$ for
each pair of harmonic-cutoff times
$(t_{\mathrm{hc}}^{\phantom{a}},t^{\prime}_{\mathrm{hc}})$ are given, as in
(16), by

$\displaystyle\Omega_{\mathrm{hc}}=\frac{\partial S_{V}}{\partial
t}(t_{\mathrm{hc}}^{\phantom{a}},t^{\prime}_{\mathrm{hc}})\ \text{and}\
A_{\mathrm{hc}}=\frac{1}{2}\frac{\textrm{d}^{3}S_{V}}{\textrm{d}t^{3}}(t_{\mathrm{hc}}^{\phantom{a}},t^{\prime}_{\mathrm{hc}}),$
(48)

with the third derivative $\frac{\textrm{d}^{3}}{\textrm{d}t^{3}}$ understood
in the constrained sense of Section I.6.999An added wrinkle can appear if the
prefactor, and particularly the dipole moments it includes, have singularities
at the solutions of (47), which is often the case in the regular SPA. If this
occurs, a regularization scheme like the one described in Ref.
Popruzhenko2014,  will be required, but this is likely to extend relatively
cleanly.  Our implementation of this approximation, in the Wolfram Language,
is available in the RBSFA software package RBSFA .

---

(a)

(b)

(c)

Figure 9:  Approximations to (a) the pure action term $e^{-iS}$, (b) the
action term with wavepacket diffusion, $\tau^{-3/2}e^{-iS}$, and (c) the full
harmonic dipole, via the Saddle-Point Approximation for the long and short
trajectories (green and blue lines) and their coherent combination (dashed
gray lines), the Uniform Approximation (black lines), and the Harmonic-Cutoff
Approximation (red lines). The HCA is quantitatively accurate at the harmonic
cutoff, and reasonably qualitatively accurate in the upper plateau, but it
only requires solving a single saddle-point equation, instead one per harmonic
frequency over a tight grid in $\Omega$.

We show the results of this approximation, compared to the SPA and the UA (as
described in Refs. Figueira2002, ; Milosevic2002, ), in Fig. 9, plotted as in
Fig. 2c. When the prefactor is fully ignored, as in Fig. 9a, the HCA is
extremely accurate at the harmonic cutoff, and it retains a good qualitative
accuracy as $\Omega$ descends from $\Omega_{c}$ and into the plateau: it
shifts vertically from the UA in the lower plateau and the predicted period
for the QPI beatings is too long, but the QPI contrast is mostly well
reproduced.

The qualitative accuracy of the HCA, particularly regarding the QPI contrast,
gets degraded when the wavepacket-diffusion dilution factor of $\tau^{-3/2}$
is included, as shown in Fig. 9b. This happens because the $\tau^{-3/2}$
factor affects the long trajectories (green line) much more than it does the
short ones (blue line), bringing them closer together and allowing for better
interference between them; the HCA, with its information coming from a single
point, is blind to this effect, which only acts as a global vertical shift on
the curve. The same is true for the full harmonic dipole
$\mathbf{D}(\Omega)$,101010Here we use the ground-state wavefunction and
dipole transition matrix element for the $s$ state of a short-range-potential,
as described in Ref. Pisanty2017, .  which we plot in Fig. 9c – the HCA has
good quantitative accuracy at the cutoff, but it becomes more of a qualitative
estimation for the middle and lower plateau.

However, this weakness of the HCA is also its biggest strength: precisely
because it is able to describe the harmonic yield using only information
coming from a single solution of a set of saddle-point equations (47), it can
be calculated in a small fraction of the time taken to produce SPA or UA
calculations of the harmonic spectrum, since those require the solution of one
system of saddle-point equations – the system in (6) – for each harmonic
frequency $\Omega$ of interest, and these will typically form a grid with
hundreds of instances.

In this sense, the HCA sits at the far end of the tradeoff spectrum between
numerical accuracy and computational complexity, which is typically understood
as running from full simulations of the time-dependent Schrödinger equation
(TDSE) Scrinzi2014 , passing through explicit numerical integration of the SFA
dipole (5), to the SPA and UA quantum-orbit approaches. Normally, the quantum-
orbit methods are considered to be fast enough that no further optimization is
necessary, but this is not always the case when large scans are required
spanning a high-dimensional parameter space—particularly if the waveforms
involved cause nontrivial behaviour in the quantum orbits—, such as when
optimizing the length or strength of the harmonic plateau and cutoff
Chipperfield2009 ; Chipperfield2010 ; Haessler2014 . Moreover, when the
primary focus is on the harmonic cutoff, the HCA will typically be
sufficiently accurate.

##  III Parameter scaling and the cutoff law

As we have seen, the formalism we have presented allows us to extract both a
unique and natural definition for the (complex) times that correspond to the
harmonic cutoff, and also a natural definition of the harmonic frequency
$\Omega_{c}$ at which the cutoff occurs. This identification is a valuable
connection, since the high-harmonic cutoff is often understood as a vague
term, referring to a spectral region where the behaviour changes, instead of a
concrete point—largely because of the difficulty in pinning down a specific
frequency at which this change occurs.

This is particularly important, since the scaling of the cutoff frequency with
the field parameters is one of the key hallmarks of HHG, in terms of the so-
called ‘cutoff law’,

$\Omega_{\mathrm{max}}\approx 3.17\,U_{p}+I_{p},\phantom{\,F(I_{p}/U_{p})}$
(49)

which relates the cutoff frequency $\Omega_{\mathrm{max}}$ to the target’s
ionization potential $I_{p}$ and the ponderomotive energy
$U_{p}=F^{2}/4\omega^{2}$ of the field. This relationship was first uncovered
in numerical simulations Krause1992  and subsequently explained using the
classical three-step model Corkum1993 ; Kulander1993 : the numerical factor of
$3.17$ arises as the maximal kinetic energy that can be achieved at the return
to the origin by an electron released at zero velocity, with the $I_{p}$ term
representing quantum energy conservation at the recombination, with an
essentially heuristic justification.

The fully-quantum quasiclassical theory Lewenstein1994  re-derives this cutoff
law, giving also a quantum correction to the $I_{p}$ term of the form

$\Omega_{\mathrm{max}}\approx
3.17\,U_{p}+1.32\,I_{p}.\phantom{\,F(I_{p}/U_{p})}$ (50)

To obtain this formulation, the cutoff frequency is defined as the maximum of
the recollision kinetic energy, $\real(\mathbf{v}(t_{r})^{2})$, taken under
the restriction that its imaginary part $\imaginary(\mathbf{v}(t_{r})^{2})$
vanish. The dynamics are then analyzed using a systematic expansion on
$I_{p}/U_{p}$, using the purely-classical case at $I_{p}=0$ for the leading
$3.17U_{p}$ term, after which the cutoff law can be derived in the exact form

$\Omega_{\mathrm{max}}=3.17\,U_{p}+F(I_{p}/U_{p})\,I_{p},$ (51)

in terms of a universal function $F(I_{p}/U_{p})$ which can be calculated
numerically as $F(0)=1.32$.

This prescription can be computed exactly for linearly-polarized monochromatic
fields, and it can be extended to elliptical polarizations Milosevic2000 ;
Flegel2004 , but it is laborious to apply for polychromatic combinations
Figueira1999  (where the complex waveform makes analytical calculations
challenging) and for tailored field polarizations Milosevic1996 , especially
for field shapes whose recolliding orbits do not ionize at zero velocity, as
required by the classical theory. In any case, this definition of the harmonic
cutoff has not seen wide adoption in the literature.

Instead, a wide variety of other methods have been used in its place,
including relaxing the $\imaginary(\mathbf{v}(t_{r})^{2})=0$ condition to
$\imaginary(t_{r})=0$ Milosevic2000bicircular , focusing only on the classical
recollision velocity Chipperfield2009 ; Chipperfield2010 , examining the
change in the intensity scaling of the harmonic yield at a fixed harmonic
order Lewenstein1995 , the use of graphical methods based on tangents to the
optical waveform Figueira1999 , and extracting the cutoff from the Bessel-
function expansion of the oscillatory integral Averbukh2001 ; Avetissian2014 ,
as well as direct numerical methods used to extract the cutoff from the
results of TDSE simulations Neufeld2019timescales . None of these definitions,
however, is particularly satisfactory, and – like the original definition –
none of the analytical approaches have been used very widely in the
literature.

---

(a)

(b)

(c)

(d)

Figure 10:  Scaling of the harmonic-cutoff frequency $\Omega_{\mathrm{hc}}$
with respect to the field parameters. (a) The affine dependence of the cutoff
law $\real(\Omega_{\mathrm{hc}})\approx 3.17U_{p}+I_{p}$ is reproduced well;
the dots at the axis mark the position of $I_{p}$, and the mismatch to the
lines comes from the quantum corrections. (b) The imaginary part
$\imaginary(\Omega_{\mathrm{hc}})$ decays with $U_{p}$ as
$I_{p}^{3/2}/U_{p}^{1/2}$. These scalings indicate a relationship of the form
$\Omega_{\mathrm{hc}}=3.17U_{p}+F(\gamma)I_{p}$, with the real and imaginary
parts of $F(\gamma)$, shown in (c) and (d), suggesting that it is a power
series in $i\gamma$.

As we argued earlier, the real part $\Omega_{c}$ of our complex-valued
$\Omega_{\mathrm{hc}}$ forms a natural candidate as a precise definition of
the harmonic-cutoff frequency $\Omega_{\mathrm{max}}$, with the added
advantage that it can be trivially adapted to complicated waveforms and to
tailored polarizations without significantly complicating the calculation, a
useful feature as the complexity of the field shapes under consideration
continues to increase Neufeld2019 ; Javsarevic2020 . To strengthen this
identification, though, it is important to verify that it reproduces the
cutoff law, including the quantum corrections coming from the SFA as captured
by the original definition. We show this in Fig. 10a: the variation of
$\real(\Omega_{\mathrm{hc}})$ with $U_{p}$ is linear with an $I_{p}$-dependent
offset, which closely follows the quasiclassical quantum-corrected cutoff law
(50).

However, our formalism allows us to go beyond this level, since
$\Omega_{\mathrm{hc}}$ also has an imaginary part, whose scaling with $U_{p}$
and $I_{p}$ is shown in Fig. 10b. As rough behaviour,
$\imaginary(\Omega_{\mathrm{hc}})$ decreases with $U_{p}$ and increases with
$I_{p}$, and simple testing shows that the leading asymptotic component is the
behaviour

$\imaginary(\Omega_{\mathrm{hc}})\sim{I_{p}^{3/2}}\big{/}{U_{p}^{1/2}}.$ (52)

At first glance, this shows very different behaviour to the real part, which
follows (51). However, once a factor of $I_{p}$ has been set apart, as we did
for the real part, the scaling for $\imaginary(\Omega_{\mathrm{hc}})$ is
simply the Keldysh parameter, $\gamma=\sqrt{I_{p}/2U_{p}}$. This is therefore
indicative that the parameter in the quantum scaling law (51) should be
amended from $I_{p}/U_{p}$ to its square root, so it should read

$\Omega_{\mathrm{hc}}=3.17\,U_{p}+F(\gamma)\,I_{p}.$ (53)

In this form, (53) essentially acts as a definition for
$F(\gamma)=(\Omega_{\mathrm{hc}}-3.17U_{p})/I_{p}$, but this definition can
only work if that value is independent of what combination of $I_{p}$ and
$U_{p}$ gives rise to $\gamma$. This is indeed the case, as we show in Figs.
10c and 10d for the real and imaginary parts, respectively: plotting
$F(\gamma)$ by changing $I_{p}$ reveals the same universal curve for a range
of different values of $U_{p}$. Moreover, we can extract the low-$\gamma$
behaviour here to obtain

$\displaystyle F(\gamma)$ $\displaystyle\approx 1.323+i\,0.07361\,\gamma$ (54)
$\displaystyle\quad\ -0.068\,\gamma^{2}-i\,0.025\,\gamma^{3}+\cdots,$

with the first two terms shown as the gray dashed line in Figs. 10c and 10d.
This understanding coincides with the existing theory (so, in particular, Fig.
10c coincides with Fig. 5 of Ref. Lewenstein1994,  when plotted as a function
of $\gamma^{2}$), but it also helps extend it and clarify its structure,
though a formal analysis of the existence and convergence of the power series
in $i\gamma$ for $F(\gamma)$, as defined here for $\Omega_{\mathrm{hc}}$, is
still lacking. We summarize our parameter-scaling results, and their
relationship to the existing theories, in Table 1.

Theory |  | Cutoff law
---|---|---
Corkum Corkum1993  |  | $\Omega_{\mathrm{max}}\approx 3.17U_{p}+I_{p}$
Kulander Kulander1993  |
Lewenstein Lewenstein1994  |  | $\begin{aligned} \Omega_{\mathrm{max}}&=3.17U_{p}+F(I_{p}/U_{p})I_{p}\\\ F(0)&=1.32\end{aligned}$
This work |  | $\begin{aligned} \Omega_{\mathrm{hc}}&=3.17U_{p}+F(\gamma)I_{p}\\\ F(\gamma)&=1.323+i\,0.07361\,\gamma\\\ &\quad\ -0.068\,\gamma^{2}-i\,0.025\,\gamma^{3}+\cdots\end{aligned}$
Table 1:  Scaling of the harmonic cutoff with $I_{p}$ and $U_{p}$, for the
different understandings of the cutoff.

##  IV Nontrivial quantum-orbit topology for bicircular fields

To conclude our (current) explorations of the role of the complex harmonic-
cutoff times in HHG, we return to the classification scheme for separating
saddle-point solutions into quantum orbits, which we first showcased in Fig. 1
and constructed in detail in Section I.4. One crucial aspect of this
classification scheme is the orientation of the separatrix, which indicates
the direction that the saddle points take in their avoided crossing at the
cutoff: i.e., whether the short-trajectory saddle, coming in from the left,
will go up into positive (or down into negative) imaginary time when it meets
the long-trajectory saddle. Our construction shows that this direction is
given, via (23), by

$\displaystyle\delta t_{\mathrm{sep}}=\sqrt{-i\eta/A_{\mathrm{hc}}},$ (55)

and thus that it has a sensitive dependence on the sign of

$\eta=\frac{\partial S_{V}}{\partial
t}(t_{\mathrm{hc}}^{\phantom{a}},t_{\mathrm{hc}}^{\prime}),$ (56)

the imaginary part of the time derivative of the Volkov action at the
harmonic-cutoff time.

This imaginary part is essentially controlled by the internal details of the
action and, as such, its sign can be either positive or negative depending on
the precise particulars of the optical waveform of the driving laser. If the
shape of the driving laser pulse remains essentially constant (say, when doing
an intensity scan as in Fig. 10) then the sign of the imaginary part is
unlikely to change (as was seen to be the case in Fig. 10b).

However, if one has a more complicated driving field with a nontrivial shape
parameter that affects the waveform of the pulse, and thus the details of the
Volkov action, then the sign of $\eta$ will change, and the direction of the
separatrix, $\delta t_{\mathrm{sep}}$, will switch to one
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$ away. This will, in turn, switch the
direction of the avoided crossing, and with that the identity (say, as ‘short’
or ‘long’ trajectories) of the post-cutoff branches of the orbit. This
behaviour is generic and commonplace: as a simple example, it can be observed
in monochromatic fields once an envelope is introduced, where changes in the
carrier-envelope phase can produce this type of topological change in the
quantum-orbit layout.

In this section we showcase a concrete example of this behaviour, and we
explore the nontrivial topologies that it induces in the set of quantum-orbit
saddle points. We focus, in particular, on tailored polarization states known
as ‘bicircular’ fields Fleischer2014 ; Kfir2015 ; Eichmann1995 ; Milosevic1996
; Milosevic2000bicircular ; Baykusheva2016 ; Dorney2017 ; JimenezGalan2018 ;
Pisanty2014 ; Pisanty2017 ; Milovsevic2019 , formed by the combination of two
counter-rotating circularly-polarized strong laser pulses at frequencies
$\omega$ and $2\omega$,

$\mathbf{F}(t)=F_{1}\begin{pmatrix}\cos(\omega t)\\\ \sin(\omega
t)\end{pmatrix}+F_{2}\begin{pmatrix}\phantom{-}\cos(2\omega t)\\\
-\sin(2\omega t)\end{pmatrix}.$ (57)

These fields have been the focus of widespread interest in recent years
because they are subject to a spin selection rule Fleischer2014 ; Pisanty2014
which ensures that the harmonics of the driver are circularly polarized
Kfir2015 . For our purposes, however, we select them as an example of a
reasonably well-understood waveform with a nontrivial shape parameter, namely,
the relative intensity between the two pulses; this is normally kept at unity,
but the effects of its variation have seen some exploration
Milosevic2000bicircular ; Baykusheva2016 ; Dorney2017 ; JimenezGalan2018 . If
we keep the total intensity constant, at
$I_{\mathrm{tot}}=$2\text{\times}{10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$$
for concreteness, then the shift in intensity is best described by the mixing
angle $\theta$, which we use to define the individual field amplitudes as

$\displaystyle F_{1}=F\cos(\theta),\qquad F_{2}=F\sin(\theta),$ (58)

where $F=\sqrt{2}\times$0.053\text{\,}\mathrm{{a.u.}}$$ is the peak field
strength.

---

(a)

(b)

(c)

Figure 11:  Topological transition and reconnection in the quantum orbits of a
bicircular field as the mixing angle $\theta$ changes. Here we show a detail
of the quantum-orbit layout for the short and long trajectories, plotted as in
Fig. 1c, with the triangles marking the harmonic-cutoff times
$t_{\mathrm{hc}}$ and the coloured separatrices marking the regions occupied
by the different orbits, as described in Section I.4. In (b) we show the
topological transition, where
$\eta=\imaginary(\tfrac{\textrm{d}S_{V}}{\textrm{d}t})$ vanishes and the
separatrix becomes undefined, which we show by highlighting the
$t_{\mathrm{hc}}$ with a black border and plotting both options (using an
artificial $\eta=\pm 1$) as the gray cross. In this case, both choices of
separatrix are acceptable for saddle classification, since the saddles fully
coalesce.

We show in Fig. 11 a representative topological reconnection transition,
between the short and long orbits (blue and green curves, respectively) of the
bicircular field. Before the transition, in Fig. 11a, the short orbit connects
up to the positive-imaginary saddle, while after the transition, in Fig. 11c,
it connects to negative imaginary time. Between these two there must always
lie a transition, which we show in Fig. 11b: here the saddles fully coalesce,
producing a monkey-saddle landscape like the ones shown in Fig. 7, with a
purely-real HCA Airy function corresponding to enhanced quantum path
interference between the two orbits.

|
---
|
---

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 12: Topological transitions and reconnections in the quantum orbits of
a bicircular field as a function of the mixing angle $\theta$, plotted as in
Fig. 11. As $\theta$ varies over a wider range, the various pairs of
neighbouring quantum orbits go through reconnection transitions with each
other (with two separate transitions, in (b) and (f), for the first pair of
orbits), so that all of the quantum orbits are connected into a single
topology, which we show below in 3D Fig. 2. At each panel the inset at bottom
right shows the shape of the electric field at the relevant mixing angle.

More importantly, however, this transition separates the instances with
$\eta>0$ from those with $\eta<0$, so at the transition itself $\eta$ must
vanish – and, as with the exact coalescences we studied in Section I.7, no
(nonzero) unique separatrix direction $\delta t_{\mathrm{sep}}$ can be found.
That said, this failure is fairly benign since, at the transition, either of
the separatrices obtained by the artificial choices of $\eta=\pm 1$ will work
correctly. At the transition the missed approach becomes a full coalescence,
and any identification of pre- and post-cutoff orbits is artificial, so both
signs of $\eta$ will work equally well. In Fig. 11b, we mark the transition by
showing both possible separatrices in gray, taking an arbitrary choice between
the two as to which one to use for the actual classification. We also
highlight the triangle at the harmonic-cutoff time $t_{\mathrm{hc}}$ with a
black border for further clarity.

On a more practical footing, the fact that
$\eta=\imaginary(\tfrac{\textrm{d}S_{V}}{\textrm{d}t})$ vanishes at the
transition is especially useful, since it can be used to look for the
parameters where the transitions happen, by looking for zeros of $\eta$ as a
function of $\theta$.

That said, the most important aspect of these topological transitions becomes
apparent when we survey the quantum-orbit landscape on a wider perspective, in
terms of the multiple quantum orbits present in the dynamics, as well as for a
wider interval in the mixing angle, which we show in Fig. 12. Here we see that
all of the pairs of quantum orbits undergo one or more reconnection
transitions with each other, so that no pair of quantum orbits, however
apparently distant, can actually be separated: in the same way that the
evanescent orbits of Fig. 11b have ‘confused’ identities between the short and
long orbits, the short orbit ‘mingles’ at low $\Omega$ with the short-time
orbit shown in cyan, at two separate transitions (shown in Fig. 12b and Fig.
12f), the long orbit connects with the second return (with the transition
shown in Fig. 12d), and so on.

---
List of tdfigures 2  Unified surface formed by the quantum orbits in Fig. 12,
when the $\theta$ dependence is unfolded as a third dimension. (Thus, a cut
through the surface at constant $\theta$ will return a panel from Fig. 12.) At
the transitions, the horizontal parts of the surface (the plateau harmonics)
change from going up to going down and vice versa, but from a global
perspective, the saddles form a single surface with a unified topology, and
the separation of this surface into individual regions corresponding to the
different quantum orbits is somewhat artificial. This highlights the fact that
the various quantum orbits, as inverse images under the multivalued inverse of
an analytical function, are simply different branches of a single Riemann
surface. An interactive version and a 3D-printable model of this plot is
available at imaginary-harmonic-cutoff.github.io SupplementaryMaterial .

|
---|---
|
|
|

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 13:  Energy-time relations (a, c, e, g) and indicative SPA harmonic
spectra, including the factors from the action and the wavepacket dilution (b,
d, f, h) for four of the quantum-orbit transitions from Fig. 12. For
bicircular fields, attention is generally focused on the short quantum orbit
(blue curve), since at equal intensities
($\theta=$45\text{\,}\mathrm{\SIUnitSymbolDegree}$$) it dominates the spectrum
(as shown in (b), which closely follows Fig. 8 of Ref.
Milosevic2000bicircular, ). However, this is no longer the case at mixing
angles with a high intensity for the $2\omega$ field. At the long-short
topological transition (e, f), where the two become indistinguishable at the
cutoff, the latter is intensified by this effect, while at large mixing
angles, as shown in (g, h), the contribution of the short trajectory plummets.

In other words, the quantum orbits here form a single unified topology, and
they should be understood as such. We show this topology in 3D Fig. 2, by
unfolding the $\theta$ dependence into a third dimension for the plot: this
reveals that the quantum orbits form a single surface, with discontinuous
changes in colour where the (apparent) reconnection transitions force us to
choose where to split this unified surface into individual components.

The fundamental topological object here, as we explored in depth in Section
I.5, is the Riemann surface formed by the saddle points: that is, the surface
defined by the equation $\frac{\textrm{d}S_{V}}{\textrm{d}t}(t_{s})=\Omega$,
encoding the multiple inverses of the action’s derivative, whose intersection
with the real $\Omega$ axis gives the quantum orbits themselves. Within that
perspective, the difference between the two types of connections between the
quantum orbits boils down to what side of the real $\Omega$ axis the branch
point (i.e., the harmonic-cutoff time $t_{\mathrm{hc}}$) falls on.

For fixed field shapes, this branch-cut structure can essentially be ignored,
if so desired, since we are only looking for the inverse images for $\Omega$
over the real axis, and those inverse images will normally be well separated.
In 3D Fig. 2, however, as well as in any other situations where the field
shape changes over a control parameter, we see how these changes to the
position and shape of the Riemann surface bring these once-separate images
into collision (and reconnection) with each other, and this structure can no
longer be ignored.

In more general terms, the reconnection phenomenon we have demonstrated in
this section entails that, ultimately, attempting to attach labels to
individual quantum orbits is fundamentally impossible, since they are
essentially just different instances of one and the same object. Saddle-point
tagging schemes of this type, in terms of multi-indices typically notated as
$\alpha\beta m$, have been used widely in the literature, both for
monochromatic fields Chipperfield2007 ; Milosevic2002 ; Odzak2005 ;
Chipperfield2005 ; Milosevic2017chapter  as well as more complex waveforms
Hasovic2016 ; Milosevic2019xray ; Milovsevic2019quantum ; Milosevic2016 . The
topological features we have explored here imply that these schemes should be
regarded as labelling the principal branches of the quantum-orbit Riemann
surface.

These schemes can, in principle, be applied to changing driver waveforms with
nontrivial connections between the various branches involved, but then those
branch changes must be explicitly tracked – in which case, the branch-point
machinery we have developed here (both for finding the branch points
$t_{\mathrm{hc}}$ as well as the parameters where their branch cuts cross the
real $\Omega$ axis) is essential.

As a final note here, it is important to remark that, while some of the
dynamics we have explored here occurs for quantum orbits with weak or
negligible contributions to the spectrum, that is not always the case. Indeed,
as we show in Fig. 13, the equal-intensities viewpoint that the spectrum is
completely dominated by the short-trajectory contribution breaks down at large
mixing angles, where the $2\omega$ contribution is significant, and where many
of the topological transitions we have discussed in this section take place.

##  V Discussion and conclusions

As we have shown, the problem of saddle-point classification can be solved in
a robust and flexible way by using the harmonic-cutoff times $t_{\mathrm{hc}}$
to locate the center of the missed approach and implement a separatrix that
lies between the two quantum orbits that perform an avoided crossing at each
cutoff. These harmonic-cutoff times can be easily found as the zeros of a
second-order saddle-point equation. In a clear sense, the $t_{\mathrm{hc}}$’s
form the centerpiece of the quasiclassical structure at the cutoff, and they
are a useful and flexible tool to understand a range of questions about the
quantum-orbit dynamics.

The strongest implication of this is that the harmonic-cutoff times can be
used to provide a natural identification for the energy of the cutoff,
$\Omega_{\mathrm{hc}}=\frac{\partial S_{V}}{\partial t}(t_{\mathrm{hc}})$.
This now takes on a complex value, with the imaginary part controlling the
strength of quantum-path interference between the quantum orbits that meet at
that cutoff, as well as the closeness and direction of the approach between
them.

Similarly, our approach provides an efficient method for finding the position
of the cutoff as well as the harmonic yield there, which will work uniformly
and efficiently for a broad range of optical waveforms for the driving laser.
This extends to a full estimation of the spectrum, the Harmonic-Cutoff
Approximation, which uses only quantities local to $t_{\mathrm{hc}}$, and
which accurately captures the cutoff as well as the qualitative shape of the
spectrum down into the middle of the plateau.

On a more abstract note, our search for structures that can be used to
classify the solutions of the saddle-point equations into individual quantum
orbits yields a fresh perspective on the quasiclassical theory of strong-field
phenomena. The quantum orbits are thus revealed as the
$\imaginary(\frac{\partial S_{V}}{\partial t})=0$ contour of the time
derivative of the action, with the rest of its contour map holding crucial
information – most notably via its saddles, the harmonic-cutoff times
$t_{\mathrm{hc}}$. Likewise, the usual saddle points are re-understood as
individual branches of a larger Riemann surface, which encompasses all of the
quantum orbits. The branch points that separate these branches are again the
harmonic-cutoff times, which can thus be used to watch for changes in the
quantum-orbit topology as the driving field’s waveform changes.

The location of the cutoff at a zero of the second derivative of the action
also admits a rather more pedestrian interpretation, though. If the emission
is modelled using only the classical-trajectory level, then the energy of
return $E_{\mathrm{return}}=\frac{\partial S}{\partial t}$ is a function of
the return time, and if we want to find its extrema, then we simply need to
solve the equation $\frac{\partial^{2}S}{\partial t^{2}}=0$. When the theory
is upgraded to the quasiclassical formalism, that equation might seem to lose
its meaning, since the action and the trajectories are complex, as would be
the solutions to the extremum equation $\frac{\partial^{2}S}{\partial
t^{2}}=0$. Our solution in this work is entirely in line with the spirit of
the quasiclassical formalism for strong-field physics: we solve this extremum
equation in the same form it takes classically, allowing for complex values
wherever necessary, and then look at the underlying oscillatory integral for
the correct interpretation of those complex quantities.

More physically, the key structure at play is the fact that the harmonic
cutoff is a caustic: the quantum-orbit analysis of the underlying matter-wave
dynamics is exactly analogous to the geometric-optics analysis of wave optics
in terms of rays, and in this analogy the harmonic cutoff corresponds to
caustics where two families of rays interfere and then meet and fold into each
other, giving way to evanescent-wave behaviour analogous to the post-cutoff
decay in the harmonic spectrum. In one dimension, this caustic can be
precisely localized to a point, known as the ‘fold’ (more technically, the
‘bifurcation set’), which marks the transition between the two regimes, and
our harmonic-cutoff times are the precise embodiment of that understanding of
the caustic.

This view of the harmonic cutoff as a caustic has been used for some time
Raz2012 ; Goulielmakis2012 ; Austin2014 , but recent years have seen a marked
increase in interest in that perspective on the dynamics Facciala2016 ;
Facciala2018 ; Hamilton2017 ; Birulia2019 ; Uzan2020 ; Kelvich2017 , as laser
sources achieve better control over polychromatic combinations with high
relative intensities (which is slated to continue to increase Mitra2020 ).
This has opened the door to the observation of higher-order catastrophes in
strong-field observables, involving higher-dimensional bifurcation sets;
however, the theoretical understanding remains somewhat behind, and a
quantitative understanding that reflects those structures has not yet
followed. In this work we have precisely pinned down the embodiment of the
bifurcation set in the quasiclassical formalism for the ‘fold’ catastrophe,
and this can then be used to analyze in detail the more complicated
configurations involved in newer experiments.

Similarly, our work suggests that the harmonic-cutoff times we have discussed
here should have analogous structures in ionization experiments (most notably
high-order above-threshold ionization Figueira2019 ), where the
experimentally-measurable variable, momentum, has multiple dimensions. This
allows greater freedom to the theory (while at the same time substantially
complicating its analysis), and this in turn allows for higher-dimensional
singular matter-wave structures to be contained in the experimental results.
The tools we have demonstrated here for understanding the one-dimensional
caustic formed at the harmonic cutoff should then provide a useful basis for
understanding those configurations.

###### Acknowledgements.

We thank Przemyslaw Grzybowski for essential and generous help in
understanding the structures presented here. E.P. acknowledges Cellex-ICFO-MPQ
Fellowship funding. We acknowledge funding from the Spanish Ministry MINECO
(National Plan 15 Grant: FISICATEAMO No. FIS2016-79508-P, SEVERO OCHOA No.
SEV-2015-0522, FPI), European Social Fund, Fundació Cellex, Fundació Mir-Puig,
Generalitat de Catalunya (AGAUR Grant No. 2017 SGR 1341, CERCA/Program), ERC
AdG NOQIA, EU FEDER, European Union Regional Development Fund – ERDF
Operational Program of Catalonia 2014-2020 (Operation Code: IU16-011424),
MINECO-EU QUANTERA MAQS (funded by The State Research Agency (AEI)
PCI2019-111828-2 / 10.13039/501100011033), and the National Science Centre,
Poland, Symfonia Grant No. 2016/20/W/ST4/00314.

## Author ORCIDs

  * Emilio Pisanty: 0000-0003-0598-8524

  * Marcelo F. Ciappina: 0000-0002-1123-6460

  * Maciej Lewenstein: 0000-0002-0210-7800

## References

  * (1) F. Krausz and M. Ivanov.  Attosecond physics. _Rev. Mod. Phys._ 81 no. 1, pp. 163–234 (2009).
  * (2) P. B. Corkum and F. Krausz.  Attosecond science. _Nat. Phys._ 3 no. 6, pp. 381–387 (2007).
  * (3) P. B. Corkum.  Plasma perspective on strong field multiphoton ionization. _Phys. Rev. Lett._ 71 no. 13, pp. 1994–1997 (1993).
  * (4) K. C. Kulander, K. J. Schafer and J. L. Krause.  Dynamics of short-pulse excitation, ionization and harmonic conversion.  In B. Piraux, A. L’Huillier and K. Rzążewski (eds.), _Super-Intense Laser Atom Physics_ , vol. 316 of _NATO Advanced Studies Institute Series B: Physics_ , pp. 95–110 (Plenum, New York, 1993).
  * (5) A. Scrinzi.  Time-dependent Schrödinger equation.  In T. Schultz and M. Vrakking (eds.), _Attosecond and XUV Physics: Ultrafast Dynamics and Spectroscopy_ , pp. 257–292 (Wiley-VCH, Weinheim, 2014).
  * (6) 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 no. 3, pp. 2117–2132 (1994).
  * (7) K. Amini, J. Biegert, F. Calegari et al.  Symphony on strong field approximation. _Rep. Progr. Phys._ 82 no. 11, p. 116001 (2019).  arXiv:1812.11447.
  * (8) P. Salières, B. Carré, L. Le Déroff et al.  Feynman’s path-integral approach for intense-laser-atom interactions. _Science_ 292 no. 5518, pp. 902–905 (2001).
  * (9) G. G. Paulus, F. Grasbon, A. Dreischuh et al.  Above-threshold ionization by an elliptically polarized field: Interplay between electronic quantum trajectories. _Phys. Rev. Lett._ 84 no. 17, pp. 3791–3794 (2000).
  * (10) R. Kopold, W. Becker and D. B. Milošević.  Quantum orbits: a space-time picture of intense-laser-induced processes in atoms. _J. Mod. Opt._ 49 no. 12, pp. 1987–1999 (2002).
  * (11) M. Ivanov and O. Smirnova.  Multielectron high harmonic generation: simple man on a complex plane.  In T. Schultz and M. Vrakking (eds.), _Attosecond and XUV Physics: Ultrafast Dynamics and Spectroscopy_ , pp. 201–256 (Wiley-VCH, Weinheim, 2014).  arXiv:1304.2413.
  * (12) A. Nayak, M. Dumergue, S. Kühn et al.  Saddle point approaches in strong field physics and generation of attosecond pulses. _Phys. Rep._ 833, pp. 1–52 (2019).  U Szeged eprint.
  * (13) N. G. de Bruijn.  _Asymptotic methods in analysis_ (Dover Publications, New York, 1981).
  * (14) N. Bleistein and R. A. Handelsman.  _Asymptotic Expansions of Integrals_ (Dover, New York, 1975).
  * (15) U. H. Gerlach.  _Linear Mathematics in Infinite Dimensions: Signals, Boundary Value Problems and Special Functions_ (Ohio State University, 2017), pp. 370–380.  2nd Beta Edition. Available online at http://people.math.osu.edu/gerlach.1/math/BVtypset. Retrieved 12 February 2019.
  * (16) M. Lewenstein, P. Salières and A. L’Huillier.  Phase of the atomic polarization in high-order harmonic generation. _Phys. Rev. A_ 52 no. 6, pp. 4747–4754 (1995).
  * (17) A. Zaïr, M. Holler, A. Guandalini et al.  Quantum path interferences in high-order harmonic generation. _Phys. Rev. Lett._ 100 no. 14, p. 143902 (2008).
  * (18) C. Figueira de Morisson Faria, H. Schomerus and W. Becker.  High-order above-threshold ionization: The uniform approximation and the effect of the binding potential. _Phys. Rev. A_ 66 no. 4, p. 043413 (2002).  UCL eprint.
  * (19) D. B. Milošević and W. Becker.  Role of long quantum orbits in high-order harmonic generation. _Phys. Rev. A_ 66 no. 6, p. 063417 (2002).
  * (20) L. E. Chipperfield.  _High Harmonic Generation with Few-Cycle Pulses_. PhD thesis, Imperial College London (2007).
  * (21) D. Hoffmann.  _High Harmonic Generation Using Multicolour Fields_. PhD thesis, Imperial College London (2011).
  * (22) T. Das.  _Quantum-orbit analysis of laser-matter interactions in intense orthogonally polarised fields_. PhD thesis, University College London (2017).
  * (23) C. Manzoni, O. D. Mücke, G. Cirmi et al.  Coherent pulse synthesis: towards sub-cycle optical waveforms. _Laser Photonics Rev._ 9 no. 2, pp. 129–171 (2015).  PoliMi eprint.
  * (24) A. Fleischer, O. Kfir, T. Diskin, P. Sidorenko and O. Cohen.  Spin angular momentum and tunable polarization in high-harmonic generation. _Nat. Photon._ 8 no. 7, pp. 543–549 (2014).  arXiv:1310.1206.
  * (25) O. Kfir, P. Grychtol, E. Turgut et al.  Generation of bright phase-matched circularly-polarized extreme ultraviolet high harmonics. _Nat. Photon._ 9 no. 2, pp. 99 – 105 (2015).
  * (26) H. Eichmann, A. Egbert, S. Nolte et al.  Polarization-dependent high-order two-color mixing. _Phys. Rev. A_ 51 no. 5, pp. R3414–R3417 (1995).
  * (27) D. B. Milošević and B. Piraux.  High-order harmonic generation in a bichromatic elliptically polarized laser field. _Phys. Rev. A_ 54 no. 2, pp. 1522–1531 (1996).
  * (28) D. B. Milošević, W. Becker and R. Kopold.  Generation of circularly polarized high-order harmonics by two-color coplanar field mixing. _Phys. Rev. A_ 61 no. 6, p. 063403 (2000).
  * (29) D. Baykusheva, M. S. Ahsan, N. Lin and H. J. Wörner.  Bicircular high-harmonic spectroscopy reveals dynamical symmetries of atoms and molecules. _Phys. Rev. Lett._ 116 no. 12, p. 123001 (2016).  ETH eprint.
  * (30) K. M. Dorney, J. L. Ellis, C. Hernández-García et al.  Helicity-selective enhancement and polarization control of attosecond high harmonic waveforms driven by bichromatic circularly polarized laser fields. _Phys. Rev. Lett._ 119 no. 6, p. 063201 (2017).  JILA eprint.
  * (31) A. Jiménez-Galán, N. Zhavoronkov, D. Ayuso et al.  Control of attosecond light polarization in two-color bicircular fields. _Phys. Rev. A_ 97 no. 2, p. 023409 (2018).  arXiv:1805.02250.
  * (32) E. Pisanty, S. Sukiasyan and M. Ivanov.  Spin conservation in high-order-harmonic generation using bicircular fields. _Phys. Rev. A_ 90 no. 4, p. 043829 (2014).  arXiv:1404.6242.
  * (33) E. Pisanty and A. Jiménez-Galán.  Strong-field approximation in a rotating frame: high-order harmonic emission from $p$ states in bicircular fields. _Phys. Rev. A_ 96 no. 6, p. 063401 (2017).  arXiv:1709.00397.
  * (34) D. B. Milošević.  Quantum-orbit analysis of high-order harmonic generation by bicircular field. _J. Mod. Opt._ 66 no. 1, pp. 47–58 (2019).
  * (35) T. Poston and I. N. Stewart.  _Catastrophe Theory and its Applications_ (Pitman, London, 1978).
  * (36) O. Raz, O. Pedatzur, B. D. Bruner and N. Dudovich.  Spectral caustics in attosecond science. _Nat. Photon._ 6 no. 3, p. 170 (2012).
  * (37) E. Goulielmakis.  Attosecond photonics: Extreme ultraviolet catastrophes. _Nat. Photon._ 6 no. 3, p. 142 (2012).
  * (38) D. R. Austin and J. Biegert.  Attosecond pulse shaping using partial phase matching. _New J. Phys._ 16 no. 11, p. 113011 (2014).
  * (39) D. Faccialà, S. Pabst, B. D. Bruner et al.  Probe of multielectron dynamics in xenon by caustics in high-order harmonic generation. _Phys. Rev. Lett._ 117 no. 9, p. 093902 (2016).  handle:11311/996104.
  * (40) D. Faccialà, S. Pabst, B. Bruner et al.  High-order harmonic generation spectroscopy by recolliding electron caustics. _J. Phys. B: At. Mol. Opt. Phys._ 51 no. 13, p. 134002 (2018).
  * (41) K. R. Hamilton, H. W. van der Hart and A. C. Brown.  Pulse-shape control of two-color interference in high-order-harmonic generation. _Phys. Rev. A_ 95 no. 1, p. 013408 (2017).  arXiv:1701.01640.
  * (42) V. A. Birulia and V. V. Strelkov.  Spectral caustic in two-color high-order harmonic generation: Role of Coulomb effects. _Phys. Rev. A_ 99 no. 4, p. 043413 (2019).  arXiv:1901.10518.
  * (43) A. J. Uzan, G. Orenstein, Á. Jiménez-Galán et al.  Attosecond spectral singularities in solid-state high-harmonic generation. _Nat. Photon._ 14 no. 3, pp. 183–187 (2020).  arXiv:1812.02498.
  * (44) S. A. Kelvich, W. Becker and S. P. Goreslavski.  Caustics and catastrophes in above-threshold ionization. _Phys. Rev. A_ 96 no. 2, p. 023427 (2017).
  * (45) E. Pisanty.  RB-SFA: High Harmonic Generation in the Strong Field Approximation via Mathematica.  GitHub, https://github.com/episanty/RB-SFA, v2.1.2 (2016).
  * (46) E. Pisanty.  Figure-maker code for The imaginary part of the high-harmonic cutoff.  Zenodo:3692563, doi:10.5281/zenodo.3692563 (2020).
  * (47) Supplementary Material available online at imaginary-harmonic-cutoff.github.io, or for local use at doi:10.5281/zenodo.3758483 (2020).
  * (48) E. Pisanty.  _Electron dynamics in complex time and complex space_. PhD thesis, Imperial College London (2016).
  * (49) E. Pisanty and M. Ivanov.  Slalom in complex time: emergence of low-energy structures in tunnel ionization via complex time contours. _Phys. Rev. A._ 93 no. 4, p. 043408 (2016).  arXiv:1507.00011.
  * (50) T. Keil, S. V. Popruzhenko and D. Bauer.  Laser-driven recollisions under the Coulomb barrier. _Phys. Rev. Lett._ 117 no. 24, p. 243003 (2016).  arXiv:1608.03844.
  * (51) F. W. J. Olver, D. W. Lozier, R. F. Boisvert and C. W. Clark (eds.).  _NIST Handbook of Mathematical Functions_ (Cambridge University Press NIST, Cambridge, New York, 2010).  Available online as the Digital Library of Mathematical Functions.
  * (52) M. V. Berry and C. Upstill.  Catastrophe optics: morphologies of caustics and their diffraction patterns.  In E. Wolf (ed.), _Progress in Optics_ , vol. 18, chap. IV, pp. 257–346 (Elsevier, 1980).
  * (53) M. Frolov, N. Manakov, T. Sarantseva and A. F. Starace.  Analytic formulae for high harmonic generation. _J. Phys. B: At. Mol. Opt. Phys._ 42 no. 3, p. 035601 (2009).  UNL eprint.
  * (54) M. V. Frolov, N. L. Manakov, A. A. Silaev and N. V. Vvedenskii.  Analytic description of high-order harmonic generation by atoms in a two-color laser field. _Phys. Rev. A_ 81 no. 6, p. 063407 (2010).
  * (55) R. Wong.  _Asymptotic approximations of integrals_ (SIAM, Philadelphia, 2001), pp. 366–372.
  * (56) M. V. Frolov, N. L. Manakov, T. S. Sarantseva and A. F. Starace.  High-order-harmonic-generation spectroscopy with an elliptically polarized laser field. _Phys. Rev. A_ 86 no. 6, p. 063406 (2012).  UNL eprint.
  * (57) T. Sarantseva, M. Frolov, N. Manakov, M. Y. Ivanov and A. F. Starace.  Harmonic generation spectroscopy with a two-colour laser field having orthogonal linear polarizations. _J. Phys. B: At. Mol. Opt. Phys._ 46 no. 23, p. 231001 (2013).  UNL eprint.
  * (58) Y. Okajima, O. I. Tolstikhin and T. Morishita.  Adiabatic theory of high-order harmonic generation: One-dimensional zero-range-potential model. _Phys. Rev. A_ 85 no. 6, p. 063406 (2012).
  * (59) M. A. Khokhlova and V. V. Strelkov.  Phase properties of the cutoff high-order harmonics. _Phys. Rev. A_ 93 no. 4, p. 043416 (2016).  arXiv:1511.02133.
  * (60) H. Schomerus and M. Sieber.  Bifurcations of periodic orbits and uniform approximations. _J. Phys. A: Math. Gen._ 30 no. 13, p. 4537 (1997).  arXiv:chao-dyn/9701022.
  * (61) S. Popruzhenko.  Keldysh theory of strong field ionization: history, applications, difficulties and perspectives. _J. Phys. B: At. Mol. Opt. Phys._ 47 no. 20, p. 204001 (2014).
  * (62) L. E. Chipperfield, J. S. Robinson, J. W. G. Tisch and J. P. Marangos.  Ideal waveform to generate the maximum possible electron recollision energy for any given oscillation period. _Phys. Rev. Lett._ 102 no. 6, p. 063003 (2009).
  * (63) L. E. Chipperfield, J. W. Tisch and J. P. Marangos.  Derivation of the perfect wave for optimising the strong-field driven electron recollision energy. _J. Mod. Opt._ 57 no. 11, pp. 992–998 (2010).
  * (64) S. Haessler, T. Balčiunas, G. Fan et al.  Optimization of quantum trajectories driven by strong-field waveforms. _Phys. Rev. X_ 4 no. 2, p. 021028 (2014).
  * (65) 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 no. 24, pp. 3535–3538 (1992).  zenodo:1233893.
  * (66) D. B. Milošević.  Cut-off law for high-harmonic generation by an elliptically polarized laser field. _J. Phys. B: At. Mol. Opt. Phys._ 33 no. 13, p. 2479 (2000).
  * (67) A. Flegel, M. Frolov, N. Manakov and A. F. Starace.  Cutoffs of high-energy plateaux for atomic processes in an intense elliptically polarized laser field. _J. Phys. B: At. Mol. Opt. Phys._ 38 no. 1, p. L27 (2004).  UNL eprint.
  * (68) C. Figueira de Morisson Faria, M. Dörr, W. Becker and W. Sandner.  Time-frequency analysis of two-color high-harmonic generation. _Phys. Rev. A_ 60 no. 2, pp. 1377–1384 (1999).
  * (69) V. Averbukh, O. E. Alon and N. Moiseyev.  High-order harmonic generation by molecules of discrete rotational symmetry interacting with circularly polarized laser field. _Phys. Rev. A_ 64 no. 3, p. 033411 (2001).
  * (70) H. K. Avetissian, B. R. Avchyan and G. F. Mkrtchian.  Generation of harmonics via multiphoton resonant excitation of hydrogenlike ions in an x-ray free-electron-laser field. _Phys. Rev. A_ 90 no. 5, p. 053812 (2014).
  * (71) O. Neufeld, A. Fleischer and O. Cohen.  High-order harmonic generation of pulses with multiple timescales: selection rules, carrier envelope phase and cutoff energy. _Mol. Phys._ 117 no. 15-16, pp. 1956–1963 (2019).
  * (72) O. Neufeld, D. Podolsky and O. Cohen.  Floquet group theory and its application to selection rules in harmonic generation. _Nat. Commun._ 10 no. 1, pp. 1–9 (2019).
  * (73) A. Jašarević, E. Hasović, R. Kopold, W. Becker and D. B. Milošević.  Application of the saddle-point method to strong-laser-field ionization. _J. Phys. A: Math. Gen._ In press (2020).
  * (74) S. Odžak and D. B. Milošević.  High-order harmonic generation in the presence of a static electric field. _Phys. Rev. A_ 72 no. 3, p. 033407 (2005).
  * (75) L. E. Chipperfield, L. N. Gaier, P. L. Knight, J. P. Marangos and J. W. G. Tisch.  Conditions for the reliable production of attosecond pulses using ultra-short laser-generated high harmonics. _J. Mod. Opt._ 52 no. 2–3, pp. 243–260 (2005).
  * (76) D. B. Milošević.  Strong-field approximation and quantum orbits.  In D. Bauer (ed.), _Computational Strong-Field Quantum Dynamics: Intense Light–Matter Interactions_ , pp. 203–225 (De Gruyter, Berlin, 2017).
  * (77) E. Hasović, W. Becker and D. B. Milošević.  Electron rescattering in a bicircular laser field. _Opt. Express_ 24 no. 6, pp. 6413–6424 (2016).
  * (78) D. B. Milošević and W. Becker.  X-ray harmonic generation by orthogonally polarized two-color fields: Spectral shape and polarization. _Phys. Rev. A_ 100 no. 3, p. 031401 (2019).
  * (79) D. B. Milošević.  Quantum-orbit analysis of high-order harmonic generation by bicircular field. _J. Mod. Opt._ 66 no. 1, pp. 47–58 (2019).
  * (80) D. B. Milošević and W. Becker.  Improved strong-field approximation and quantum-orbit theory: Application to ionization by a bicircular laser field. _Phys. Rev. A_ 93 no. 6, p. 063418 (2016).
  * (81) S. Mitra, S. Biswas, J. Schötz et al.  Peak suppression in two-colour high harmonic generation. _J. Phys. B: At. Mol. Opt. Phys._ 53 no. 13, p. 134004 (2020).
  * (82) C. Figueira de Morisson Faria and A. S. Maxwell.  It is all about phases: ultrafast holographic photoelectron imaging. _Rep. Progr. Phys._ 83, p. 034401 (2019).  arXiv:1906.11781.