Modelling intensity volume averaging in ab initio calculations of high harmonic generation

We present an approach to assess the survival of single-atom effects in the macroscopic high-harmonic generation (HHG) spectrum, by accounting for focal volume averaging. We apply this technique to R-matrix with time-dependence (RMT) studies, which are designed to include the full multielectron response of an atom. Such an approach allows the assessment of which features of an experimentally-measured HHG spectrum of diffuse gases may be traced directly to single-atom effects and vice-versa. While accounting for the phase of harmonics produced at different locations in the focal volume gives the most accurate results, a simplified approach, using a smaller number of RMT calculations, is found to provide comparable conclusions. We apply these approaches to compute intensity-averaged harmonic spectra in two different experimental regimes.


Introduction
Attosecond physics studies the dynamics induced by lightmatter interaction on ultrashort timescales. Central to the field of attosecond physics is the process of high harmonic generation (HHG) which has been used extensively both as the means of generating ultrafast laser pulses [1][2][3][4] and as a spectroscopic tool to probe the electronic structure and dynamics of atoms and molecules [5,6]. The process of HHG can be understood using the semi-classical three step model [7,8]. In this model, the laser field (i) causes tunnel ionisation of the electron, (ii) launching the electron on a classical trajectory driven * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. by the laser field before (iii) the electron returns to the parent ion, emitting a high-energy photon upon recombination. All information about the electron's interaction with the parent ion is encoded in the emitted photon: multielectron interference, nuclear dynamics and other spatial or temporal effects can be extracted by carefully analysing the spectrum [9,10]. The three step model predicts the cutoff energy for the harmonics, given by where I p is the ionisation energy and U p is the ponderomotive energy (the cycle-averaged electron kinetic energy in the laser field). The ponderomotive energy is U p ∝ I/4ω 2 , where I and ω are the laser field intensity and frequency, respectively. HHG experiments are inherently macroscopic, involving many atoms responding coherently when illuminated by a focused laser beam. We may exploit this macroscopic response to measure HHG as the coherent buildup of many single atom-emissions [11]. However, if the goal of the HHG experiment is to understand the electronic response, very often we want to decouple the single-atom behaviour from the macroscopic effects. The two main macroscopic effects are propagation through the medium, and the intensity volume effect. Both are caused by the relative geometries of the laser focus and the target medium, which cause each atom in an ensemble to experience a different electric field. Understanding the macroscopic build-up of individual harmonics, including their phases and their interplay with the wavefront of the generating laser pulse, is key to increasing both the HHG conversion efficiency [12] and the intensity of attosecond pulses generated through HHG [4]. A review of the issues involved in optimising HHG in gases is presented in [12]. However, considerations beyond the focal volume effect are outside the scope of the present paper.
Experiments can benefit from the support of ab initio calculations which usually model a single atom irradiated by a laser at a single intensity, providing an insight into the microscopic detail [13][14][15][16]. It is common to attribute features in experimentally-observed HHG spectra to single-atom effects [10,17,18], and often it is possible to verify or support experimental results with ab initio calculations [13,15,19]. As experimental techniques mature, the increased resolution they provide will require ever-more detailed calculations for prediction and interpretation of results. However, the predictive power of single-atom, ab-initio calculations is limited by the omission of macroscopic effects, while the interpretive power of experiment can be limited by their inclusion. In order to bridge this gap between simulation and experiment, we require some means of assessing which features in a single-atom spectrum may 'survive' experimental conditions to be measured and, on the other hand, which measured features may be traced to effects at the single-atom level.
The most straightforward way to replicate the intensity volume effect with theoretical simulations is to run calculations at thousands of intensities and thereby compute the gross response. Whilst this method is suitable for simple model simulations which are computationally cheap, it is not feasible for detailed, computationally demanding ab initio methods. A more appropriate technique for averaging observables obtained from ab initio calculations is to use a small sample of intensities, weight the contribution from each intensity in the sample and then sum them to obtain an intensity-averaged result. This approach has been demonstrated previously for calculating an incoherent intensity average of observables such as photoelectron spectra [20] and momentum distributions [21]. Because HHG relies on the coherent response of each atom in an ensemble, any technique to calculate an intensity average must account for the phase of the harmonic response. Moreover, the harmonic phase depends on the experimental geometry [4], so this must also be accounted for in the averaged result.
Many papers have been written on various intensity averaging and/or phase-matching methods for ab initio calculations of HHG in atoms, molecules and solids [22][23][24][25][26] which focus on the finer details of macroscopic harmonic spectra. In this paper, however, we are not concerned with accurately reproducing the individual features of harmonic peaks or other fine details of experimental HHG spectra. Rather, we seek a method which can give insight into which qualitative aspects of the microscopic HHG spectra are observed in experiment and which aspects may be obscured due to the macroscopic response.
Here, we outline two techniques for approximating the intensity volume effect on high harmonic spectra obtained from ab initio calculations: a 'full' model which considers both the geometry of the laser focus and the phase of the harmonic response and a 'simple' model which neglects the relative phase of harmonics arising from different points in the laser focus. These methods can then be used to identify which features of the microscopic response can be disentangled from the macroscopic spectrum. We use the ab initio R-matrix with time (RMT) method [27] to test and demonstrate this approach to intensity averaging. Previously, calculating an intensity average of RMT results was not feasible due to the computational cost. However, this is now possible due to recent gains in efficiency of the RMT method and the low computational cost of the proposed intensity averaging technique. To validate the method, we look at two different experimental schemes with different atomic targets: xenon irradiated by a 1800 nm laser pulse and argon irradiated by a 780 nm pulse.

Intensity averaging
Experimentally measured, 'far-field', harmonics are not detected at the exit face of the gas medium, but rather at some distance, d, away on a detector. Because the harmonics generated at different locations in the laser focus interfere with each other, we must account for their phase. Also, to account fully for the geometry of the experiment, the radius, R, of the detector at which the far-field harmonics are measured must also be considered.
In order to model this experimental geometry we first consider the gas jet containing the target atoms to be sufficiently tenuous to neglect the propagation effects. This is a reasonable assumption, as most HHG studies designed to measure singleatom effects implement this experimental condition as far as possible. This then allows us to treat the laser beam focus as 2-dimensional. Having made these assumptions, we may then adopt a cylindrical coordinate system (r, θ, z), wherein harmonic radiation is generated at the point (r, ϕ, 0) and is measured at point (R, 0, d). The Fresnel diffraction formula [28] relates the spectrum of far-field harmonics, E f (ω, R)| z=d , to the near-field harmonics, E ′ (ω, r)| z=0 , generated at a distance r from the centre of the laser focus: Here, we use the Fresnel diffraction formula in a cylindrical coordinate system and assume that the transverse wave of the driving laser beam is perfectly circular, so that there is cylindrical symmetry and E f (ω, R) does not depend on ϕ. Thus, rather than integrating over the entire focal volume, the Fresnel diffraction formula, equation (2), becomes an integral over r only. The distance travelled, D, between the point of emission and the point of detection is: Now, we may assume that d ≫ R, r and thus we obtain If we consider d the reference pathlength, and ∆d the change in pathlength, then the change in phase, ∆φ, associated with a change in pathlength at wavelength λ is given by: By integrating the phase factor over ϕ, we obtain a prefactor which can be used to calculate the harmonic spectrum at a detector, a distance R away from the axis of rotation [29]: with where j l are spherical Bessel functions. The magnitude of −2π Rr dλ determines how many terms are needed in the summation. For example, for d = 24 cm, R = r = 0.1 mm and λ = 10 nm we required 26 terms in the summation for convergence.
For theoretical simulations, it is more appropriate to consider the intensity of the laser, I, rather than the radius of the laser focus, r, and thus we rewrite equation (6) as a function of I. For a 2D Gaussian focal spot with beam waist, r 0 (the intersection of the gas jet and the laser beam), the radial distance, r, at which an intensity I is observed is given by, where I max is the peak intensity of the laser. Thus, changing from r to I in equation (6), we obtain, Following [30], in calculations which solve the TDSE, the electric field E(ω, I) is proportional to the dipole acceleration,d(ω, I). Substituting the calculated dipole acceleration at a single intensity, I into equation (9), and taking the square modulus, we obtain the harmonic spectrum at detector radius R, given by, To obtain the full, focal-volume averaged harmonic spectrum, we take a discrete summation over the radius, R, of the detector. In theory, the inclusion of the phase accumulated during the propagation of harmonics from source to detectorthat is, the term A(R, I)-is absolutely crucial for the accurate determination of the measured HHG spectrum. However, as we shall see, it is useful to consider an alternative approach wherein we neglect the relative phase of harmonics originating a different points in the focus, or arriving at different points on the detector. Here, we only use the weight for an intensity I in a 2D-Gaussian focal spot, giving the simplified expression In equation (11), the beam waist r 0 is common to all intensities and can be considered a scaling factor for the total harmonic yield: In practice, we replace the integral over intensities in equations (10) and (12) with a simple quadrature, taking a weighted sum over a discrete set of intensities. This requires the selection of a suitable, minimum intensity, I min and intensity increment ∆I. These choices will be discussed in more detail in sections 3 and 4.

Time-dependent R-matrix theory
The RMT technique is a non-perturbative, ab initio method for solving the time-dependent Schrödinger equation (TDSE) taking full account of the multi-electron effects, obtaining an accurate picture of strong-field attosecond processes [31][32][33][34][35]. The method employs the well known R-matrix approach of dividing space into two distinct regions [36]. The inner region (close to the nucleus) describes a multielectron system taking into full account electron-electron interactions while the outer region, in which the ejected electron is far from the residual ion, describes an effective single active electron with long-range potentials, neglecting the effect of electronexchange. The TDSE is solved separately in the inner region and the outer region, using the most appropriate numerical method in each region. In the inner region, a B-spline basis is used to determine the multielectron wavefunction, whereas the outer region employs a finite difference method. The wavefunction is then matched explictly at the boundary between the two regions. This approach allows efficient use of computer resources: keeping the region with complicated, computationally demanding calculations small and employing a simpler approach in the large outer region further from the atom [27].
In RMT calculations, the dipole acceleration cannot be computed easily (as it is prohibitively sensitive to the description of atomic structure) [30] but, by using the relationship between acceleration, velocity and displacement, the harmonic spectrum can be expressed in terms of the dipole velocity or the dipole length: Consistency between length and velocity form spectra is used as a test for accuracy, and in this paper the velocity form is shown. The time-dependent expectation value of the dipole velocity operator isḋ where Ψ(t) is the time-dependent wavefunction of the system. The single-atom harmonic spectrum is then: where ω is the photon energy andḋ(ω) is the Fourier transform ofḋ(t).
We may then replace the dipole acceleration term in equations (10) and (12) with the frequency-scaled singleatom dipole-velocity, equation (15), to compute the intensityaveraged harmonic spectrum.

Xenon irradiated by a 1800 nm pulse
To validate the approach, we first consider the experimental scheme outlined in [17] to determine how the singleatom/single-intensity harmonic spectra relate to the experimental and intensity averaged results. In this scheme, xenon is irradiated by a linearly-polarised, 1800 nm laser field of duration 8 fs and peak intensity 1.9 × 10 14 W cm −2 . The experiment uses a tenuous gas jet in order to minimise propagation effects. We therefore can restrict our consideration of the macroscopic effects to the intensity volume, and treat the focal volume as a focal spot, with a corresponding 2-dimensional distribution of intensities (as employed in equation (8)).
In the RMT calculations, the interaction region comprises an inner region of 25 atomic units (a.u.) and an outer region of 1145 a.u.. An absorbing boundary is implemented beginning at 800 a.u. from the origin to prevent unphysical reflections of the wavefunction from the outer-region boundary. The finitedifference grid spacing in the outer region is 0.08 a.u. and we propagate the wavefunction using an eighth order Arnoldi scheme with a time step of 0.01 a.u.. The non-relativistic atomic description of Xe includes all available electron emission channels up to a maximum angular momentum of L max = 180 coupled to three possible configurations of the residual Xe + ion-4d 10 5s 2 5p 5 , 4d 10 5s 1 5p 6 and 4d 9 5s 2 5p 6 .

Results
Before calculating an intensity averaged result, it is important to first determine both the region of interest in the harmonic spectrum, and the range of intensities required to capture the behaviour in this region. The calculated cutoff energy for the laser parameters given is 196 eV. We note that the purpose of this experimental setup is to examine the giant resonance in xenon, and therefore in this instance, we are most interested in the harmonic yield at photon energies around 100 eV. Examining the single (peak) intensity RMT result, figure 1 shows good agreement between experimental results and the single intensity calculation up to around 100 eV, beyond which the results differ. If we restrict our concern to energies above 100 eV, then it is apparent from figure 1 that the lowest intensity which contributes at 100 eV is I min = 0.5I max . Further calculations show that in fact, including intensities above I min = 0.6I max is sufficient to capture the influence of focal volume averaging in the energy range of interest. Figure 2 shows RMT results obtained at a single (peak) intensity and two different focal volume averaging methods. The 'naive model', equation (12), simply accounts for the range of intensities in the laser focus. The 'full model', equation (10), also accounts for the relative phase of harmonics originating at different points on the focus and arriving at different points on the detector. Previously, it was thought that the discrepancy between the RMT, single-intensity result and the experiment was due to limitations in the atomic structure description used in the RMT calculations [37]. Calculation of a coherent intensity average improves the agreement with experiment at energies beyond 100 eV, suggesting instead that this difference may be due to the macroscopic intensity effects in experiment. Figure 2 shows that the main effect of performing an intensity average of the RMT results is to reduce the yield near the harmonic cutoff. The 'full method' (with I min = 0.6I max ) gives the best agreement with experiment, which is to be expected as this model is the most physically complete. However, when considering the naive model, this intensity range (I min = 0.6I max ) actually worsens agreement with experiment. There are two potential reasons for this: first, the weights for a 2D Gaussian focal spot used to calculate the intensity average (that is, the term in square brackets in equations (10) and (12)) favour the lower intensity contributions, and so as more intensities are included in the average, they contribute more. Second, in experiment, the focusing geometry is often set up such that the short trajectory dominates whereas the single-atom calculations include both long and short trajectories [39,40]. This long trajectory contribution, which dominates at lower photon energies, may explain the higher yield seen in the theoretical results.
To obtain a coherent average that had converged with respect to the intensity increment, it was required to use ∆I = 3 × 10 11 W cm −2 . For this increment and with I min = 0.6I max , 254 calculations are required which is computationally demanding and not often feasible. However, if we are Shown is experimental data [17] (thick black line) alongside data obtained from an RMT calculation for several different IR intensities, expressed as a fraction of the peak intensity. Shown is data obtained from an RMT calculation at the peak intensity (red line), experimental data [17] (thick black line) and data which has been averaged coherently using both the 'full method', equation (10), and the 'naive method', equation (12). To aid comparison, the RMT data sets have been normalised to match experimental yield at the cutoff. The data presented may be found at [38]. seeking only a picture of what differences between theoretical and experimental results are due to the laser intensity, it may be possible to include fewer intensities in the average to identify the basic behaviour. Figure 2 also shows focal volume averages calculated using a much smaller range of intensities (I min = 0.95I max ), from which we still observe the influence of the intensity averaging effect. This smaller range (and ∆I = 3 × 10 11 W cm −2 ) requires the calculation of only 32 individual intensities, which has a significantly smaller computational cost. For this intensity range, the 'naive method' gives good agreement with experiment across the spectrum. However, the 'full method' does not. In particular, there appears to be a 'suppression' of the harmonic yield at lower photon energies. This behaviour is due to the exclusion of the lower IR intensities: the yield for low IR intensities is much larger at low photon energies than at the cutoff (as demonstrated by figure 1) and the weights used in the averaging are inversely proportional to the intensity. Thus lower intensity contributions dominate at low energies.
However, we are most interested in the impact of taking an intensity average at the cutoff, and it is clear that in the energy range of interest, the 'naive model' results agree well with the 'full model', even when only a small range of intensities is included (I min = 0.95I max ). We reiterate that our goal is not to provide quantitative predictions of experimental spectra, and estimating this gross effect of intensity averaging at the cutoff even with the crude, approximate method we propose is valuable.
Another factor that must be considered is the increment in intensity, ∆I. Including more intensities (small ∆I) increases the computational cost. However, if ∆I is too large, then the harmonic phase may vary significantly between successive intensities used in the average and impact the result. Figure 3 shows coherently averaged spectra using both the naive and full model for various increments in intensity. In this specific set up, when using the naive model, the spectrum converges with respect to ∆I for a step size of 1 × 10 11 W cm −2 . However, figure 3 demonstrates that we can increase the intensity increment and still obtain spectra which capture the effect of using a coherent intensity average. This is especially important when calculations are computationally demanding as fewer will be required. The full model, equation (10) is much more sensitive to the intensity increment than the naive model due to the additional phase factor involved.
The proposed methods for modelling the intensity-volume effect account differently for the geometry of the experimental set-up. Crucially, the beam radius at the focus, r 0 must be considered. Often this is reported in the description of the experimental set-up. Figure 4 shows that the choice of r 0 has a significant impact on the coherent average obtained using the full model, equation (10). Here, for an intensity range I min = 0.6I max , a beam radius of r 0 = 0.0564 mm (corresponding to a focal spot of area 0.01 mm 2 ) gives the best agreement with experiment in the cutoff region. However, if a different range of intensities was used, then a different value of the beam radius may be more suitable. For example, in figure 2, a beam radius of r 0 = 0.1262 mm (corresponding to a focal spot of area 0.05 mm 2 ) was chosen for the full method using an intensity range of I min = 0.95I max . By contrast, for the naive model, equation (12), r 0 may be considered a scaling factor common to all intensities, and thus does not change the results, especially once they are normalised. The 'naive method' thus provides a good picture of the impact on experimental results without requiring the finer details of the experimental setup to be known, and removes the need for varying multiple parameters (the beam waist and intensity range) in lockstep. This may make the naive method more useful for making predictions.

Argon irradiated by a 780 nm pulse
As a second demonstration, we consider another atomic target with a different laser regime: argon irradiated by an eight cycle, 780 nm laser pulse at several peak intensities. Experimental data from [9] is used as a comparison. Shown is data obtained from an RMT calculation at the peak intensity (red line), experimental data [17] (thick black line) and data which has been averaged coherently using the full model, equation (10), for different values of the beam radius, r 0 . (I min = 0.6Imax, ∆I = 3 × 10 11 W cm −2 ).
In the RMT calculations, the argon atom is described with an inner region of 20 a.u. and an outer region of 1227 a.u. The finite-difference grid spacing in the outer region is 0.08 a.u. and the wavefunction is propagated with a time step of 0.01 a.u.. The atomic description of argon includes all 3s3p 6 ϵℓ and 3s 2 3p 5 ϵℓ ionisation channels up to a maximum total angular momentum of L max = 120. The inner region is constructed using a set of 40 B splines.

Results
The Ar experiment was set up ostensibly to minimise the effect of phase mismatch and reabsorption of high harmonic radiation, which enabled the observation of the single atom response [9]. Our approach allows a comparison of the experimental results with both single intensity and coherently averaged results to determine which features of the single-atom/single-intensity spectra can be seen in the macroscopic experimental results. The cutoff energy calculated using the experimental laser parameters is 78 eV but appears (see figure 5) much higher at around 88 eV. To aid comparison, the peak intensity used in the RMT calculations is increased by 6 × 10 13 W cm −2 from that reported in experiment, increasing the cutoff energy to 88 eV. Figure 5 shows RMT results obtained at a single intensity and intensity averaged results. The intensity averaged results were calculated over two ranges: a range of intensities down to I min = 0.95I max and a range of intensities down to I min = 0.7I max , both with ∆I = 1 × 10 12 W cm −2 and a beam radius of r 0 = 0.1784 mm. These two intensity ranges were chosen to demonstrate the compromise between a low computational cost (I min = 0.95I max ) and clarity of results in the intended experimental observable (the Cooper minimum at around 53 eV).
In this regime, the 'naive method' with I min = 0.95I max does not cause any drastic change to the single-intensity spectrum. In particular, the shape and position of the Cooper minimum are largely preserved in the averaged spectrum. With our stated goal of predicting which single-atom effects may 'survive' intensity averaging, this result is an apparent endorsement of this simple tool. Indeed, using the 'full' model to calculate the intensity average improves agreement with experiment at lower energies and at the Cooper minimum for both intensity ranges considered, but the effect of the intensity averaging is well captured by even the smallest set of calculations. For this small set of intensities (with I min = 0.95I max and ∆I = 1 × 10 12 W cm −2 ) only 22 calculations are required: a notably smaller computational cost than I min = 0.7I max which requires 125 calculations.
The experiment was set up specifically to investigate the Cooper minimum in the harmonic spectrum. Considering this, we are most concerned with agreement between simulation and experiment at energies around the Cooper minimum. The interpretation of the minimum in the experimental results is already well supported by the single-atom, singleintensity results. However, if instead we were predicting the appearance of the minimum based on simulations alone, it would be useful to know that this feature would not be 'washed out' in more realistic experimental conditions.

Conclusion
The information encoded in the radiation emitted during HHG necessarily couples microscopic (single-atom/molecule) and macroscopic processes. If the goal of an experiment is to determine electronic dynamics at the single-atom level, the relevant signatures in the harmonic spectrum may be obscured by the two key macroscopic processes: propagation through the medium and the intensity volume effect. Whilst experimental methods have been developed to minimise the propagation effects (for instance, using tenuous gas jets), there is no method of producing a laser beam at a single uniform intensity. Since ab initio methods model a single atom irradiated by a laser beam at a single intensity, the microscopic details of these harmonic spectra are accessible. In this paper, we have outlined and demonstrated a method which can be used to help decode the single atom response from the macroscopic spectrum and/or predict the appearance of single-atom features in the observed, experimental spectrum.
We have outlined a method of approximating the intensity volume effect without performing thousands of single-atom calculations over the full-range of intensities in the laser focal spot. While the full physical model must include the frequency-dependent phase of harmonics spanning a large range of intensities (between the peak intensity, I max and I min ≈ 0.3I max in our calculations) and also consider the geometry of the detector, we have found that the effect is approximated well using a significantly smaller range of intensities (I min = 0.95I max ) and neglecting the relative phase of harmonics.
The simplified method we propose extends an already well-established method for focal volume averaging of photoelectron emission spectroscopy calculations [20]. Exploratory RMT calculations have also shown that the method is appropriate for the determination of other experimental observables of interest in attosecond science, such as the absorption spectrum which, like the HHG spectrum, depends on the time-dependent expectation value of the atomic/molecular dipole [13,41]. Although this suggests the method is quite general, we note that the possible range of experimental schemes-invoking laser pulses and target properties over a vast parameter space-necessitates care in the application of the simplified approach. In particular, the use of a small number of sample intensities with a relatively large intensity increment was made possible in our calculations because the harmonic phase changes slowly with intensity in the cutoff region. Different wavelengths or peak intensities may require a more fine-grained approach. Additionally, the energy range of interest must be considered, as this may impact the range of intensities required in the focal volume average. If we are interested in the cutoff, then the 'naive' method with a small range of intensities (I min = 0.95I max ) is suitable because the lower intensities will not contribute at high energies. However, if the feature of interest lies at low energies, then many more intensities may contribute.
We reiterate that the goal here is not to accurately reproduce experimental results with high-fidelity-although this has been used as an assessment of validity here. Rather, these tools may be used to assess which single-atom effects may 'survive' in the experimentally measured signal or, conversely, to determine which features of measured spectra can be attributed to single-atom processes. In this paper we have applied the intensity averaging methods to RMT calculations. However, in principle, the methods could be used to calculate a focal volume average of any ab initio results, where the computational cost prohibits performing a full intensity scan. We note that the underlying ab initio calculations do not change whether or not the simplified (neglecting phase) or full (including phase) approach is used, and so the computational expense involved in exploring the parameter space to determine appropriate values of I min and ∆I is mitigated.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.5281/ zenodo.7886420.