of high

The demanding experimental access to the ultrafast dynamics of materials challenges our understanding of their electronic response to applied strong laser fields. For this purpose, trapped ultracold atoms with highly controllable potentials have become an enabling tool to describe phenomena in a scenario where some effects are more easily accessible and twelve orders of magnitude slower. In this work, we introduce a mapping between the parameters of attoscience platform and atomic cloud simulators, and propose an experimental protocol to access the emission spectrum of high harmonic generation, a regime that has so far been elusive to cold atom simulation. As we illustrate, the benchmark offered by these simulators can provide new insights on the conversion efficiency of extended and short nuclear potentials, as well as the response to applied elliptical polarized fields or ultrashort few-cycle pulses.

The demanding experimental access to the ultrafast dynamics of materials challenges our understanding of their electronic response to applied strong laser fields.For this purpose, trapped ultracold atoms with highly controllable potentials have become an enabling tool to describe phenomena in a scenario where some effects are more easily accessible and twelve orders of magnitude slower.In this work, we introduce a mapping between the parameters of attoscience platform and atomic cloud simulators, and propose an experimental protocol to access the emission spectrum of high harmonic generation, a regime that has so far been elusive to cold atom simulation.As we illustrate, the benchmark offered by these simulators can provide new insights on the conversion efficiency of extended and short nuclear potentials, as well as the response to applied elliptical polarized fields or ultrashort few-cycle pulses.
Over the last three decades, progress in laser technologies has led to significant advances in our ability to manipulate and understand electron dynamics on their natural attosecond (10 −18 s) timescale [1][2][3][4].This has triggered the development of a huge range of tools for probing and controlling matter, which includes high harmonic spectroscopy [5], laser-induced electron diffraction [6,7], photoelectron holography [8,9], attosecond streaking [10,11], and reconstruction of attosecond harmonic beating by interference of two-photon transitions [12,13], to name only a small selection.One of the most prominent processes underlying some of these successes is high harmonic generation (HHG), a highly non-linear phenomenon where a system absorbs many photons of the driving laser and emits a single photon of much higher energy.
The experimental realization and interpretation of these complex processes have been guided most often by simplified theoretical descriptions that still capture the main properties of the dynamics.In the field of attoscience, simplifications such as only considering one or two active electrons, disregarding the interaction of the ionized electron with the parent ion during its propagation in the continuum, or performing saddle-point approximations have provided valuable quantitative predictions concerning HHG [14][15][16] and other phenomena [17][18][19][20].Additional experimental regimes do, however, require a more complete description of the system, including those where multielectronic processes [16], * javier.arguello@icfo.euor Coulomb nuclear potentials play a key role [21] (see Refs. [9,22] for reviews).This has motivated an intense development of analytical and numerical methods aimed at pushing current computing capabilities.
To circumvent this computation complexity, analog simulation has become an enabling tool to access phenomena with highly controllable devices, whose temporal and spatial scales are more favorable to measure than those naturally present in attosecond physics [23][24][25].Experimental advances in the engineering of interactions now foster the simulation of quantum chemistry phenomena, such as molecular geometries [26,27], vibronic calculations [28,29], or the presence of conical intersections [30,31].Recently, this experimental benchmark of chemical processes has accessed the response of an atomic system to strong pulses in the regime where the energy imparted by the simulated laser field is strong enough to ionize the target atoms.This is the case of experiments where the incoming field is simulated by the curved shape of a waveguide [32], or by a shaken potential applied on a neutral atomic species [33], following early proposals where the correspondence to the static frame was introduced [34][35][36].Going beyond the ionization regime, Sala et al. [36] noted that controllable Zeeman or Stark shifts can also give access to the relevant regime of HHG in atomic simulators.However, the simulation of the HHG spectrum has remained so far elusive in these experiments where photons or neutral atoms mimic the oscillating electron, as no radiation is emitted by the associated neutrally-charged simulated dipole.In this work, we show that current atomic platforms offer a unique opportunity to access and measure the emission spectrum of HHG through absorption measurements.Furthermore, it simulates the physical response of a single-atom target.This is in contrast with real experiments, where thousands of atoms are simultaneously driven to collect enough photons to resolve the spectrum, which challenges phase-matching conditions when a large ionization occurs under strong fields.The simulator thus provides an important bridge between experiments and theory, offering controllable experimental access to complex systems that otherwise can only be theoretically approximated.
The text is structured as follows.First, we present key concepts on attoscience physics and ultracold atom simulators that is of special interest for readers newly exposed to either of these areas.In particular, Sec.I discusses the status of current investigation in attosecond science and motivate the different regimes that result from the frequency and strength of the laser field that drives the process.There, we introduce the regime of HHG, and some of its distinctive features.Next, Sec.II presents current applications of analog quantum simulation, focusing our attention on the flexibility offered by atomic systems subjected to tunable light fields.There, we introduce an atomic simulator capable of mimicking the dynamics of an electron exposed to a strong oscillatory laser field, and derive a mapping between the experimental parameters of the simulator and the relevant units encountered in attosecond science.We also devise a protocol to access the simulation of the emission yield in HHG, highlighting the range of parameters where this correspondence is valid.As an illustrative example, we show how one can use the simulator to study the effect that short pulses and the ellipticity of the incoming field have on the efficiency of HHG.Sec III further presents details on the experimental choice of atomic species and laser pulses that are needed to simulate specific targets of common studies in attosecond science, and discuss the main sources of errors that the experimental implementation would encounter, which we numerically benchmark.In Sec.IV, we present an outlook of the venue for exploration that the proposed simulator opens.
Although HHG is a process that has been observed in many different targets, such as atoms [14,53], molecules [54], solid-state systems [55,56] and nanostructures [4], in this work we will focus on simulating HHG process in atoms.We shall now discuss the conditions required to generate high-order harmonics from atomic targets by introducing relevant parameters.The first one is the so-called multiquantum parameter K 0 = I p /(ℏω) [22], where I p is the ionization potential of the corresponding atom, and ω is the frequency of the driving field.This parameter provides an estimate of the minimal number of photons required to ionize the electron from the atomic ground state.Here, we are particularly interested in the low-frequency limit, where K 0 ≫ 1, indicating that ionization requires a significant number of photons.Furthermore, whether a multiphoton absorption process occurs or not also depends on the amplitude F 0 , of the applied field.This motivates the introduction of the Keldysh parameter γ = I p /(2U p ) [14,22,57], where U p = e 2 ξ 2 0 /(4mω 2 ) is the ponderomotive energy, i.e. the average kinetic energy of an electron with mass m and charge e in the presence of an oscillating laser electric field of amplitude ξ 0 [see Fig. 1(a)].
In the case K 0 ≫ 1, the multiphoton regime is observed when γ ≫ 1, indicating that the field only slightly perturbs the atomic potential [see Fig. 2 (a)].The regime of interest for HHG processes corresponds to γ ≲ 1, known as the tunneling regime, where the external field force is comparable to the atomic potential, which gets distorted and forms an effective potential barrier through which the electron can tunnel out [see Fig. 2 (b)].Finally, and for the sake of completeness, we have the regime of over-the-barrier ionization that, typically [58], happens when the electric field reaches a critical value that makes the barrier maximum to coincide with the energy level of the electron ground state [see Fig. 2 (c)].
Within the tunnelling regime, the three-step model, also referred to as the simple-man's model [39,59,60], provides a powerful picture of the underlying electron dynamics behind HHG.The steps within this model are as follows: the electron (i) tunnels out from the parent atom through the barrier formed by the Coulomb potential combined with the dipole interaction of the field, (ii) oscillates in the continuum under the influence of the laser electric field and, if it passes around the nucleus' vicinity, (iii) can recombine back to the ground state emitting harmonic radiation.The energy of the emitted radiation upon recombination depends on the electronic kinetic energy at the moment of recombination and the ionization potential of the atom.However, the maximum kinetic energy that the electron can acquire from the field during its propagation is limited, leading to a natural cutoff frequency ω c for the highest harmonic order in HHG, which is determined by ℏω c = I p + 3.17U p .The nonlinear character of the HHG process can be understood The resulting oscillation of the charge emits radiation of characteristic harmonic frequencies (inset in blue).The same emission yield can be simulated on the simulator (b), where the characteristic frequencies are retrieved through absorption images of an atomic gas (green) that is trapped by a laser potential (red), and addressed by an external magnetic gradient that is tuned over time (brown).from its main features, namely (i) a strong decrease in the low-order harmonics amplitude, (ii) a plateau, where the harmonic yield is almost constant, and (iii) a sudden cutoff, given by the above presented classical formula [39].
Based on the previous discussion, it is evident that the polarization character of the driving field can have significant consequences on the generated harmonic radiation.For instance, when an elliptically polarized driver is considered, the ionized electron may miss the parent ion, resulting in the absence of the recombination event [59].This phenomenon has been extensively studied in both experimental [61][62][63][64][65][66] and theoretical [67,68] works, demonstrating a reduced HHG conversion efficiency, i.e. the ratio between the outgoing and incoming photon fluxes, as the driving laser ellipticity increases.As an alternative strategy, a combination of two drivers that have different ellipticities and frequencies can be used to generate bright phase-matched circularly-polarized high harmonics, as shown in Refs.[69][70][71][72][73].
While we will focus on HHG, strong laser-matter interactions can also give rise to other fascinating phenomena, including Above-Threshold Ionization (ATI) [3,[74][75][76][77] and Non-Sequential Double Ionization (NSDI) [59,[78][79][80][81][82].In ATI, an electron is ionized by the strong-laser field, surpassing the ionization threshold of the corresponding atom by absorbing more photons than the ones required for ionization.The typical observable measured in ATI is the photoelectron spectrum, which exhibits distinctive peaks at electron kinetic energies separated by the energy of a single photon of the driving field [74,76].In the non-perturbative (tunneling) regime, these peaks form a plateau that extends over electron energies on the order of 10 U p [76,83].On the other hand, NSDI occurs when an ionized electron undergoes rescattering with its parent ion, resulting in the ionization of a second electron.This phenomenon is reflected in the ion yield, which exhibits a knee structure at a specific intensity threshold.The presence of this distinctive feature signifies a sudden change in the energy distribution of the emitted electron pairs, indicating a transition from sequential to non-sequential ionization processes [78,79].
Here, we demonstrate the capability of analog simulators to accurately replicate the key characteristics of the HHG processes in atoms.Specifically, it recovers the main features of a typical HHG spectrum -a plateau extending for few harmonic orders followed by a cutoffwith an harmonic yield that reduces for increasing values of the ellipticity.We study how the conversion efficiency of the harmonics depends on the values of K 0 and γ, and discuss how the simulated spectrum can be measured in practice for analog simulators, also providing estimates on relevant quantities towards feasible experimental implementations.

II. ANALOG QUANTUM SIMULATION
The numerical simulation of chemical problems generally requires to describe many electrons that interact with external fields, the nuclei, and among themselves through Coulomb interactions.Even if one considers the nuclear positions fixed due to their larger mass (the Born-Oppenheimer approximation [84,85]), this is an extremely challenging task, as the associated Hilbert space grows exponentially with the number of electrons.
Over the last few decades, an alternative route to study electronic problems has emerged, based on using quantum devices that can better capture the complexity of the system.This idea was first proposed by Feynman as a way of preventing the exponential explosion of resources of quantum many-body problems [86], and later formalized by Lloyd [87].Complementary to current efforts in digital simulation [88][89][90] (where the problem is first encoded as qubits and gates addressed by a generalpurpose quantum device), simulators based on ultracold atoms have become an enabling tool, already addressing quantum matter phenomena that the most advanced classical computers cannot compute [91,92].
At low temperatures, the interaction among atoms can be highly engineered with external lasers, which allows one to induce a rich variety of effective Hamiltonians on a highly controllable platform [25,93].Early experiments dealt with condensed matter problems [94][95][96], detecting a transition between the superfluid and Mott insulating phases of effective Hubbard models.More recently, atomic simulators have experimentally addressed problems related to high-energy physics, such as Gauge theories, both in lattice geometries [97][98][99] and the continuum [100].An exciting perspective consists of extending the benefits of analog quantum simulation to the field of chemistry and the response of atoms and molecules to external fields.Soon after the first experimental realization of bosonic gases [101,102], the mapping between degenerate atomic gases and single-electron dynamics was noticed [34,103].As compared to a real material, where many target atoms are present, here the simulated electron moves in a cleanly isolated environment.Furthermore, the typical energy-scales of these experiments are in the range of kHz-MHz, which provides a temporal magnification of the simulator, where attosecond pulses are associated to convenient timescales of µs-ms, i.e. 10 9 ∼ 10 12 times slower.The additional tunability and accessibility of these simulators can thus offer a complementary tool to investigate ultrafast phenomena.
Using the Kramers-Henneberger correspondence, the shaking of the optical trap can mimic the effect of an external force [35,36,104], which has allowed for the recent experimental simulation of ATI processes using a bosonic gas of 84 Sr [33], where the kinetic energy of the simulated ionized electrons is accessed through a time-offlight measurement (TOFM) of the atoms emitted during the shaking.Extending this approach to the HHG regime is however challenging, as the needed inertial force corresponds to strong nuclear potentials (see Appendix A), and the associated photonic spectrum is neither emitted during the oscillation of neutral atoms, nor directly revealed by the TOFM.The present contribution advances the study of HHG simulation in different directions.In particular, (i) we establish a correspondence between the experimental parameters in the simulator and the K 0 and γ parameters conventionally used in attosecond science, presenting specific configurations associated with standard choices of atomic targets and ionizing pulses; (ii) we show that the region of maximum HHG conversion efficiency can be accessed in atomic simulator platforms where the external pulses are simulated by existing Stark and Zeeman shifts; (iii) we introduce a scheme to experimentally access the generated harmonic spectra through absorption measurements of the atomic gas, and characterize the main sources of error.

The simulator
Whenever multielectron processes can be neglected, the underlying physics of strong field processes can be described by a single-active electron.The dynamics of this electron is dictated by the Hamiltonian which accounts for its kinetic energy, the attractive nuclear potential, and the incident laser field, respectively.Following the dipole approximation, the interaction with a laser field linearly polarized along the x axis writes as, V F (r, t) = F (t)x, where F (t) = F 0 f (t) sin (ωt + φ) and f (t) = sin 2 [ωt/(2n c )] is the carrier-envelope for a pulse with n c cycles and wave phase φ.Here, F 0 = eξ 0 is the maximum force imparted by the field on the electron.In this proposal, the dynamics of the electron inside the nuclear potential is simulated by an atomic Bose gas optically trapped by a laser beam, whose spatial profile can be highly engineered, even dynamically, with the use of spatial light modulators [106,107] or digital mirror devices [108,109].Working in the single-active electron approximation, this can allow one to engineer an effective potential for the valence electron based on an ab initio calculation of the multielectronic species, which can even provide quantitative agreement on the emitted spectrum, as theoretically studied in Refs.[110,111].As a first experimental step for the simulator, here we will consider the natural Gaussian profile of a laser beam of waist r 0 and depth V 0 in a one-dimensional system Similarly to the widely used atomic units, it becomes convenient to define the natural units of this system as In these units (that we denote as [•] along the text) the nuclear potential is fully characterized by the width of the trap, as it follows from the relation [V 0 ] = [r 0 ] −1 , and approximates a quantum harmonic oscillator whenever the associated zero-point motion of the oscillator is smaller than the waist of the beam, [r 0 ] ≫ 1.
In Fig. 3(a), one can see that the simulated Gaussian potential (dashed blue line) matches at short distances (x/r 0 ≪ 1) the widely used soft-core nuclear potential of the form, V sc (x, r 0 ) = −1/ x 2 + r 2 0 (red line, see Appendix B for further details) [17,112].To highlight this connection between the simulator and the soft-core potential, in Fig. 3(b) we calculate the ionization energy associated with both nuclear potentials for different values of r 0 , observing discrepancies smaller than 10% for [r 0 ] > 1. Focusing on the average width, ⟨x 2 ⟩, of the ground-state, the long-range polynomial scaling of the soft-core potential leads to more extended eigenstates than the short-range Gaussian trap [113,114] (see inset), which can be experimentally explored by further shaping the laser beam.As compared to a real target with a fixed nuclear potential, the tunability offered by the simulator allows one to benchmark how the predictive power of conventional numerical methods is influenced by the range of the studied nuclear field.For example, this is relevant for the strong field approximation, where the Coulomb potential along the excursion time of the electron is often numerically disregarded [14].
To mimic the effect of an incoming electric field under the dipole approximation, the atomic cloud is subjected to a time-dependent linear energy gradient, V F (r, t), which can be created by an optical Stark shift that is proportional to the intensity of an applied off-resonant laser field.In current platforms, a linearly-varying intensity can be created with a spatial light-modulator [91], an acousto-optical device, or simply by using the slope of a Gaussian beam whose intensity and position can be dy- namically adjusted.For atomic levels that are sensitive to magnetic fields, one can alternatively rely on Zeeman shifts induced by linear magnetic field gradients created and modulated with current-carrying coils [115][116][117], as represented in Fig. 1(b).

Simulation of HHG emission yield
While there is not an analog equivalent to the photons emitted during HHG, its emission yield can be accessed through the time-dependent dipole moment d(t) = ⟨ψ(t)| x |ψ(t)⟩, or its associated time-dependent dipole acceleration, d a (t) = ⟨ψ(t)| ẍ |ψ(t)⟩, where |ψ(t)⟩ denotes the atomic state at time t.In HHG experiments, the spectrum of energies for photons emitted over the duration of the laser pulse is characterized by the Fourier transform, d a (Ω) = dt d a (t)e −iΩt [118].
In Fig. 4 we use a time-dependent Schrödinger equation (TDSE) to calculate the dipole acceleration d a (Ω) in one spatial dimension.From Fig. 3(b), we choose parameters compatible with the Hydrogen ionization potential [r 0 ] = 1.26, so that the corresponding dissociation energy is [I p ] = 0.5, and a 6-cycle pulse with an associated laser field of wavelength 800 nm and intensity 10 14 W/cm 2 .There, we see that the short-range correspondence between the soft-core (red line) and Gaussian potentials (blue) translate into a qualitative agreement of the resulting harmonic spectrum.
To access this quantity in the simulator, one option consists of measuring the atomic spatial density and following the Ehrenfest theorem [119] where ∇V n (x, r 0 ) denotes the gradient of the trapping potential.For the Gaussian trap expressed in Eq. ( 2), this gradient can be approximated as, an absorption measurement allows one to access d a (t) by quantifying the asymmetry in the atomic population for positions at a characteristic separation r 0 / √ 2 from the center of the trap.By repeating the measurement at different times, one can recover the emission spectrum by Fourier transforming d a (t).
As an alternative, one can also access the timedependent dipole velocity, d v (t) = ⟨ψ(t)| ẋ |ψ(t)⟩ , whose spectra is related to the emission yield as, d v (Ω) = id a (Ω)/Ω, in the case of finite pulses [120].Interestingly, the velocity components of the atomic cloud at a given time can be conveniently accessed in the simulator through a TOFM, where the nuclear trap is suddenly released and the gas expands ballistically.Once the gas has expanded for a time τ beyond the initial size of the cloud, the velocity component, v, of the state of interest is associated to an absorption detection of the cloud at a distance x = vτ from the initial trap.Given that these separations are much larger that the size of the original cloud, this approach improves the accuracy of the reconstructed emission yield for a given spatial resolution.
In Fig. 5(a) we calculate the conversion efficiency [121] η for one of the harmonics in the plateau region, a fixed nuclear potential [I p ] = 0.5, and different values of In order to interpret the observed regimes, both K 0 and γ describe a complete set of parameters that characterizes the response of the simulated system to the oscillatory field.For example, expressed in natural units, the last harmonic Ω c = ω c /ω that becomes accessible below the cutoff energy writes as As expected, a higher yield appears above the continuous line 5 ○ where the condition, Ω ≤ Ω c , is satisfied.To enhance the conversion efficiency, it is also desirable not to be too deep in the tunneling regime.An optimal situation occurs when the maximum tilt of the nuclear potential, F 0 r 0 , is comparable to and we observe that the region of the largest conversion efficiency follows this heuristic scaling (red line 1 ○).Regarding the mapping to attoscience, the dipole approximation followed in Eq. ( 1) requires that the magnetic component of the incoming field can be disregarded.This imposes an upper bound 2 ○ on the intensity of the field, which writes as [122] Preventing relativistic velocities on the accelerated particle (U p ≪ mc 2 ) also induces a lower bound on the Keldysh parameter  ○ follows the region where a higher efficiency is expected from Eq. ( 6), dashed lines 3 ○ and 4 ○ describe the HHG region (K0 ≥ 1 and γ ≤ 1, respectively). 5 ○ follows the frontier where the selected harmonic is below the cutoff frequency in Eq. 5.The dipole approximation 2 ○, associated to the upper bound in Eq. ( 7), lies above the represented range of values for K0.Arrows point towards the region of validity for each of these conditions.Circled red marker indicates the Gaussian configuration in Fig. 4, while the crossed marker indicates the point K0 = 6.31 and γ = 0.81, associated to [F0] = 0.171 and [ω] = 0.108.The simulated emission spectrum of the latter configuration is illustrated in panel (b) for incoming pulses with a number of cycles running from nc = 6 (red) to nc = 1 (green).To improve visibility, an arbitrary vertical shift between the curves has been applied.
which is less demanding than the dipole approximation along this region of greatest conversion efficiency in Eq. ( 6) for γ ≪ 1.In the natural units of the simulator, one can see that [c] ∼ 10 11 for the parameters studied in Fig. 5(a), which places this upper bound above the represented range of values for K 0 .
One should note that the harmonic radiation measured in in actual attoscience HHG experiment results from a collective phenomena, involving 10 10 ∼ 10 12 atoms that emit coherent radiation in phase.When comparing theory and experiment, it is therefore mandatory to include macroscopic propagation effects, which is a formidable computational task that is only addressed by a limited number of models [123].In this simulator, however, all atoms in the bosonic gas contribute to magnify the effects manifested by a simulated single electron.The resulting harmonic spectrum is thus unaffected by the additional phase-matching condition that is often encountered in attosecond science experiments, providing clean access to the simulated single-atom dipole acceleration.Furthermore, the more favorable energy, spatial and temporal scales in the simulator offer a benchmark that is less affected by the energy resolution, the uncertainty in the laser intensity and duration, or the limited dynamic range of spectroscopic measurements in ultrafast experiments [58].
The high tunability of the induced incoming pulse is also one of the advantages of the simulator, as compared with the fast high-intensity lasers that are typically needed in attoscience.For example, this allows one to easily simulate the response of the system to ultrashort pulses, a configuration that is otherwise demanding to access in real attosecond experiments.In Fig. 5(b), we show the simulated emission yield for pulses with different number of cycles.As the pulses get shorter, we observe that the plateau structure vanishes and the last harmonics disappear, even though the theoretical cutoff frequency ω c only depends on the ionization energy and ponderomotive energy, which are the same for all the curves.Intuitively, fewer interference processes can take place when the number of cycles of the pulse decreases, which reduces the number of harmonics that are visible on the emission spectrum as one approaches n c = 1 (green line).Focusing on this latter case, in Fig. 6 we simulate the response of the system to a single-cycle pulse (n c = 1), for different values of K 0 and γ.When focusing on one of the lowest harmonics (the 5th one), we observe that the region of largest emission still remains well described by the different conditions introduced in Fig. 5(a).
In addition to the linearly polarized fields considered so far, the simulator also allows one to induce oscillations in a second axis, when we extend the systems to two dimensions.By controlling the ratio between the two amplitudes, i.e. the ellipticity ε = F 0y /F 0x , one can induce elliptic [ε ∈ (0, 1)] and circularly-polarized fields (ε = 1).In HHG, the introduction of ellipticity in the laser beam leads to the deflection of the returning electron, causing it to deviate from its intended path toward the parent nucleus.This results in a decrease of the overlap between the wavepackets of the returning electron and the initial bound state, as it has been experimentally observed [61].In Fig. 7 we calculate the change in the conversion efficiency of different harmonics as we change from a linear field (ε = 0) to an elliptically-polarized field.As expected, we observe a decline in the intensity of harmonics as the ellipticity of ○ follows the region where a higher efficiency is expected from Eq. ( 6), dashed lines 3 ○ and 4 ○ describe the HHG region (K0 ≥ 1 and γ ≤ 1, respectively). 5 ○ follows the frontier where the selected harmonic is below the cutoff frequency in Eq. 5. Arrows point towards the region of validity for each of these conditions.The crossed marker indicates the point K0 = 6.31 and γ = 0.81, whose emission spectrum is depicted in green in Fig. 5 the laser beam increases, which is more pronounced for the lowest harmonics.

III. EXPERIMENTAL IMPLEMENTATION
At this point, it is worth exploring the feasibility of the experimental parameters associated with this implementation.For a fixed geometry [r 0 ] = 1.26 (associated to V 0 (K) 7 Li 84
[I p ] = 0.5) the needed nuclear potential scales as In the case of 84 Sr and a Gaussian beam with waist r 0 = 1 µm, the configuration K 0 = 6.31 and γ = 0.81 [illustrated in Fig. 5(b)] corresponds to parameters, F 0 = 1.55 • 10 −26 N, ω = 94.7 Hz, and V 0 = 7.22 nK • k B .Using moderate conditions of 1 Watt of 532 nm light shaped to give a linear intensity gradient across a 100 µm × 100 µm area, it is possible to achieve the needed values of F 0 for 84 Sr atoms, and even produce forces two orders of magnitude stronger.However, one can observe that the associated trap depth V 0 ∼ 10 nK, is well below the trap depth of around 10 µK used in previous attoscience simulators with Sr [33].
The trap depth sets an upper limit to the temperature of atomic clouds that can be trapped by the laser potential.To increase the range of allowed temperatures for the experiment, Eq. ( 9) indicates that more relaxed cooling conditions would benefit from choosing narrower beams and lighter atomic species.This is the case of 7 Li and r 0 ∼ 1 µm, where the parameters needed to simulate the same attoscience configuration are: F 0 = 1.86 • 10 −25 N, ω = 1.14 kHz, and V 0 ∼ 0.1 µK • k B .The associated critical temperature is then consistent with state-of-the-art experiments [124].As mentioned in the previous section, for magnetic ground states, the oscillatory force can be applied by a tunable magnetic field gradient.In the case of the |m F | = 2 hyperfine state of 7 Li, easily attainable gradients of 50 G/cm translate into the needed range of forces F 0 , or even one order of magnitude larger.
In Fig. 8, we show the scaling of the trap depth with the beam waist for 84 Sr (green) and 7 Li (blue lines).We also show that more relaxed cooling conditions would appear in the simulation of weaker ionization energies, as we illustrate with dashed lines for the Sodium first ionization energy, [I p ] = 0.19.Depending on the choice of atomic species and width of the trap, we observe a range of 3 orders of magnitude on the associated trap depth, which illustrates the flexibility that these simulators offer to access the regime of interest.

Further experimental considerations
The reconstruction of the dipole acceleration is also affected by experimental imperfections and limitations that need to be accounted for.Here, we numerically benchmark the main sources of error for the previous configuration K 0 = 6.31, γ = 0.81 and n c = 6, where a larger conversion efficiency is expected [crossed marker in Fig. 5(a)].To quantify the effect of these imperfections, we calculate the relative error in the determination of the local maxima of the spectra for frequencies below the cutoff frequency, which we highlight with orange dotted markers along Figs.9-12.FIG. 9. Emission yield extracted from the discrete Fourier transform of Nt = 300 temporal intervals uniformly distributed along the duration of the pulse (red dotted line), as compared to the converged spectrum (Nt = 3000, blue line).Orange dots indicate the positions of local maxima before the cutoff energy.In the inset, we show the relative error of these local maxima for increasing number of temporal divisions, Nt, where the choice Nt = 300 is highlighted with a red cross, and the rest of parameters coincide with Fig. 5(b).See Appendix C for details on the numerical calculation.
In the experiment, the time-dependent dipole acceleration is measured at a finite set of times, d a (t n ) n=1...Nt .The accuracy of the discrete Fourier transform then depends on the number of time points used in the reconstruction, that we consider to be uniformly distributed along the duration of the pulse.In the inset of Fig. 9, we calculate the relative error in the local maxima of the emission yield for different values of temporal divisions.We observe that a moderate number of time in- tervals N t ∼ 300 provide an error of ∼ 1% that is enough to resolve the cutoff energy (see red dashed line in Fig. 9).
Experimental errors in the time points used to measure, ∆t n = |t meas − t n |, would also manifest in the reconstructed yield.In the inset of Fig. 10, we calculate the relative error in the local maxima of the emission yield for different values of noise with standard deviation ∆t around the time intervals t n .We observe that moderate values ω • ∆t ∼ 10 −3 provide an overall relative error ∼ 10% in the local maxima that is enough to resolve the cutoff energy (see red dotted line in Fig. 10).Expressed in experimental units, this translates into a correct control of time in the order of ∼ 10 µs.Overall, we see that the highest harmonics are the most sensitive ones to errors in the timing and number of time intervals, as they need to capture the rapidly oscillating response of the emission spectrum.
To quantify the time-dependent dipole acceleration from absorption measurements one also needs to characterize the needed accuracy on the determination of d a (t).In the inset of Fig. 11, we calculate the relative error in the local maxima of the emission yield for different values of relative error in the measurement of the time-dependent dipole acceleration, ∆d a (t)/d a , where d a = T 0 |d a (t)|/T is the temporal average throughout the pulse.We observe that ∆d a (t)/d a ∼ 1% provides an accuracy in the emission yield that is enough to resolve the cutoff energy (see red dotted line in Fig. 11).We also observe that the highest harmonics, which have the smallest intensity, are the most affected by a limited resolution.Following the TOFM approach to access d a (t), typical imaging lengths, 0.1 − 1.0 mm, associated to expansion times τ ∼ 10 − 100 ms [125], can provide  the needed accuracy of 1% on the retrieved velocities for an atomic cloud that has an initial width of 1 − 10 µm before its ballistic expansion.Also, an error in the measured density of order ∆|ψ| 2 ∼ 10 −4 can be tolerated to reach the desired relative error of 1% in the extracted time-dependent dipole (see Appendix C).As the atomic cloud spreads over its different velocity components, this results in N ev ∼ 10 6 single-atom events to reach the needed accuracy on velocities and atomic densities during the TOFMs used to extract each instantaneous dipole acceleration.
Considering all these demands, one can estimate the overall experimental time needed to reconstruct the timedependent dipole acceleration.For an atomic cloud with N at ∼ 10 4 atoms [33] and an estimated running time of τ ∼ 1 s per experiment, collecting enough statistics translates into a reasonable running time τ • N ev • N t /N at ∼ 10 h.The major demand for resources comes from the repetition rate of the experiment and the required number of events needed to resolve the small variations in the absorption images.First experimental realizations are then specially suitable for the simulation of the first harmonics in the plateau region, as they are associated with larger conversion efficiencies, reducing the number of independent measurements that are required and, with that, the overall experimental time.
As a final remark, we have seen that thousands of atoms in an atomic cloud are used to simulate the state of a single electron.Atoms do, however, experience interactions at distances characterized by the scattering length, a s .In a mean-field approach, the effect of these interactions can be described by the Gross-Pitaevskii effective Hamiltonian that depends on the atomic density at each point of space [126] Ĥint (r, t) = Ĥ(r, t) where, g = 4πℏ 2 a s /m.In the inset of Fig. 12, we calculate the relative error in the retrieved spectrum for increasing values of interaction, observing an error below 10 % for [gN at ] ∼ 0.2.As compared to the previously discussed sources of error, we observe in Fig. 12 that the presence of interactions distort the entire spectrum, and not only the largest harmonics.Expressed in experimental parameters, for a gas of N at ∼ 10 4 atoms of 84 Sr, this value of interaction corresponds to a s = r 0 [gN at ]/(4πN at [r 0 ]) ∼ 10 −12 m, which is three orders of magnitude shorter than the scattering length of the experiment in Ref. [33].Using magnetic atoms such as Li, one could further rely on Feshbach resonances to engineer the required small values of scattering length [127].

IV. OUTLOOK
In this work, we have shown that the HHG emission yield can be simulated in current analog experiments using atomic clouds.We have characterized the response of the simulator to elliptic potentials and short pulses, as well us the main sources of errors for experimental implementations.Based on our calculations, we see that simulating attosecond dynamics by means of HHG is indeed possible with current analog simulators.They also offer a richer tunability on both the simulated incoming fields and the effective nuclear potential, where spatially modulated light can engineer potentials [128,129] that mimic the long-range coulomb attraction that ionized electrons experience during their excursion time.The flexibility of this approach offers several opportunities for the simulation of utltrafast processes: • As non-classical light sources become intense enough to drive the process of HHG [130], it reveals the limitations of semi-classical approximations [131], such as simple man's model that treats the driving field purely classically.This implies new challenges for semi-classical methods to accurately predict the HHG spectrum [130] or the harmonic coherence properties [132] for such cases.Interestingly, even if photons do not have a counterpart in this simulator, the response of the system to nonclassical light (e.g.squeezed pulses) can be approximated by an effective photon-statistics force [133] that can potentially be included in the simulator thanks to the high tunability of the effective incoming field.
• Further extensions of this platform would focus on studying multielectronic processes.For example, using shaken potentials in the ATI regime and the refined control provided by optical tweezers, each electron can be mapped to an individual fermionic atom in the trap.Extended atom-atom repulsive interactions provided by e.g., magnetic dipoles [134], Rydberg dressing [135,136] or mediating particles [137,138], allow one to directly tailor Coulomb corrections for a tunable degree of electronic repulsion.Going beyond the single-active electron approximation for HHG, this offers the opportunity to simulate processes such as NSDI, where the electron-electron correlation plays a pivotal role and needs to be included in the theoretical models to satisfactorily describe the experimental observations [139][140][141][142][143].
• In this direction, the creation of tweezer arrays of controllable geometries can also assist with the understanding of rescattering processes that take place in the crystalline structure of materials [144,145].In this context, conventional theoretical analyses of HHG processes in solid-state media result in quasicontinuum harmonic spectra, lacking the clear presence of harmonic peaks that would be expected under ideal conditions.To make these harmonics distinctly visible, exceedingly small dephasing times (typically on the order of ∼ 1 fs) are numerically employed [146], which may conflict with experimental observations [147,148] where the observed dephasing times are in the range of 15-50 fs.While recent theoretical studies have suggested that achieving numerical convergence is crucial to obtain well-resolved harmonic spectra without requiring extremely low dephasing times [149], others have tried to explain this phenomenon by considering the effect of propagation [150], or destructive interference between spatially long trajectories caused by the spatial distribution of laser fields [151].In this regard, using a cleaner environment compared to real materials can offer a higher temporal and spatial resolution that can help discern the relevance of these interference effects.Furthermore, it can also be utilized to study the connection between the delocalized recombination of electrons [144,145,152] and the reduced dependence of the harmonic yield with the driving field's ellipticity [55].
Overall, the tunability and accessibility of these experiments can be of significant value as a complementary tool to benchmark recently accessed experimental regimes, and stimulate the development of novel numerical techniques and theoretical models, helping to reach a more profound understanding of the electronic response of atoms and materials to ultrafast and intense laser fields.
where the operators H x (t) = V (x) + F 0 cos(ωt)x and H p = p 2 /(2m) are diagonal in real and momentum space respectively.This simplifies their exponentiation by appropriately Fourier transforming the state from the real to the momentum space, where either of the trotterized evolutions is diagonal.Along the Figures of this work, we considered τ 0 ω = 0.01.In Fig. C1, we show the velocity components of the atomic cloud (corresponding to the squared wavefunction in momentum space) at a given instant of the driving.This curve can be measured in the experiment through TOFM, where the momentum is mapped to different positions of the atomic cloud after a ballistic expansion of the gas.In the inset, we show that the relative error of the extracted dipole acceleration scales polynomially with the noise on the atomic density detection, providing an error smaller than 1% for ∆|ψ| 2 ∼ 10 −4 .This indicates a desired accuracy for experimental fluorescence measurements.

FIG. 1 .
FIG. 1. Schematic representation of HHG in atoms (a), and the proposed analog simulator (b).In (a) an ultrafast incoming field (brown) accelerates an electron (green), that is originally trapped by the nuclear Coulomb potential of the atom (red).The resulting oscillation of the charge emits radiation of characteristic harmonic frequencies (inset in blue).The same emission yield can be simulated on the simulator (b), where the characteristic frequencies are retrieved through absorption images of an atomic gas (green) that is trapped by a laser potential (red), and addressed by an external magnetic gradient that is tuned over time (brown).

FIG. 2 .
FIG. 2. Schematic representation of the different ionizationmechanisms at the instant of maximum tilt for different values of field strength.In (a) γ ≫ 1 such that it slightly perturbs the atomic potential, potentially leading to multiphoton ionization mechanisms.In (b) γ ≲ 1, such that the field is strong enough to allow for tunneling ionization.In (c), the electric field has reached a critical value that makes the maximum of the barrier to coincide with the energy level of the electron, leading to an over-the-barrier ionization mechanism.Note that in all these cases, we are working in the regime K0 ≫ 1.
FIG. 3. (a) Comparison between the Coulomb potential, −1/|x| (black line), the soft-core one, Vsc(x, r0) (red) and the Gaussian potential Vn(x, r0) (blue) for [r0] = 1.(b) Ionization energy associated to the Gaussian and soft-core potential (same colour code as before), for different values of r0.Black dashed line follows the Hydrogen-like condition, [Ip] = 0.5.Inset indicates the average ground-state width.See Appendix C for details on the numerical calculation.

FIG. 5 .
FIG. 5. (a)Conversion efficiency, η, of the 5th harmonic for the Gaussian potential in Fig.4, a six-cycle pulse (nc = 6), and different values of Keldysh (γ) and multiquantum parameter (K0).Red line 1 ○ follows the region where a higher efficiency is expected from Eq. (6), dashed lines 3 ○ and 4 ○ describe the HHG region (K0 ≥ 1 and γ ≤ 1, respectively). 5 ○ follows the frontier where the selected harmonic is below the cutoff frequency in Eq. 5.The dipole approximation 2 ○, associated to the upper bound in Eq.(7), lies above the represented range of values for K0.Arrows point towards the region of validity for each of these conditions.Circled red marker indicates the Gaussian configuration in Fig.4, while the crossed marker indicates the point K0 = 6.31 and γ = 0.81, associated to [F0] = 0.171 and [ω] = 0.108.The simulated emission spectrum of the latter configuration is illustrated in panel (b) for incoming pulses with a number of cycles running from nc = 6 (red) to nc = 1 (green).To improve visibility, an arbitrary vertical shift between the curves has been applied.

FIG. 6 .
FIG.6.Conversion efficiency, η, of the 5th harmonic for the Gaussian potential in Fig.4, a single-cycle pulse, nc = 1, and different values of Keldysh (γ) and multiquantum parameter (K0).Red line 1 ○ follows the region where a higher efficiency is expected from Eq. (6), dashed lines 3 ○ and 4 ○ describe the HHG region (K0 ≥ 1 and γ ≤ 1, respectively).5○ follows the frontier where the selected harmonic is below the cutoff frequency in Eq. 5. Arrows point towards the region of validity for each of these conditions.The crossed marker indicates the point K0 = 6.31 and γ = 0.81, whose emission spectrum is depicted in green in Fig.5(b).

FIG. 7 .
FIG.6.Conversion efficiency, η, of the 5th harmonic for the Gaussian potential in Fig.4, a single-cycle pulse, nc = 1, and different values of Keldysh (γ) and multiquantum parameter (K0).Red line 1 ○ follows the region where a higher efficiency is expected from Eq. (6), dashed lines 3 ○ and 4 ○ describe the HHG region (K0 ≥ 1 and γ ≤ 1, respectively).5○ follows the frontier where the selected harmonic is below the cutoff frequency in Eq. 5. Arrows point towards the region of validity for each of these conditions.The crossed marker indicates the point K0 = 6.31 and γ = 0.81, whose emission spectrum is depicted in green in Fig.5(b).

FIG. 10 .
FIG.10.mission yield extracted from the discrete Fourier transform where the position of the temporal intervals randomly deviate from a uniform distribution with a standard deviation, ω∆t ∼ 10 −3 (red dotted line), as compared to the perfect spectrum (blue line).Orange dots indicate the positions of local maxima before the cutoff energy.In the inset, we show the relative error of these local maxima for increasing values of standard deviation, ∆t, where the choice illustrated before is highlighted with a red cross.The rest of parameters coincide with Fig.5(b).

FIG. 11 .
FIG.11.Emission yield extracted from the discrete Fourier transform of the temporal dipole acceleration, whose values are randomly distorted with standard deviation, ∆da(t)/ da = 0.02 (red dotted line), as compared to the perfect spectrum (blue line).Orange dots indicate the positions of local maxima before the cutoff energy.In the inset, we show the relative error of these local maxima for increasing values of standard deviation, ∆da(t), where the choice illustrated before is highlighted with a red cross.The rest of parameters coincide with Fig.5(b).

FIG. 12 .
FIG.12.Emission yield associated to an atomic cloud with interatomic interaction strength, [gNat] = 0.2 (red dotted line), as compared to the noninteracting case (blue line).Orange dots indicate the positions of local maxima before the cutoff energy.In the inset, we show the relative error of these local maxima for increasing values of interacting strength, where the choice illustrated before is highlighted with a red cross.The rest of parameters coincide with Fig.5(b).
FIG. C1.Velocity components of the atomic cloud, |ψ(v)| 2 after half of an incoming 6-cycle pulse is applied for the same parameters as in Fig. 5(b).The inset shows the relative error on the reconstructed dipole acceleration, da(t = T /2) as one introduces Gaussian noise of standard deviation ∆|ψ| 2 on the velocity components of the main Figure.The markers indicate the average over 10 random realizations.