The imaginary part of the high-harmonic cutoff

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.

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.
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 [1,2], 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 [3,4].
This emission can be modelled using a wide range of approaches, from classical heuristics [3,4] to intensive numerical computations [5], but the quantitative models that most closely follow the overall intuition are quasiclassical methods [6,7], 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 [8][9][10], i.e., quasiclassical trajectories whose start and end times are complex [11,12].
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 steepestdescent method known as the saddle-point approximation (SPA), to discrete contributions coming from the saddle points of the integrand [13][14][15]. These saddle points represent the quasiclassical trajectories, and they typically come in pairs -most notably the 'short' and 'long' trajectories [16,17]. Each pair of trajectories approaches each other over the harmonic plateau and then performs a Stokes transition at the harmonic cutoff [18,19], 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 [20][21][22], particularly since the saddles tend to move quickly, and approach each other very closely, at the harmonic cutoff. * emilio.pisanty@icfo.eu In this work we construct a quantum-orbit classification scheme based on a natural notion of complex-valued harmonic-cutoff times t 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. 1(b) can be cleanly organized programmatically into families of trajectories, as shown in Fig. 1(c), 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 iden- Solving the HHG saddle-point equations returns a discrete set of saddle points as a function of the harmonic photon energy Ω, 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 hc (triangles) and the separatrices through them, which allow the cloud of saddle-point solutions in (b) to be organized into individual trajectories. We model neon in a field of wavelength 800 nm and intensity 2 × 10 14 W/cm 2 . tified, 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 Ω max ∼ I p + 3.17U p with the ionization potential I p and the ponderomotive energy U p of the field (previously derived both classically [3,4] and via a systematic expansion in I p /U p [6]). However, our method extends trivially to higher harmonic content as well as to tailored polarizations.
Moreover, the information at t 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 [23].
This understanding of the high-harmonic cutoff Ω 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 quantumpath interference (QPI) [17] between the pair of orbits.
In particular, the zeros of the imaginary part of the harmonic cutoff energy Im(Ω 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, Im(Ω 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 counterrotating circularly-polarized 'bicircular' fields [24][25][26][27][28][29][30][31][32][33][34] 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 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 [35], the harmonic cutoff is a spectral caustic, in the 'fold' class of diffraction catastrophes [36]. In this sense, the complex harmonic-cutoff energy Ω 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 [36][37][38][39][40][41][42][43][44], and our approach provides a useful tool for exploring and describing the higher-dimensional bifurcation sets at the heart of the higherorder 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 A and constructing a model action with a cubic action in Section I B, which we explore in depth in Sections I C and I D, where we construct the classification algorithm. 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 Ω 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 [45], and our specific implementation is available from Ref. 46.

A. 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 |g or in a laser-driven continuum where the effect of the atomic potential is negligible. The calculation [6,7,11,12] then gives the harmonic dipole in the form × e −iS V (p,t,t )+iΩt + c.c., (1) where the three integrals over the times of ionization and recollision, t and t, and the intermediate momentum, p, form a reduced Feynman path integral summing over a restricted set of relevant orbits [8]. This Feynman sum is modulated by ionization and recollision amplitudes given by the transition dipole matrix element d(k) = k|r|g and the scalar function Υ(k) = (I p + 1 2 k 2 ) k|g , evaluated at the ionization and recollision velocities: A(t) is the vector potential of the laser and p is the canonical momentum, so v(t) = p + A(t) is the recollision velocity and m e v(t) = v(t) is the kinematic momentum. (We use atomic units with m e = = 1 throughout, and we consider neon in a linearly-polarized field F(t) = F 0êz cos(ωt) of wavelength 800 nm and intensity 2 × 10 14 W/cm 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, 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 ω, the action scales with the ratio z = U p /ω of the field's ponderomotive energy U p = F 2 0 /4ω 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 [13][14][15], also known as the saddlepoint 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 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 α(t) = t t v(τ )dτ to the ion at the time of recollision. 1 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 p guarantees the unique solution for any arbitrary pair of ionization and recollision times (t, t ). Once the momentum integral has been performed in this way, the expression for the harmonic dipole takes the form where the fractional power of the excursion time τ = t − t represents the dispersion of the released wavepacket in position space. 2 We notate S(t, t ) = S(p s (t, t ), t, t ) 1 We 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. 2 If the time integrals are performed numerically, this factor needs to be regularized at τ → 0, to account for a breakdown of the saddle-point approximation in that limit [47]. If the time integrals are also evaluated via saddle-point methods, however, this is not required, as that limit is not used. 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, and the saddle-point equations for the amended action read and they can only be solved numerically. 3 A typical set of solutions for these saddle-point equations is shown in Fig. 1(b): the solutions form a discrete set of points, which shift when the harmonic frequency Ω changes. These discrete points thus trace out individual curves on the complex t and t planes, which form the individual quantum orbits.
Once the saddle points have been found, the SPA gives the harmonic dipole as 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. 2(c) 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 [16,17].
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 steepestdescent integration contour, a property that can be nontrivial to determine. Fig. 2(b) 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 [18][19][20]. At the anti-Stokes transition, which is defined as the point where Im(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 Re(S) for the two orbits is equal and then changes order. This is an important change, as the steepestdescent method requires integration contours to follow the contour lines of Re(S): thus, the change in the ordering of Re(S(t s , t s )) 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 Ω, 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 D(Ω) 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 [18,19] (UA), which is shown in Fig. 2(c).
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. 4 In practice, this is a high-friction point in the calculation, requiring expensively tight grid spacings in Ω to accurately resolve the avoided crossings, or the additional design of an adaptive grid to increase the energy resolution there [20][21][22]. Moreover, this problem is compounded in more complex polychromatic fields, where the saddlepoint 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. 2(b) along the diagonal into the two curves shown. B. 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. 2(b) 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 (For simplicity, we restrict our attention for now to actions with only one variable, t, which corresponds to solving (6b) first for t s = t s (t), giving S(t) = S(t, t s (t)) and then examining (6a). We shall lift this restriction later.) Moreover, we use the clear symmetry evident in Fig. 2(b), 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 which can be integrated to find 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 and then compare the linear term, e +iAt 2 s t , with the Fourier kernel e +iΩ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 2 s = Ω − (Ω c + iη) into a variable and a constant part, and thus define Here the constant part Ω hc = Ω c +iη has a nonzero imaginary part +iη, as we do not have any guarantees that At 2 s is real, but its real part Ω c can be set to zero if desired, since it simply acts as an offset for the variable Ω. Here the functional dependence on Ω 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 (Ω) that result from this identification form a good model for the avoided-crossing behaviour in Fig. 2 In addition, if we now separate S(t) = S V (t)−Ωt as we did in (2) above, the saddle-point equation (9) can now be rephrased as the requirement that with Ω 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 of the derivative of the model Volkov action S V (and also of the full action S, since they differ by a real number). This insight then lights the way further: our principal task is to look for objects between the two curves in Fig. 2(b) that will help us separate them, and the contour map of Im dS V dt 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. 2(b) 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 hc for this model. Like all saddle points, this can be found as a zero of the Contour map of Im ∂S V ∂t (t, t s (t)) 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 is real-valued. The harmonic-cutoff points of Fig. 1(c) (triangles) are, correspondingly, the saddle points of this derivative; the coloured lines show the contours at these saddle points.
derivative of the contour map's original function, or in other words, via In the particular case of our model action (14a), for which 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 Im dS V dt 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 [48,49], in the context of tunnel ionization.

C. 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: with the second derivative vanishing by construction. 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 2 s = Ω−(Ω c +iη) as set above. In particular, this term directly encodes the saddle-point solutions, and we can recover these explicitly by inverting that relationship: Here the square in At 2 s 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 ± 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 = t s (t), then (6a) can be expressed in the form where it explicitly involves the inverse of the many-to-one, complex-valued function S V (t, t s (t)), which is precisely the type of problem encoded by Riemann surfaces. 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 We exemplify this directly in Fig. 4, showing the contour map of Im ∂S V ∂t (t, t s (t)) 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. 1(c). Similarly, the harmoniccutoff times t hc used there are found as the saddle points of the ∂S V ∂t (t, t s (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 Ω, but the same is true for the contours of Im dS V dt , 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 dS V dt (t hc ) is simply the frequency offset Ω. This offset is carried from the saddle point t hc to the quantum-orbit lines by means of the contour lines of Re dS V dt . 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 as the frequency of the harmonic cutoff.
Having made this leap, we must now confront the fact that dS V dt (t hc ) as given in (16a) also has a nonzero imaginary part, which similarly needs to be addressed. This imaginary part has already played a part, and it is explicitly shown in Fig. 3 as the height of the dS V dt saddle with respect to the zero contour where the quantum orbits lie. As such, the imaginary part η 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 A below, this then controls the strength of quantum-path interference between them.
Moreover, the sign of η 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 hc . Generally, this imaginary part of dS V dt is fixed by the problem and it can be either positive or negative, but, as we shall see in Section IV, sign changes in η 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 Ω to acquire complex values, i.e., by considering the analytical continuation of D(Ω). Doing this amounts to directly affecting the value of η in (14a), and if Ω is taken at the complex harmonic cutoff itself, Ω c + iη, then the constant term in the saddle-point equation vanishes, leaving a double zero of the form at t = t hc . In other words, the harmonic-cutoff times t hc are the locations where the saddle points for the two quantum orbits fully coalesce into a single point.

D. 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 dS V dt and, as such, it controls the orientation of the saddle in dS V dt at t hc , as shown in Fig. 3. (Indeed, we chose a nonzero phase for A = e 10 • 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 harmoniccutoff times t hc to classify the saddle points into quantum orbits: in order to turn the clouds of points in Fig. 1(b) into the ordered curves of Fig. 1(c), 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 hc to the closest point on the Im dS V dt = 0 contour that the quantum orbits follow. As we noted above, this point occurs at Ω = Ω c , which means that we simply need to solve the saddle-point equation (14a) at that frequency, i.e., or, in terms of the explicit derivatives of the action and with an explicit relation to the saddle center t hc , .
We can then write the explicit condition for the separatrix by treating δt sep and t s − t 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 hc is the sign comparison with t hc defined as in (15) and δt sep defined as in (23). In a practical calculation, there will typically be multiple harmonic-cutoff times, alternating between nearthreshold 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. 1(c) -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.

E. 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 = t s (t) and with that the singlevariable S V (t) = S V (t, t s (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 and here the partial derivative in the second term, , vanishes by the definition of t s (t). However, if we now turn to the second derivative, the procedure no longer works, and the equivalent calculation, returns a term that includes dt s dt (t). This cannot be evaluated explicitly within this system, as t s (t) itself is only defined implicitly.
To resolve this, we use the implicit definition of t s (t), and then differentiate with respect to t, to obtain which can then be solved for the implicit derivative as Finally, this can be substituted into (26) to give the form Here we have dropped the right-hand side evaluations at (t, t s (t)) for simplicity, but it is also possible to understand the right-hand side as a function of (t, t ) 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, using the substitution rule (29) at each step to retain explicit formulations.

F. 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. 1(c), this is not the case for the linearly-polarized field we use above, since half of the harmonic-cutoff saddles t 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 Ω c + iη = I p , and the threshold harmonics at Ω = I p always involve an exact saddle coalescence. In other words, for threshold harmonics, action landscape's saddles are not gaussian, but exact monkey saddles as shown in Fig. 6(a) below; this was noticed as early as Ref. 6, but has not drawn much attention since then. 5 To understand this behaviour, we return to the saddlepoint equations (6), which for linearly-polarized fields can be rephrased in the simpler, scalar form Here the first equation is the crucial one, since at Ω = I p it tells us that i.e., the recollision velocity v r (t, t ) = p s (t, t ) + 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, (30), whose zeros determine the harmonic-cutoff times. The partial derivatives involved in (30) can be evaluated by differentiating (31a) with respect to t and t , which yields Both of these derivatives share a common factor of v r (t, t ) = p s (t, t ) + A(t), and this then means that the second derivative d 2 S V dt 2 (t) will always vanish when evaluated at solutions of the first-order saddle-point equations (31) at Ω = 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 dS V dt (t hc ) also vanishes, so the direction vector δt sep from (23) also vanishes. In practice, this can cause numerical instability, with Im dS V dt (t hc ) 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 Im dS V dt (t hc ) to be used there, as either direction will work well for the separatrix. (For the plots in Fig. 4 and Fig. 1(c), we used a small ellipticity of = 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 Ω = 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 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 hc itself.
A. 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), over an integration contour C which should start at the ionization time t (or some suitably compatible valley of Im(S) to the left of t hc ) and end at t → ∞. This integral is essentially in Airy form, and it is almost structurally identical to the Airy function's integral representation [50, Eq. (9.5.E4)], This is expected, since the harmonic cutoff is a 'fold' catastrophe, whose associated diffraction-catastrophe integral is the Airy function [51].
To bring our representation (34) into the canonical form (35), the core transformation is to eliminate the coefficient in front of the cubic term, In a sense, this is relatively simple, as it boils down to a change in integration variable tot = 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, for k ∈ {0, 1, 2}, which requires dedicated attention. The variable change itself is essentially trivial, and (34) transforms under (37) to This is essentially in explicit Airy form, so long as the contour is correct, and thus we can substitute in the Airy function as where k needs to be chosen such that the altered con-tourC = ie 2πik/3 A 1/3 C is compatible with the standard contour in (35), which we depict in Fig. 5.
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. 6, 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π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. 6, A can lie very close to commonly-used choices for the branch cut of the radical.) 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 − π 2 < arg(t) < − π 6 and end at π 6 < arg(t) < π 2 , i.e., that it go down the valleys that surround e −iπ/3 and e iπ/3 .
In practice, if there are known constraints on the behaviour of the integral (in particular, exponential decay at Ω > Ω 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 Ω > Ω c requires e 2πik/3 A 1/3 to have a negative real part.
Our result for the model integral, (39), now allows us to have a closer look at the role of the imaginary part of the complex harmonic-cutoff energy Ω 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. 7. The HCA approximates the integral (34) as an explicit Airy function in Ω (red line), which transitions from oscillatory to exponential-decay behaviour at Ω c = Re(Ω 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 η = Im(Ω hc ) is nonzero. To increase the interference features, thus, we can artificially reduce (blue line) or eliminate (green line) this imaginary part; conversely, increasing η (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 [11,52,53], 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 [18,19,54] (as well as the equivalent constructions in Refs. [55][56][57], 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  Fig. 4, taking a complex Ω equal to Ωc + iη at the various harmonic-cutoff times t 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 byC = ie 2πik/3 A 1/3 C, choosing the branch index k so that the integration contour matches the canonical Airy contours of 5. The phase arg(A) is inset in (b) and (c), and it equals −87.97 • , −0.011 • and −0.155 • , resp., for the three monkey saddles shown in (a), consistently with their different orientations.
suited for efficient evaluation and clear analysis of the harmonic strength (as well as phase properties [58]) at the cutoff.

B. The full HHG integral
We now return to the full two-dimensional integral for HHG, (5), and transform it into the model form (34) so that the HCA can be used to estimate it. For simplicity, we write the integral in explicit prefactor-action form, To transform this into the model form (34), we need to reduce the integral to a single dimension, and then translate the origin to the relevant harmonic-cutoff time t hc . 6  (39), for Ω hc and A corresponding to the first-return cutoff in Fig. 1. The coloured lines have the imaginary part η = Im(Ω 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 Im( ∂S V ∂t ) contour to follow in 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πik/3 A 1/3 with −|A| 1/3 . The gray vertical line is at Ω = Ωc.
The first is achieved, as mentioned earlier, by doing a saddle-point approximation on the t integral only, and this returns the harmonic dipole in the form with t s (t) defined implicitly via (6b). After this, we perform a Taylor expansion of the Volkov action S V (t, t s (t)) at the harmonic-cutoff time t hc up to third order, where A and Ω hc are the third and first derivatives of S V (t, t s (t)), as set in (16) and evaluated as constructed in Section I E; here we have also allowed the integration contour C to vary so that it can pass through t hc as apthe Uniform Approximation [18] (as well as its earlier analogues in a semiclassical context [59]), the coordinate separation done here is in essence an application of the splitting lemma of catastrophe theory [35], 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 B 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.
propriate. Finally, we move the integration origin to t hc , keeping the explicit t =t + t hc for notational simplicity.
To apply the Harmonic-Cutoff Approximation, we now assume that the prefactor varies slowly enough at t 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. 6), 7 which gives The integral is now in the form of (34), with an added functional dependence on Ω coming from the term e iΩt hc . This term results from the translation of the time origin, and it can have a sizeable effect on the harmonic yield if Im(t hc ) is nonzero. That said, we can now use the result (39) for the model form directly. We obtain thus our final result for the HCA, where the harmonic-cutoff time pairs (t hc , t hc ) are found by solving the simultaneous equations and where the coefficients Ω hc and A hc for each pair of harmonic-cutoff times (t hc , t hc ) are given, as in (16), by with the third derivative d 3 dt 3 understood in the constrained sense of Section I E. 8 Our implementation of this approximation, in the Wolfram Language, is available in the RBSFA software package [45].
We show the results of this approximation, compared to the SPA and the UA (as described in Refs. 18,19), in Fig. 8, plotted as in Fig. 2(c). When the prefactor is fully ignored, as in Fig. 8(a), the HCA is extremely accurate at the harmonic cutoff, and it retains a good qualitative accuracy as Ω descends from Ω 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 τ −3/2 is included, as shown in Fig. 8(b). This happens because the τ −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 D(Ω), 9 which we plot in Fig. 8(c) 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 a set of saddle-point equations (46), 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 saddlepoint equations -the system in (6) -for each harmonic frequency Ω 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) [5], 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 highdimensional 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 [61][62][63]. 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 Ω 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', which relates the cutoff frequency Ω max to the target's ionization potential I p and the ponderomotive energy U p = F 2 /4ω 2 of the field. This relationship was first uncovered in numerical simulations [64] and subsequently explained using the classical three-step model [3,4]: 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 [6] re-derives this cutoff law, giving also a quantum correction to the I p term of the form To obtain this formulation, the cutoff frequency is defined as the maximum of the recollision kinetic energy, Re(v(t r ) 2 ), taken under the restriction that its imaginary part Im(v(t r ) 2 ) vanish. The dynamics are then analyzed using a systematic expansion on I p /U p , using the purelyclassical case at I p = 0 for the leading 3.17U p term, after which the cutoff law can be derived in the exact form 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 linearlypolarized monochromatic fields, and it can be extended to elliptical polarizations [65,66], but it is laborious to apply for polychromatic combinations [67] (where the complex waveform makes analytical calculations challenging) and for tailored field polarizations [27], 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 Im(v(t r ) 2 ) = 0 condition to Im(t r ) = 0 [28], focusing only on the classical recollision velocity [61,62], examining the change in the intensity scaling of the harmonic yield at a fixed harmonic order [16], the use of graphical methods based on tangents to the optical waveform [67], and extracting the cutoff from the Bessel-function expansion of the oscillatory integral [68,69], as well as direct numerical methods used to extract the cutoff from the results of TDSE simulations [70]. 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.
As we argued earlier, the real part Ω c of our complexvalued Ω hc forms a natural candidate as a precise definition of the harmonic-cutoff frequency Ω 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 [71,72]. 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. 9(a): the variation of Re(Ω hc ) with U p is linear with an I p -dependent offset, which closely follows the quasiclassical quantum-corrected cutoff law (49).
However, our formalism allows us to go beyond this level, since Ω hc also has an imaginary part, whose scaling with U p and I p is shown in Fig. 9(b). As rough behaviour, Im(Ω hc ) decreases with U p and increases with I p , and simple testing shows that the leading asymptotic component is the behaviour At first glance, this shows very different behaviour to the real part, which follows (50). However, once a factor of I p has been set apart, as we did for the real part, the scaling for Im(Ω hc ) is simply the Keldysh parameter, γ = I p /2U p . This is therefore indicative that the parameter in the quantum scaling law (50) should be amended from I p /U p to its square root, so it should read In this form, (52) essentially acts as a definition for F (γ) = (Ω 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 γ. This is indeed the case, as we show in Figs. 9(c) and 9(d) for the real and imaginary parts, respectively: plotting F (γ) by changing I p reveals the same universal curve for a range of different values of U p . Moreover, we can extract the low-γ behaviour here to obtain with the first two terms shown as the gray dashed line in the figure. This understanding coincides with the existing theory (so, in particular, Fig. 9(c) coincides with Fig. 5 of Ref. 6 when plotted as a function of γ 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γ for F (γ), as defined here for Ω hc , is still lacking. We summarize our parameter-scaling results, and their relationship to the existing theories, in Table I.

Theory Cutoff law
Corkum [3] Ωmax ≈ 3.17Up + Ip Kulander [4] Lewenstein [6] Ωmax = 3.17Up + F (Ip/Up)Ip This work 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 D. 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 and thus that it has a sensitive dependence on the sign of 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. 9) then the sign of the imaginary part is unlikely to change (as was seen to be the case in Fig. 9(b)).
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 η will change, and the direction of the separatrix, δt sep , will switch to one 45 • 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 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 [24][25][26][27][28][29][30][31][32][33][34], formed by the combination of two counter-rotating circularly-polarized strong laser pulses at frequencies ω and 2ω, These fields have been the focus of widespread interest in recent years because they are subject to a spin selection rule [24,32] which ensures that the harmonics of the driver are circularly polarized [25]. 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 [28][29][30][31]. If we keep the total intensity constant, at I tot = 2 × 10 14 W/cm 2 for concreteness, then the shift in intensity is best described by the mixing angle θ, which we use to define the individual field amplitudes as where F = √ 2 × 0.053 a.u. is the peak field strength. We show in Fig. 10 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. 10(a), the short orbit connects up to the positive-imaginary saddle, while after the transition, in Fig. 10(c), it connects to negative imaginary time. Between these two there must always lie a transition, which we show in Fig. 10(b): here the saddles fully coalesce, producing a monkey-saddle landscape like the ones shown in Fig. 6, with a purely-real HCA Airy function corresponding to enhanced quantum path interference between the two orbits.
More importantly, however, this transition separates the instances with η > 0 from those with η < 0, so at the transition itself η must vanish -and, as with the exact coalescences we studied in Section I F, no (nonzero) unique separatrix direction δt 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 η = ±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 η will work equally well. In Fig. 10(b), we mark the transition by showing both possible separatrices in Topological transition and reconnection in the quantum orbits of a bicircular field as the mixing angle θ changes. Here we show a detail of the quantum-orbit layout for the short and long trajectories, plotted as in Fig. 1(c), with the triangles marking the harmonic-cutoff times t hc and the coloured separatrices marking the regions occupied by the different orbits, as described in Section I D. In (b) we show the topological transition, where η = Im( dS V dt ) vanishes and the separatrix becomes undefined, which we show by highlighting the t hc with a black border and plotting both options (using an artificial η = ±1) as the gray cross. In this case, both choices of separatrix are acceptable for saddle classification, since the saddles fully coalesce.
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 hc with a black border for further clarity.
On a more practical footing, the fact that η = Im( dS V dt ) 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 η as a function of θ.
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. 11. 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. 10(b) have 'confused' identities between the short and long orbits, the short orbit 'mingles' at low Ω with the short-time orbit shown in cyan, at two separate transitions (shown in Fig. 11(c) and Fig. 11(b)), the long orbit connects with the second return (with the transition shown in Fig. 11(g)), and so on.
In other words, the quantum orbits here form a single unified topology, and they should be understood as such. We show this topology in Fig. 12, by unfolding the θ 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 mentioned in Section I C is the Riemann surface formed by the saddle points: these are defined, in essence, as the discrete solutions of an equation of the form which means that they are the multivalued inverse of an analytical function evaluated over the real axis on Ω. This is precisely the definition of a Riemann surface, with the various quantum orbits forming different branches which we perceive as separate objects. (A more explicit embodiment of this branch choice is apparent in (17), where the simplified model makes the branch structure easier to visualize.) 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 Ω over the real axis, and those inverse images will normally be well separated. In Fig. 12, however, as well as in any other situations where the field shape changes over a control parameter, we see how the change in the field shape alters the Riemann surface, bringing these once-separate images into collision, and reconnection, with each other, so that this structure can no longer be ignored.
(In any case, just to state this clearly, the surface shown in Fig. 12 is not the Riemann surface of the quantum orbits. Each individual mixing angle has a Riemann surface, which forms a manifold of (complex) dimension one, embedded in the four-dimensional space of complex (t, Ω) pairs. The precise topology of the Riemann surfaces for quantum orbits in HHG is at present an open question.) Now, if we are identifying the different quantum orbits of the problem as different branches of a single Riemann surface, then it behooves us to identify the branch points where the arbitrary cuts between these different branches must be made. However, we have already identified these branch points: they are precisely the harmonic-cutoff times t hc we began with, i.e., the double zeros at complex Ω where the two different branches meet and merge.
In more general terms, the reconnection phenomena 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 αβm, have been used widely in the literature, both for monochromatic fields [19,20,[73][74][75] as well as more complex waveforms [76][77][78][79]. 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 hc as well as the parameters where their branch cuts cross the real Ω 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 FIG. 11. Topological transitions and reconnections in the quantum orbits of a bicircular field as a function of the mixing angle θ, plotted as in Fig. 10. As θ 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 Fig. 12. At each panel the inset at bottom right shows the shape of the electric field at the relevant mixing angle.
FIG. 12. Unified surface formed by the quantum orbits in Fig. 11, when the θ dependence is unfolded as a third dimension. (Thus, a cut through the surface at constant θ will return a panel from Fig. 11.) 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.
FIG. 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. 11. For bicircular fields, attention is generally focused on the short quantum orbit (blue curve), since at equal intensities (θ = 45 • ) it dominates the spectrum (as shown in (b), which closely follows Fig. 8 of Ref. 28). However, this is no longer the case at mixing angles with a high intensity for the 2ω 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.
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ω 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 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 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 harmoniccutoff times can be used to provide a natural identification for the energy of the cutoff, Ω hc = ∂S V ∂t (t 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 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 Im( ∂S V ∂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 hc . Likewise, the 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 return = ∂S ∂t is a function of the return time, and if we want to find its extrema, then we simply need to solve the equation ∂ 2 S ∂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 ∂ 2 S ∂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 postcutoff 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 [36][37][38], but recent years have seen a marked increase in interest in that perspective on the dynamics [39][40][41][42][43][44], as laser sources achieve better control over polychromatic combinations with high relative intensities (which is slated to continue to increase [80]). This has opened the door to the observation of higher-order catastrophes in strong-field observables, involving higherdimensional 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 harmoniccutoff times we have discussed here should have analogous structures in ionization experiments (most notably high-order above-threshold ionization [81]), 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 higherdimensional 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.