Spacetime resolved simulation of femtosecond nonlinear lightmatter interactions using a holistic quantum atomic model : Application to nearthreshold harmonics
Abstract
We introduce a new computational approach for femtosecond pulse propagation in the transparency region of gases that permits full resolution in three space dimensions plus time while fully incorporating quantum coherent effects such as highharmonic generation and strongfield ionization in a holistic fashion. This is achieved by utilizing a onedimensional model atom with a deltafunction potential which allows for a closedform solution for the nonlinear optical response due to groundstate to continuum transitions. It sidesteps evaluation of the wave function, and offers more than one hundredfold reduction in computation time in comparison to direct solution of the atomic Schrödinger equation. To illustrate the capability of our new computational approach, we apply it to the example of nearthreshold harmonic generation in Xenon, and we also present a qualitative comparison between our model and results from an inhouse experiment on extreme ultraviolet generation in a femtosecond enhancement cavity.
College of Optical Sciences, University of Arizona, Tucson, AZ 85721, U.S.A.
Department of Physics, Constantine the Philosopher University, Nitra, Slovakia
kolesik@acms.arizona.edu
(320.0320) Ultrafast optics; (320.7110) Ultrafast Nonlinear Optics; (190.4610) Multiharmonic generation; (190.5940) Selfaction effects; (320.2250) Femtosecond phenomena;
References
 [1] A. Couairon and A. Mysyrowicz, “Femtosecond filamentation in transparent media,” Phys. Rep. 441, 47–189 (2007).
 [2] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, and P. B. Corkum, “Theory of high harmonics generation by lowfrequency laser fields,” Phys. Rev. A 49, 2117–2132 (1994).
 [3] L. Bergé, S. Skupin, R. Nuter, J. Kasparian, and J. Wolf, “Ultrashort filaments of light in weakly ionized optically transparent media,” Rep. Prog. Phys. 70, 1633–1713 (2007).
 [4] M. Gaarde, J. Tate, and K. Schafer, “Macroscopic aspects of attosecond pulse generation,” J. Phys. B: At. Mol. Opt. Phys. 41, 132001 (2008).
 [5] M. Nurhuda, A. Suda, and K. Midorikawa, “Generalization of the Kerr effect for high intensity, ultrashort laser pulses,” New J. Phys. 10, 053006 (2008).
 [6] E. A. Volkova, A. M. Popov, and O. V. Tikhonova, “Nonlinear polarization response of an atomic gas medium in the field of a highintensity femtosecond laser pulse,” JETP Lett. 94, 519–524 (2011).
 [7] I. Christov, “Propagation of ultrashort pulses in gaseous medium:breakdown of the quasistatic approximation,” Opt. Express 6, 34–39 (2000).
 [8] E. Lorin, S. Chelkowski, and A. D. Bandrauk, “The WASP model: A MicroMacro system of WaveSchrödingerPlasma equations for filamentation,” Commun. Comput. Phys. 9, 406–440 (2011).
 [9] Q. Su and J. Eberly, “Model atom for multiphoton physics,” Phys. Rev. A 44, 5997–6008 (1991).
 [10] G. P. Arrighini and M. Gavarini, ‘‘Ionization of a model atom by strong and superstrong electric fields,” Lett. Nuovo Cimento 33, 353–358 (1982).
 [11] W. Elberfeld and M. Kleber, “Tunneling from an ultrathin quantum well in a strong electrostatic field: A comparison of different methods,” Z. Phys. B: Cond. Matt. 73, 23–32 (1988).
 [12] R. M. Cavalcanti, P. Giacconi, and R. Soldati, “Decay in a uniform field: an exactly solvable model,” J. Phys. A: Math. Gen. 36, 12065–12080 (2003).
 [13] H. Uncu, H. Erkol, E. Demiralp, and H. Beker, “Solutions of the Schrödinger equation for dirac delta decorated linear potential,” C. Eur. J. Phys. 3, 303–323 (2005).
 [14] G. Álvarez and B. Sundaram, “Perturbation theory for the stark effect in a double quantum well,” J. Phys. A: Math. Gen. 37, 9735–9748 (2004).
 [15] V. M. Villalba and L. A. GonzálezDíaz, ‘‘Particle resonance in the dirac equation in the presence of a delta interaction and a perturbative hyperbolic potential,” Eur. Phys. J. C 61, 519–525 (2009).
 [16] S. Geltman, “Ionisation dynamics of a model atom in an electrostatic field,” J. Phys. B: At. Mol. Phys. 11, 3323–3337 (1978).
 [17] G. V. Dunne and C. S. Gauthier, “Simple soluble molecular ionization model,” Phys. Rev. A 69, 053409 (2004).
 [18] Q. Su, B. P. Irving, C. W. Johnson, and J. H. Eberly, “Stabilization of a onedimensional shortrange model atom in intense laser fields,” J. Phys. B: At. Mol. Opt. 29, 5755–5764 (1996).
 [19] T. Dziubak and J. Matulewski, “Stabilization of onedimensional softcore and singular model atoms,” Eur. Phys. J. D 59, 321–327 (2010).
 [20] A. Sanpera, Q. Su, and L. RosoFranco, “Ionization suppression in a veryshortrange potential,” Phys. Rev. A 47, 2312–2318 (1993).
 [21] S. Geltman, “Shortpulse modelatom studies of ionization in intense laser fields,” J. Phys. B  At. Mol. Opt. 27, 1497–1514 (1994).
 [22] Z. X. Zhao, B. D. Esry, and C. D. Lin, “Boundaryfree scaling calculation of the timedependent Schrödinger equation for laseratom interactions,” Phys. Rev. A 65, 023402 (2002).
 [23] A. Teleki, E. M. Wright, and M. Kolesik, “Microscopic model for the higherorder nonlinearity in optical filaments,” Phys. Rev. A 82, 065801 (2010).
 [24] J. M. Brown, E. M. Wright, J. V. Moloney, and M. Kolesik, “On the relative roles of higherorder nonlinearity and ionization in ultrafast lightmatter interactions,” Opt. Lett. 37, 1604–1606 (2012).
 [25] J. Brown, A. Lotti, A. Teleki, and M. Kolesik, “Exactly solvable model for nonlinear lightmatter interaction in an arbitrary timedependent field,” Phys. Rev. A 84, 063424 (2011).
 [26] J. Lee, D. R. Carlson, and R. J. Jones, “Optimizing intracavity high harmonic generation for xuv fs frequency combs,” Opt. Express 19, 23315–23326 (2011).
 [27] M. Kolesik and J. V. Moloney, “Nonlinear optical pulse propagation simulation: From Maxwell’s to unidirectional equations,” Phys. Rev. E 70, 036604 (2004).
 [28] A. BideauMehu, Y. Guern, R. Abjean, and A. JohanninGilles, “Measurement of refractive indices of neon, argon, krypton and xenon in the 253.7—140.4 nm wavelength range. Dispersion relations and estimated oscillator strengths of the resonance lines,” J. Quant. Spect. Rad. Transfer 25, 395–402 (1981).
 [29] P. B. Corkum, “Plasma perspective on strong field multiphoton ionization,” Phys. Rev. Lett 71, 1994–1997 (1993).
 [30] K. J. Schafer, B. Yang, L. F. DiMauro, and K. C. Kulander, “Above threshold ionization beyond the high harmonic cutoff,” Phys. Rev. Lett. 70, 1599–1602 (1993).
 [31] F. Schapper, M. Holler, T. Auguste, A. Zaïr, M. Weger, P. Salières, L. Gallmann, and U. Keller, “Spatial fingerprint of quantum path interferences in high order harmonic generation,” Opt. Express 18, 2987–2994 (2010).
 [32] S. M. Teichmann, D. R. Austin, P. Bates, S. Cousin, A. Grün, M. Clerici, A. Lotti, D. Faccio, P. DiTrapani, A. Couairon, and J. Biegert, “Trajectory interferences in a semiinfinite gas cell” Laser Phys. Lett. 9, 207–211 (2012).
 [33] A. L’Huillier, K.J. Schafer, and K.C. Kulander, “Highorder harmonic generation in Xenon at 1064 nm: The role of phase matching,” Phys. Rev. Lett. 66, 2200–2203 (1991).
 [34] P. Salières, A. L’Huillier, and M. Lewenstein, “Coherence control of highorder harmonics,” Phys. Rev. Lett. 74 3776–3779 (1995).
 [35] E. Constant, D. Garzella, P. Breger, E. MÃ©vel, C. Dorrer, C. Le Blanc, F. Sali, and P. Agostini, “Optimizing high harmonic generation in absorbing gases: Model and experiment,” Phys. Rev. Lett. 82, 1668–1671 (1999).
 [36] R. J. Jones, K. D. Moll, M. J. Thorpe, and J. Ye, “Phasecoherent frequency combs in the vacuum ultraviolet via highharmonic generation inside a femtosecond enhancement cavity,” Phys. Rev. Lett. 94, 193201 (2005).
 [37] C. Gohle, T. Udem, M. Herrmann, J. Rauschenberger, R. Holzwarth, H. A. Schuessler, F. Krausz, and T. W. Hansch, “A frequency comb in the extreme ultraviolet,” Nature 436, 234–237 (2005).
 [38] D. R. Carlson, J. Lee, J. Mongelli, E. M. Wright, and R. J. Jones, ‘‘Intracavity ionization and pulse formation in femtosecond enhancement cavities,” Optics Letters 36, 2991–2993 (2011).
Appendix A Introduction
The goal of this paper is to present a new computational approach for nonlinear lightmatter interactions in atomic and molecular gases for incident femtosecond pulses whose spectra lie in the transparency region, the approach opening the door to fully resolved propagation in three spatial dimensions plus time, or (3+1) dimensions. For sufficiently low input intensities when ionization of the constituent atoms is negligible and the atoms remain dominantly in their ground state, suitable spacetime resolved computational approaches already exist and incorporate the linear dispersive properties of the medium, both refractiveindex and absorption, along with a model for the bound state electronic nonlinearity, typically a Kerrtype nonlinearity possibly with saturation, and also a timedelayed nonlinear response to model Raman effects for molecules [1]. For higher intensities when atomic ionization becomes relevant a number of effects turn on that involve ground state to continuum state transitions. These transitions can greatly modify field propagation via nonlinear generation of freed electrons and associated nonlinear losses and refraction, and high harmonic generation (HHG). Typically the nonlinear optical response arising from ground state to continuum transitions are dealt with in a piecemeal fashion. For example, the polarization source term for high harmonics is often calculated using the strongfield approximation [2] whereas the rate of ionization is treated in a separate calculation using Keldysh theory or one of the subsequent variants such as ADK theory (see [3] for a useful compendium of different ionizationrate formulas). Such an approach has met with considerable success to date for simulations of HHG (ref. [4] has a readable yet detailed description of the method) and also light filament propagation in gases [1], but it would clearly be advantageous for the advancement of these fields if the HHG and all ionization related effects could be modeled in a holistic manner so that as parameters are varied the relative roles of highharmonic generation and ionization related effects are microscopically consistent without having to tune the parameters of separate models. An obvious candidate for such a quantum model is to solve the threedimensional (3D) atomic Schrödinger equation [5, 6] to obtain the quantum averaged dipole moment at each point in space, and use this as a source in Maxwell’s equation for the field propagation. However, due to the restrictive computing requirements needed to solve the 3D Schrödinger equation it seems clear that this will not be a viable option for the foreseeable future.
The quantum model we employ in the current work is onedimensional, and it is perhaps the simplest system which has both the energetic continuum and a bound state — the minimal set of ingredients crucial to describe lightmatter interactions at femtosecond time scales and ionizationlevel intensities. However, the model has the virtue that it can be implemented to efficiently simulate spacetime resolved field propagation in (3+1) dimensions. To the best of our knowledge Christov [7] was the first to perform simulations of femtosecond pulse propagation in one space dimension plus time, or (1+1) dimensions, that incorporated the medium response via direct numerical solution of a timedependent 1D atomic model. In a series of recent papers Bandrauk and coworkers [8] have presented simulations of 1D quantum models coupled to field propagation in both (1+1) and (3+1) dimensions, and addressed a range of topics of current interest including HHG and filamentation in gases. Onedimensional atomic models have played an important role as model systems for elucidating the physics underlying nonlinear lightmatter interactions. For example, Su and Eberly [9] introduced the quasiCoulomb potential in 1D as a model for multiphoton physics. Here we use the potential atomic model. This model has one free parameter, the ionization potential, and its spectrum consists of one bound state, the ground state, plus the continuum, so it is an ideal model system in which to explore the role of ground state to continuum transitions. Indeed, this system has been studied and applied in various contexts for many years. The attraction of the model as a toy system to investigate various aspects of dynamics in quantum systems stems from the fact that many quantities can be obtained exactly [10, 11]. Its generalizations to higher dimensions [12], for more complex potentials [13, 14], and other dynamic equations [15] also exist. Besides its utility in the AMO and strong field physics [16, 17, 18, 19], the 1D potential has also been applied in semiconductors [11, 13], and used as a test system for numerical and approximate methods [20, 21, 22]. We have applied the potential model to explore the relative roles of the higherorder Kerr effect and plasma induced defocusing in filament propagation [23, 24]. Furthermore, it has recently been shown by two of the present authors (MK and JMB) that the potential model yields a closed form solution for the quantum averaged nonlinear current for this model atom in an arbitrary timedependent light field [25]. (The code computing the nonlinear current is available as opensource upon request to the corresponding author (MK).) This exact solution side steps the need for direct evaluation of the timedependent atomic wave function and therefore represents an enormous savings in computational effort and time. For example, one does not need to worry about griding issues for the solution of the 1D wave equation and the spurious effects of numerical boundary conditions. Our new approach is computationally efficient since it employs the closed form for the nonlinear current as a source in the (3+1) field propagation equation, and we find gains of more than one hundred in computation time in comparison to the approach using integration of the Schrödinger equation.
The remainder of this paper is organized as follows: In Sec. 2 we describe our computational approach and discuss the model equations. In Sec. 3 we include numerical simulations, and for this initial paper we utilize nearthreshold HHG in Xenon as an illustrative example. We chose nearthreshold HHG, with concomitant lowerorders in the harmonic spectrum, to showcase the capability of the model to cope with an extreme spectral bandwidth in a completely unified treatment. Unlike the strongfield approximation, our approach does not rely upon identification of the quantum paths that dominate the HHG spectrum, and is valid also for very low frequencies. The first part of Sec. 3 describes simulations of HHG in a Xenon gas jet to demonstrate that our model can deal with different incident focusing conditions and hence phasematching regimes, and that it reproduces expected qualitative features. In this respect we do not claim these results are new, particularly since corresponding simulations of HHG in (3+1) already exist that employ the strongfield approximation, but rather are intended to highlight the capabilities of the approach. Sec. 3 also includes what we believe are new results to model nearthreshold HHG in a femtosecond enhancement cavity (fsEC). Here a new feature is that the HHG occurs in the presence of the residual plasma that accumulates within the Xenon gas jet due to the circulating pulse, and we compare the qualitative features of the harmonic spectrum with that of the inhouse experiment that has recently been used to produce record power levels of extreme ultraviolet (XUV) [26].
Appendix B Computational approach and model
To contrast the model of this paper to the current state of the art in computer simulation of highharmonic generation and/or femtosecond filamentation, we briefly recall the main features of the standard approach. The propagation of the optical and highharmonic fields are described in oneway (i.e. directional) propagation equations, into which the various medium responses are coupled. For this purpose, we use the Unidirectional Pulse Propagation Equation (UPPE) solver. In the spectral domain the UPPE for propagation along the zaxis takes the general form [27]
(1) 
where the zcomponent of the optical wavevector is given by
(2) 
The UPPE describes the evolution of the spatially resolved spectral amplitudes of the scalar electric field , where is the frequency, and denotes transverse wavenumbers, so that the spacetime dependent electric field may always be reconstructed using a 3D Fourier (or Hankel) transform over . Note that the assumption of a linearly polarized field is consistent with our adoption of a 1D atomic model. Nonlinear effects enter the propagation equations through either the polarization () or currentdensity () terms. These are both evaluated as functionals of the electric field history at each spatial point, and then Fouriertransformed to the spectral representation, e.g. , that appears in the above propagation equation. The linear wavelengthdependent properties of a medium are represented through the susceptibility.
In the standard model, there are a number of independent contributions both in the polarization and in the current density. They include the optical Kerr effect, ionization in strong fields and a freedelectron density evolution equation, separately modeled losses due to ionization, avalanche ionization, defocusing effects of freed electrons, and losses due to freed electrons described in terms of an effective Drudeplasma model [1].
Our computational approach for simulating femtosecond nonlinear lightmatter interactions reflects the fact that the above mentioned effects are all manifestation of the single electronic system response to the strong optical field. The goal is to achieve a unified description, and reduce multiple independent model parameters. In general, our model involves three components for describing an atomic gas:

A description of the linear dispersive properties of the gas via a complexvalued frequency dependent susceptibility .

The 1D potential quantum model for modeling the nonlinear current () due to ground state to continuum transitions. This will incorporate the effects of ionization, HHG, and ionization induced absorption and refraction in a holistic fashion.

A Kerrlike nonlinearity and associated nonlinear polarization to capture the nonlinear optical response due to the ground state to excited bound state transitions. This will contribute processes such as fourwave mixing, selfphase modulation, and selffocusing. (For a molecular gas we may also add a timedelayed nonlinear response to capture Raman effects)
For this initial exposition we shall give a description of each of the three components above in the subsections below as appropriate to our illustrative example of nearthreshold HHG in Xenon gas. The three medium components are coupled into the fieldpropagation equation which encompasses all frequency components from the fundamental to high harmonics. While we use the Unidirectional Pulse Propagation Equation, we emphasize that this approach may be implemented with any pulsepropagation simulator that resolves the carrier wave of the optical field.
Linear dispersive properties
Inclusion of the linear dispersive properties of the medium over the full span of harmonics is important to properly incorporate the phasematching and reabsorption that affect the spatial evolution of the generated harmonics. In our spectral solver this does not present a problem as long as suitable data is available for the index of refraction and absorption properties in a sufficiently wide spectral range. The medium is characterized in terms of the linear complex susceptibility , and the pulse propagation method utilizes this information in the propagation constant in (2) at each frequency or wavelength resolved by the numerical simulation.
Throughout this paper we use Xenon as a representative atomic medium to showcase the modeling capabilities, and Fig. 1 shows the real (solid dark line) and imaginary (red line) parts of the susceptibility as functions of the angular frequency of the electromagnetic field. To create the data set displayed in Fig. 1 we have combined Xenon data downloaded from http://henke.lbl.gov/optical_constants/ with the refractiveindex parametrization obtained from Ref.[28] to create a tabulated representation of the linear susceptibility over a range of frequencies that spans the fundamental at 800 nm wavelength up to its harmonic. Note that with the real and imaginary parts of the susceptibility function coming from separate sources, and with no absorption data available for photon energies below 10 eV, our parametrization does not satisfy the KramersKronig relations. To our best knowledge a more complete data set is currently not available. Nevertheless, while not a perfect model of the Xenon gas for all frequencies this data set employed serves a an illustration example in which all electromagnetic frequencies are treated on the same footing by the propagation solver, which operates on the full electric field. Generally the accuracy of the linear propagation will be limited by the quality of the available data for the linear susceptibility.
Quantum atomic model for the nonlinear response
As alluded to in the introduction our computational approach employs a 1D quantum model with a shortrange potential for the nonlinear current due to ground state to continuum transitions. More specifically the corresponding timedependent Schrödinger equation in atomic units takes the form
(3) 
Here and are the space and time variables in atomic units, the strength of the attractive potential is governed by the parameter , with the ionization potential of the ground state equal to , and denotes the timedependent external field in atomic units applied to the atom.
The above Schrödinger equation would have to be solved numerically at each spatial point resolved by the optical pulse simulator. Such an approach, although feasible for this specific quantum system, is very demanding on computing resources. Fortunately, we do not need the wave function because it is only the timedependent dipole moment or current that enters as an input in the UPPE (1), and we have recently shown how an exact, closed form solution can be obtained for both for an arbitrary external field [25]. More specifically, the solution provides a means of evaluating the expectation value for the current in atomic units
(4) 
without having to go through solving Eq. 3 for the wave function . Rather, the current required for the UPPE is calculated directly, completely bypassing the wave function, which amounts to a tremendous savings in computing time while retaining the full coherent quantum dynamics of the atomic system as is pertinent to the nonlinear field propagation. Moreover, it is possible to exactly eliminate the current component which is linear in the driving field strength [25]. The quantum model can therefore be restricted to contribute only to the nonlinear response, which is advantageous when it is combined with the linear dispersive medium. For the sake of completeness the formulas for the nonlinear quantum current are summarized in the Appendix. The reader is referred to [25] for details concerning the practical numerical implementation. For those interested in using this model in their simulations, the corresponding software is available as an opensource upon request to the corresponding author (MK).
The integration with the pulse propagator solver is in principle the same as for any other nonlinear medium response. Having calculated the evolution of the optical field up to a given propagation distance, the history of the electric field at a given point in space is converted to the external field in atomic units. This drives the quantum system, and the induced current is computed as outlined in the Appendix. The nonlinear current is next converted from the atomic units, and multiplied by the number density atoms at the point in space. The resulting macroscopic current density is included in the righthandside of the UPPE equation. We remark that since we only incorporate the nonlinear current from the quantum model into the UPPE, we do not doublecount by erroneously including linear properties from the 1D atomic model. Also note that it is not an option to retain the linear part of the model’s response instead of introduced in the previous subsection. This is because in a real medium originates in virtual transitions among a large number of states, and it would be difficult to model this from first principles. The linear susceptibility arising from our 1D atomic model is far too simplistic to capture the chromatic properties of a gas to any realistic degree.
Next we would like to illustrate that our 1D potential atomic model displays features for HHG that are qualitatively similar to those obtained using the strong field approximation applied to the more exact 3D Hydrogenlike atomic model. To this end Fig. 2 shows an example of a harmonic spectrum of the nonlinear current induced in the 1D model atom for an ionization potential of 12 eV, characteristic of Xenon, and a tencycle pulse of center wavelength nm and peak intensityW/m. This harmonic spectrum exhibits the characteristic highharmonic generation plateau with a high frequency cutoff. The cutoff predicted by the formula from the standard HHG model [29, 2, 30] is indicated by the arrow in Fig. 2 for our parameters. Furthermore, the harmonic spectrum includes the harmonic which occurs near the ionization potential up to around the harmonic, and thus constitutes an example of nearthreshold HHG where the generated harmonics straddle the ionization energy. What is shown in Fig. 2 is the spectrum of the nonlinear polarization that appears in the Maxwell equations and not the spectrum of the radiation actually generated. Importantly the harmonic spectrum shown is exact, within the context of our 1D model, over the whole frequency range and in particular at low frequencies.
Finally we point out some pros and cons of our 1D model versus the strongfield approximation. The strong field approximation is often utilized in numerical atomic simulations of HHG, and describes the 3D atomic system in terms of a single bound state, the ground state, and a continuum of free electron continuum states [2]. In this sense it addresses a similar spectrum of electronic states as our 1D model, albeit in 3D. A distinct advantage of the strong field approximation is that it can in principle be used for an elliptically polarized driving field, whereas our 1D model is restricted to linear polarization. On the other hand, ours is an exact solution of a welldefined system, and as such it is valid throughout the whole frequency spectrum. Unlike the strong field approximation, it accounts for “all electron trajectories” not only those that give rise to the harmonic radiation. In particular, the lowfrequency components of the nonlinear current response affect the propagation of the driver pulse through ionization, defocusing by freed electrons, ionization losses, and the nonlinear focusing.
The nonlinear Kerr effect
The potential atomic model only incorporates the nonlinear contribution of ground state to continuum transitions. To capture the nonlinear contribution of virtual ground state to bound state electronic transitions we include a term representing the nonlinear Kerr effect of the form
(5) 
where is the nonlinear index and is the linear refractive index. This effect enters the propagation equation via the nonlinear polarization term, and gives rise to selffocusing and a cascade of lowerharmonic radiation.
Appendix C Illustrations involving nearthreshold highharmonic generation
HHG in a Xenon gas jet
The first illustrative example involves nearthreshold HHG in a Xenon gas jet. The incident pulse in each case is Gaussian in space and time with pulse duration fs and center wavelength of nm, and focused spot size . We model the Xenon gas jet as a vertical column which is tapered along the propagation or zaxis according to the jet pressure profile indicated by the dotted line in Fig. 3. This pressure or density profile has a constant region of length m between m and tapers off on each side of this region. For the examples in this paper we used 20 Torr for the maximal pressure in the gas jet. Furthermore we consider the two distinct focusing cases where the input beam is focused at m close to the entrance to the gas jet, and also the case that the input beam is focused at m close to the exit of the gas jet (We use this descriptive nomenclature of exit and entrance to the gas jet with the caveat that HHG can arise before and after the exit and entrance due to the tapers in the gas density). The variation of the onaxis intensity along the zaxis is illustrated in Fig. 3 for the cases of a) a narrow beam with m, and b) a wide beam with m.
In what follows we shall present examples of fully resolved simulations of pulse propagation and harmonic generation. As a means to present our results we have chosen to plot the angularly resolved spectra for selected harmonics as they propagate through the gas jet. Such spectra have previously been utilized in both numerical simulations and were measured experimentally [31, 32]. For example, with reference to Figs. 4,6,7 each panel shows a color encoded plot of the spectral intensity (red being maximum and blue minimum) as a function of the transverse wavenumber representing angular spread along the horizontal axis, the center corresponding to zero propagation angle, and frequency along the vertical axis for the harmonic indicated. In particular, the frequency axis in each panel is centered on the harmonic order and has a spread of one half of the harmonic spacing. Each set of results is arranged into a array of panels, the panels are arranged vertically according to the propagation distance into the gas jet, the top panel being the entrance and the bottom panel the exit of the gas jet, and the panels are arranged horizontally from left to right for harmonic orders 9,11 and 13.
Tighter focus case
For the first example we chose a focused spot size of m giving a Rayleigh range of m for the input field. In this case and we expect that the focusing of the beam will impact phasematching of the HHG via the Gouy phaseshift incurred by the Gaussian beam. The calculated angularly resolved spectra are shown in Fig. 4 for (a) focusing at the exit, and (b) focusing at the entrance of the gas jet, and we note that there is a marked difference between the angularly resolved spectra for the two cases. This is consistent with general expectations based on phasematching [33, 34, 35]. For example, onaxis phasematching requires cancellation between the accumulated Gouy phaseshift and the accumulated atomic phaseshift due to the quantum path taken by the electron. The Gouy phaseshift is a positive quantity whereas the quantum phase is proportional to the zderivative of the incident beam intensity due to the dependence of the quantum phase on the electron quiver energy, which is positive for case (a) with the beam focused at the exit, and negative for the case (b) with the beam focused at the entrance to the gas jet. Thus onaxis phase matching is not a possibility for the case (a) with the input beam focused at the exit, whereas phasematching is a possibility in case (b) with focusing at the entrance to the gas jet. The difference in the angularly resolved spectra evident for cases (a) and (b) in Fig. 4 may thus be traced to expected differences in phasematching conditions [34].
If we further examine the angularly resolved spectra at the exit of the gas jet, shown as the lower set of panels in Fig. 4 for cases (a) and (b) we observe further consequences of the different phasematching conditions: For case (a) for which strict onaxis phasematching is not possible we see that the offaxis HHG emission is more significant with respect to the onaxis emissions for each harmonic order than for case (b) where onaxis phasematching is possible. This is particularly pronounced for the harmonic where case (a) is dominated by offaxis emission whereas case (b) is dominantly onaxis emission. So the results of our simulations are compatible with expectations based on onaxis phasematching [33, 34].
To continue our discussion Fig. 5 shows the evolution of the conversion efficiency for a variety of marked harmonics during the propagation through the gas jet. (The conversion efficiency is obtained by integrating the angularly resolved spectrum over the transverse wavenumber plane for a given harmonic order.) It compares the cases of (a) focus at the exit (left), and focus at the entrance to the gas jet (right), and shows that the efficiency is higher in the latter case, as is generally accepted to be the case [34, 4]. Furthermore, Fig. 4 reveals that the evolution of the harmonics during the pulse propagation through the jet is distinctly nonmonotonic, that is, different angular portions of the harmonic fields dominate at different propagation distances. To illustrate that our simulation capability can capture the complex interplay between HHG processes and phasematching considerations due to variation in the incident field and gas jet pressure, Fig. 6 plots the evolution of the harmonic conversion efficiency for several harmonics, the upper panels showing the evolution of the angularly resolved spectrum for the harmonic at three propagation distances. At the shortest and longest propagation distances note that the spectrum has its maximum onaxis, whereas at the onaxis and offaxis become comparable in the middle. It is interesting that the turnaround in the conversion efficiency for the 9th harmonic that occurs for a propagation distance of around m is therefore dominated by a reduction in the offaxis as opposed to the onaxis emission. This highlights the point that the spatial evolution of the various harmonics involves a complex interplay between onaxis and offaxis emissions, and it is not generally sufficient to assess the conversion efficiency based on the onaxis phasematching considerations.
Weaker focus case
For this case we chose a focused spot size of m giving a Rayleigh range of m for the input field. Now, and we expect that the focusing of the beam will not impact phasematching of the HHG in a major way. The calculated angularly resolved spectra for this case are shown in Fig. 7 for the cases of (a) focusing at the exit, and (b) focusing at the entrance of the gas jet, and we note that there is no major differences difference between the angularly resolved spectra for the two cases. There are small quantitative differences though, in particular in the visibility of the ring structure (cf. the 11th harmonic farfield patterns).
Figure 8 provides a comparison of the conversion efficiency evolution versus the propagation distance for the two focusing cases. This figure shows that while the optimal position of the beam focus is still at the entrance into the gas jet, as generally accepted [34], for weaker focusing of the input field the qualitative behavior is relatively insensitive to the position of the input beam focus with respect to the gas jet.
Finally we would like to point to one other feature from our simulations that correlates with previous studies: Inspection of the high resolution image for the angularly resolved spectra for the 11th harmonic in Fig. 7 reveals that the output spectra in the bottom row display ring structures. Such structures have previously been seen experimentally and appeared in propagation simulation in (3+1) dimensions employing the strongfield approximation [34, 31]. Moreover the appearance of these anglefrequency spectral features has been associated with the interference between the short and long path contributions to HHG. Recently, Teichman et al.[32] demonstrated, experimentally and numerically, that rings in the anglefrequency spectra can be attributed to the interferences in the singleatom response in the transverse crosssection profile, and that their appearance depends on particular conditions such as intensity and specific harmonics. It is satisfying that our computational approach can reveal such structures without resort to the strong field approximation. Moreover, inspection of nearfield crosssections (not shown) of the frequencyfiltered atomic response shows that complex farfield patterns correlate with the complex spatial transverse distribution of the HHG “source.” This suggests, in keeping with [32], that dynamics which controls the transverse properties of the atomic response play equally important role in macroscopic selection through phase matching.
HHG in a femtosecond enhancement cavity
As a final illustration, we show results for harmonic generation in a xenon gas jet placed inside a passive femtosecond enhancement cavity (fsEC). This experimental configuration has been demonstrated as a means to generate XUV frequency combs with the harmonic light being generated at each cavity round trip [36, 37]. In this geometry the ionized gas jet behaves as a nonlinear medium at the focus of a 6 m ring cavity. A pair of 15 cm radius of curvature focusing mirrors lead to an intracavity spot size of m within the jet. The fsEC has a 1% input coupler that allows pulse energy enhancement over 200 times relative to the incident pulse train leading to peak intensities in excess of W/cm when the gas jet is off. The harmonics are coupled out of the cavity using a thin sapphire plate, aligned at Brewster’s angle for the fundamental pulse, and resolved spatially by reflection from an XUV grating.
For initial conditions in these simulations, we have utilized a previously calculated temporal profile of the pulse just before it enters the gas jet. This pulse profile was obtained from a 1D steadystate calculation of the fundamental pulse building up inside the fsEC in the presence of the xenon gas and includes the effects of cavity dispersion. The high nonlinearity of ionization leads to chirped pulses circulating in the cavity and limits the achievable intensities in this geometry. Details of this model and its implications for HHG in a fsEC can be found in Ref. [38].
Yet another modification of the current model consists in including the surviving plasma in the jet based on the estimated levels calculated in [38] . Because the 20 ns cavity roundtrip time is too short for the plasma to decay entirely, the pulse propagating through the jet experiences defocusing due to electrons freed during the previous passes. These electrons have been liberated from their parent atom for a sufficiently long time, and have equilibrated into a true plasma (note that the freeelectron states included in the quantum model that drives the pulse propagator are of different nature: they did not have enough time to thermalize, and their interactions is mainly with the nearby parent ion). Therefore, the influence of the plasma remnant can be included within the linear chromatic properties of the gas. Here we also assume that the diffusion resulted in a nearly homogeneous spatial distribution of free electrons, and the corresponding modification of the refractive index is constant over the crosssection of the beam. Of course, the spectral nature of our pulse solver allows to endow this susceptibility contribution with the correct frequency dependence, .
Figure 9 shows simulated and experimental spectra of harmonics from 7th to 15th. While details of relative strength of harmonics 9th to 15th depend on the exact placement of the beam focus with respect to the center of the jet (compare the two panels on the left), and on the absorption included in our simulation (compare full and dashed lines), we see that harmonics 11 and 13 are the most pronounced, and that harmonic power starts to decrease at harmonics 15th (there are many more harmonics generated, but are not discernible on the linear scale of this figure). This is compatible with the spectrum recorded in the experiment. On the other hand, the simulated 7th harmonic seems to be too strong. We think that this is mainly due to the fact that the model susceptibility of the gas has no absorption at this wavelength (this is due to limited frequency extent of the available data, see the red line in Fig. 1). The effect of medium absorption is indeed clearly visible in the 9th harmonic where our model absorption data exhibit a maximum. These results make it thus evident that it is important to obtain as realistic as possible a data set for both the index of refraction and absorption of the gas.
We certainly do not expect that the simplistic onedimensional quantum model put forward here could reproduce experimental results in a onetoone manner. However, we see the above comparison as very encouraging — it shows that the model, especially when coupled with goodquality dispersion and absorption data for the gas, can serve as a practical tool to study the qualitative trends that govern the nonlinear interactions that span the frequency range from infrared to high harmonics.
Appendix D Conclusion
In conclusion, we have presented a new computational approach for femtosecond pulse propagation in the transparency region of gases that permits full resolution in three space dimensions plus time while fully incorporating quantum coherent effects such as highharmonic generation and strongfield ionization in a holistic fashion. The key innovation that makes this possible is the use of 1D atomic model with deltafunction potential to compute the nonlinear optical response due to ground state to continuum transitions, and the closed form solution for the nonlinear response leads to significant reductions in computing time compared to a direct solution of the atomic Schrödinger equation. In spite of the approximation of a 1D atomic model we contend that our computational approach will be of great utility as a research tool in a number of forefront areas. For example, in the field of femtosecond filamentation in gases an outstanding question is whether coherent lightmatter interactions can play a significant role. The current models, based upon individual ingredients of the Kerr effect, multiphotonionization, and a Drude model for the plasma, are coming under increasing scrutiny. Our new approach answers to this need and provides a computationally viable means of performing simulations of filamentation which includes the quantum coherent atomic effects. Moreover, it eliminates the part of the filamentation model which is considered to be the weakest link, namely ionization and concomitant effects due to freed electrons. For this initial presentation we have chosen to illustrate our computational approach using the example of HHG in a gas jet as this is a distinctly quantum coherent effect and also requites an enormous spectral width. In particular, using the example of nearthreshold HHG in a Xenon gas jet we hope to have shown that our computational approach yields results in keeping with general expectations, and that it is a useful tool to study qualitative dynamic trends in a completely selfconsistent way. As a more concrete example we presented a qualitative comparison between our model and results from an inhouse experiment on extreme ultraviolet generation in a femtosecond enhancement cavity.
Appendix E Appendix
Here we summarize the formulas needed to evaluate the nonlinear component of the current induced in the system described by the timedependent Schrödinger equation (3).
An arbitrary timedependent external field enters through the classical electron trajectory (here we assume that is the direction of the optical field polarization)
(6) 
Assuming that the system was initially in its ground state, an auxiliary quantity is obtained first as a solution to the following integral equation:
(7) 
where the righthandside is
(8) 
Having calculated , the nonlinear (in the strength of the external field ) current is expressed as a sum
(9) 
with the components listed below:
(10) 
(11)  
(12)  
(13) 
Details of the derivation and description of numerical implementation can be found in Ref. [25]. An opensource implementation is also available upon request from the authors (M.K. or J.B.).
Acknowledgments: M.K. and E.M.W. acknowledge support from the U.S. Air Force Office for Scientific Research, through the MURI grant FA95501010561. The parallel pulse propagation framework used in this study was developed by J.A. with the funding from AFOSR grant FA95501110144. R.J.J, E.M.W, and D.R.C. are supported from AFOSR grant FA95501210048.