Space-time resolved simulation of femtosecond nonlinear light-matter interactions using a holistic quantum atomic model : Application to near-threshold harmonics

M. Kolesik; E. M. Wright; J. Andreasen; J. M. Brown; D. R. Carlson; R. J. Jones

doi:10.1364/OE.20.016113

1. Introduction

The goal of this paper is to present a new computational approach for nonlinear light-matter 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 space-time resolved computational approaches already exist and incorporate the linear dispersive properties of the medium, both refractive-index and absorption, along with a model for the bound state electronic nonlinearity, typically a Kerr-type nonlinearity possibly with saturation, and also a time-delayed 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 strong-field 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 ionization-rate 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 high-harmonic 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 three-dimensional (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 one-dimensional, 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 light-matter interactions at femtosecond time scales and ionization-level intensities. However, the model has the virtue that it can be implemented to efficiently simulate space-time 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 time-dependent 1D atomic model. In a series of recent papers Bandrauk and co-workers [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. One-dimensional atomic models have played an important role as model systems for elucidating the physics underlying nonlinear light-matter interactions. For example, Su and Eberly [9] introduced the quasi-Coulomb potential in 1D as a model for multi-photon 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–19], the 1D δ-potential has also been applied in semiconductors [11,13], and used as a test system for numerical and approximate methods [20–22]. We have applied the δ-potential model to explore the relative roles of the higher-order 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 time-dependent light field [25]. (The code computing the nonlinear current is available as open-source upon request to the corresponding author (MK).) This exact solution side steps the need for direct evaluation of the time-dependent 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 near-threshold HHG in Xenon as an illustrative example. We chose near-threshold HHG, with concomitant lower-orders 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 strong-field 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 phase-matching 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 strong-field approximation, but rather are intended to highlight the capabilities of the approach. Sec. 3 also includes what we believe are new results to model near-threshold 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 in-house experiment that has recently been used to produce record power levels of extreme ultraviolet (XUV) [26].

2. Computational approach and model

To contrast the model of this paper to the current state of the art in computer simulation of high-harmonic generation and/or femtosecond filamentation, we briefly recall the main features of the standard approach. The propagation of the optical and high-harmonic fields are described in one-way (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 z-axis takes the general form [27]

\partial_{z} E_{k_{⊥}} (ω, z) = i k_{z} E_{k_{⊥}} (ω, z) + \frac{i ω^{2}}{2 ε_{0} c^{2} k_{z}} P_{k_{⊥}} (ω, z, {E}) - \frac{ω}{2 ε_{0} c^{2} k_{z}} J_{k_{⊥}} (ω, z, {E}),

where the z-component of the optical wavevector is given by

k_{z} (ω, k_{⊥}) = \sqrt{ω^{2} (1 + χ (ω)) / c^{2} - k_{⊥}^{2}} .

The UPPE describes the evolution of the spatially resolved spectral amplitudes of the scalar electric field E_{k_⊥} (ω, z), where ω is the frequency, and k_⊥ denotes transverse wavenumbers, so that the space-time dependent electric field E(r,z,t) may always be reconstructed using a 3D Fourier (or Hankel) transform over k_⊥, ω. 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 (P) or current-density (J) terms. These are both evaluated as functionals of the electric field history E(t) at each spatial point, and then Fourier-transformed to the spectral representation, e.g. J_{k_⊥} (ω, z, {E}), that appears in the above propagation equation. The linear wavelength-dependent 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 freed-electron 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 Drude-plasma model [1].

Our computational approach for simulating femtosecond nonlinear light-matter 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 complex-valued frequency dependent susceptibility χ(ω).
The 1D δ-potential quantum model for modeling the nonlinear current (J) 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 Kerr-like nonlinearity and associated nonlinear polarization (P) to capture the nonlinear optical response due to the ground state to excited bound state transitions. This will contribute processes such as four-wave mixing, self-phase modulation, and self-focusing. (For a molecular gas we may also add a time-delayed 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 near-threshold HHG in Xenon gas. The three medium components χ, J, P are coupled into the field-propagation 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 pulse-propagation 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 phase-matching and re-absorption 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 k_z(ω, k_⊥) 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 refractive-index 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 35^th 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 Kramers-Kronig 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.

Fig. 1 Real (dark lines) and imaginary (red line) parts of the linear susceptibility of Xenon as employed in the numerical simulations in this paper. The vertical dashed lines mark the fundamental wavelength of the fundamental field at 800 nm, and several of its harmonic are also indicated.

Download Full Size | PDF

Quantum atomic model for the nonlinear response

As alluded to in the introduction our computational approach employs a 1D quantum model with a short-range δ-potential for the nonlinear current due to ground state to continuum transitions. More specifically the corresponding time-dependent Schrödinger equation in atomic units takes the form

[i \partial_{τ} + \frac{1}{2} \partial_{s}^{2} + B δ (s) - s F (τ)] ψ (s, τ) = 0.

Here s and τ are the space and time variables in atomic units, the strength of the attractive potential is governed by the parameter B, with the ionization potential of the ground state equal to B²/2, and F(τ) denotes the time-dependent 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 time-dependent 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

J (τ) \propto \int_{- \infty}^{\infty} d s ψ^{*} (s, τ) (i \frac{\partial}{\partial s}) ψ (s, τ),

without having to go through solving Eq. (3) for the wave function ψ. Rather, the current J required for the UPPE is calculated directly, completely by-passing 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 open-source 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 F(t) 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 right-hand-side of the UPPE equation. We remark that since we only incorporate the nonlinear current from the quantum model into the UPPE, we do not double-count 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 Hydrogen-like 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 ten-cycle pulse of center wavelength λ = 800nm and peak intensity1.5 × 10¹⁸W/m². This harmonic spectrum exhibits the characteristic high-harmonic generation plateau with a high frequency cut-off. The cut-off predicted by the formula from the standard HHG model [2, 29, 30] is indicated by the arrow in Fig. 2 for our parameters. Furthermore, the harmonic spectrum includes the 9^th harmonic which occurs near the ionization potential up to around the 40^th harmonic, and thus constitutes an example of near-threshold HHG where the generated harmonics straddle the ionization energy. What is shown in Fig. 2 is the spectrum of the nonlinear polarization P 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.

Fig. 2 High-harmonics in the spectrum of the dipole moment induced by a pulse at λ = 800nm. The intensity was 1.5 × 10¹⁸W/m², kept constant over duration of ten optical cycles. The arrow marks the location of cut-off energy calculated for these conditions.

Download Full Size | PDF

Finally we point out some pros and cons of our 1D model versus the strong-field 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 well-defined 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 low-frequency 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

P (t) = 2 ε_{0} n_{b} {\bar{n}}_{2} E^{2} (t) E (t),

where n̄₂ is the nonlinear index and n_b is the linear refractive index. This effect enters the propagation equation via the nonlinear polarization term, and gives rise to self-focusing and a cascade of lower-harmonic radiation.

3. Illustrations involving near-threshold high-harmonic generation

HHG in a Xenon gas jet

The first illustrative example involves near-threshold HHG in a Xenon gas jet. The incident pulse in each case is Gaussian in space and time with pulse duration t_p = 85 fs and center wavelength of 800 nm, and focused spot size w₀. We model the Xenon gas jet as a vertical column which is tapered along the propagation or z-axis 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 L = 200 μm between z = 100–300 μ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 z = 100 μm close to the entrance to the gas jet, and also the case that the input beam is focused at z = 300 μ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 on-axis intensity along the z-axis is illustrated in Fig. 3 for the cases of a) a narrow beam with w₀ = 7 μm, and b) a wide beam with w₀ = 15 μm.

Fig. 3 For the sake of comparison, we examine four cases of focusing geometry. Here, A and B denote geometric focus after and before the center of the gas jet. The nominal maximal intensity is kept constant.

Download Full Size | PDF

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 Fig. 4 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 (3 × 3) 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.

Fig. 4 Angularly resolved spectra of the 9th, 11th, and 13th harmonics at three different propagation distances through the gas jet for w₀ = 7 μm. Each panel shows a frequency region (vertical axis) corresponding to one half of the harmonic order, centered at the given harmonic frequency. The horizontal extent of each panel corresponds to the transverse wave-number and thus represents the angle of propagation.

Download Full Size | PDF

Tighter focus case

For the first example we chose a focused spot size of w₀ = 7 μm giving a Rayleigh range of $z_{0} = π w_{0}^{2} / λ = 192 μ m$ for the input field. In this case L ≃ z₀ and we expect that the focusing of the beam will impact phase-matching of the HHG via the Gouy phase-shift 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 phase-matching [33–35]. For example, on-axis phase-matching requires cancellation between the accumulated Gouy phase-shift and the accumulated atomic phase-shift due to the quantum path taken by the electron. The Gouy phase-shift is a positive quantity whereas the quantum phase is proportional to the z-derivative 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 on-axis phase matching is not a possibility for the case (a) with the input beam focused at the exit, whereas phase-matching 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 phase-matching 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 phase-matching conditions: For case (a) for which strict on-axis phase-matching is not possible we see that the off-axis HHG emission is more significant with respect to the on-axis emissions for each harmonic order than for case (b) where on-axis phase-matching is possible. This is particularly pronounced for the 13^th harmonic where case (a) is dominated by off-axis emission whereas case (b) is dominantly on-axis emission. So the results of our simulations are compatible with expectations based on on-axis phase-matching [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 [4, 34]. Furthermore, Fig. 4 reveals that the evolution of the harmonics during the pulse propagation through the jet is distinctly non-monotonic, 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 phase-matching 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 9^th harmonic at three propagation distances. At the shortest and longest propagation distances note that the spectrum has its maximum on-axis, whereas at the on-axis and off-axis 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 200 μm is therefore dominated by a reduction in the off-axis as opposed to the on-axis emission. This highlights the point that the spatial evolution of the various harmonics involves a complex interplay between on-axis and off-axis emissions, and it is not generally sufficient to assess the conversion efficiency based on the on-axis phase-matching considerations.

Fig. 5 Evolution of the conversion efficiency versus propagation distance in the gas jet for two different beam focusing conditions. The gas density is a smooth flat-top profile with a roughly constant pressure region between 100 and 300 micron. The pulse intensities are comparable in both cases, with their maxima approximately located at z = 300 and z = 100 micron in case a) and b), respectively.

Download Full Size | PDF

Fig. 6 On-axis and off-axis propagating harmonic components exhibit different phase-matching behavior. The top panel show the angularly resolved spectra of the 9th harmonic at the beginning, in the middle and at the end of the gas jet corresponding to case (b) from Fig. 4. Energy accumulates preferentially in the off-axis (conical) components during the propagation in the center of the jet. Later, the on-axis component becomes relatively stronger, as the pulse starts to exit from the jet and the gas pressure it experiences starts to decrease. This behavior correlates with the evolution of the energy accumulated in the 9th harmonic (shown as the red line in the lower panel). Note that this is the harmonic which suffers the strongest losses.

Download Full Size | PDF

Weaker focus case

For this case we chose a focused spot size of w₀ = 15 μm giving a Rayleigh range of $z_{0} = π w_{0}^{2} / λ = 883 μ m$ for the input field. Now, L << z₀ and we expect that the focusing of the beam will not impact phase-matching 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 far-field patterns).

Fig. 7 As expected, for a wider beam (w₀ = 15μm), the precise location of the beam focus is less important. Angularly resolved spectra are shown for the 9th, 11th, and 13th harmonics at three different propagation distances through the gas jet. Each panel shows a frequency region (vertical axis) corresponding to one half of the harmonic order, centered at the given harmonic frequency. The horizontal extent of each panel corresponds to the transverse wave-number and thus represents the angle of propagation.

Download Full Size | PDF

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.

Fig. 8 Evolution of the conversion efficiency versus propagation distance in the gas jet for two different beam focusing conditions. The gas density is a smooth flat-top profile with a roughly constant pressure region between 100 and 300 micron. The pulse intensities are comparable in both cases, with their maxima approximately located at z = 300 and z = 100 micron in case a) and b), respectively.

Download Full Size | PDF

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 strong-field approximation [31, 34]. Moreover the appearance of these angle-frequency 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 angle-frequency spectra can be attributed to the interferences in the single-atom response in the transverse cross-section 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 near-field cross-sections (not shown) of the frequency-filtered atomic response shows that complex far-field 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 w₀ = 15 μ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 1 × 10¹⁴ 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 steady-state 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 round-trip 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 free-electron 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 cross-section of the beam. Of course, the spectral nature of our pulse solver allows to endow this susceptibility contribution with the correct frequency dependence, $χ (ω) ~ - ω_{plas}^{2} / ω^{2}$ .

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.

Fig. 9 Harmonic spectrum calculated (left panels) for conditions reflecting those in the femtosecond enhancement cavity. Dashed and full line compare results for simulation with and without inclusion of losses. Left most panel: beam focus at the “entrance” into the gas jet. Middle panel: beam focus at the “exit” from the gas jet. Experimental spectrum shown in the right panel, with an account for the grating efficiency.

Download Full Size | PDF

We certainly do not expect that the simplistic one-dimensional quantum model put forward here could reproduce experimental results in a one-to-one manner. However, we see the above comparison as very encouraging — it shows that the model, especially when coupled with good-quality 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.

4. 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 high-harmonic generation and strong-field ionization in a holistic fashion. The key innovation that makes this possible is the use of 1D atomic model with delta-function 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 light-matter interactions can play a significant role. The current models, based upon individual ingredients of the Kerr effect, multi-photon-ionization, 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 near-threshold 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 self-consistent way. As a more concrete example we presented a qualitative comparison between our model and results from an in-house experiment on extreme ultraviolet generation in a femtosecond enhancement cavity.

5. Appendix

Here we summarize the formulas needed to evaluate the nonlinear component of the current induced in the system described by the time-dependent Schrödinger equation (3).

An arbitrary time-dependent external field F(t) enters through the classical electron trajectory x_cl(t) (here we assume that x is the direction of the optical field polarization)

p_{c l} (t) = - \int_{0}^{t} F (τ) d τ, x_{c l} (t) = \int_{0}^{t} p_{c l} (τ) d τ .

Assuming that the system was initially in its ground state, an auxiliary quantity A(t) is obtained first as a solution to the following integral equation:

A (t) = ψ_{R} (- x_{c l} (t), t) + \frac{i B}{\sqrt{2 π i}} \int_{0}^{t} d t^{'} \frac{e^{+ i \frac{B^{2}}{2} (t^{'} - t)}}{\sqrt{t - t^{'}}} exp [\frac{i {(x_{c l} (t) - x_{c l} (t^{'}))}^{2}}{2 (t - t^{'})}] A (t^{'}),

where the right-hand-side is

ψ_{R} (x, t) \equiv \frac{e^{+ B x}}{2} erfc (\frac{i B t + x}{\sqrt{2 i t}}) + \frac{e^{- B x}}{2} erfc (\frac{i B t - x}{\sqrt{2 i t}}) .

Having calculated A(t), the nonlinear (in the strength of the external field F) current is expressed as a sum

J^{(n l)} = J_{S S}^{(n l)} + J_{F S}^{(n l)},

with the components listed below:

J_{S S}^{(n l)} = 2 Im {\int_{0}^{t} d t_{1} \int_{0}^{t_{1}} d t_{2} \frac{{(- i)}^{\frac{3}{2}} B^{3}}{\sqrt{2 π}} \frac{e^{+ i \frac{B^{2}}{2} (t_{2} - t_{1})}}{\sqrt{t_{1} - t_{2}}} [e^{\frac{i {[x_{c l} (t_{1}) - x_{c l} (t_{2})]}^{2}}{2 (t_{1} - t_{2})}} A^{*} (t_{1}) A (t_{2}) - 1] \frac{x_{c l} (t_{1}) - x_{c l} (t_{2})}{t_{1} - t_{2}}}

\begin{array}{l} J_{F S}^{(n l)} & = & Im {i B^{3} \int_{0}^{t} d t_{1} A^{*} (t_{1}) e^{+ B x_{c l} (t_{1})} erfc (\frac{(1 + i) (B t_{1} - i x_{c l} (t_{1}))}{2 \sqrt{t_{1}}})} \\ - & Im {i B^{3} \int_{0}^{t} d t_{1} A^{*} (t_{1}) e^{- B x_{c l} (t_{1})} erfc (\frac{(1 + i) (B t_{1} + i x_{c l} (t_{1}))}{2 \sqrt{t_{1}}})} \\ - & Im {2 B^{3} \int_{0}^{t} d t_{1} x_{c l} (t_{1}) (i B erfc (\frac{(1 + i) B \sqrt{t_{1}}}{2}) - \frac{1 + i}{\sqrt{π t_{1}}} e^{- i \frac{B^{2}}{2} t_{1}})} \end{array}

Details of the derivation and description of numerical implementation can be found in Ref. [25]. An open-source 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 FA9550-10-1-0561. The parallel pulse propagation framework used in this study was developed by J.A. with the funding from AFOSR grant FA9550-11-1-0144. R.J.J, E.M.W, and D.R.C. are supported from AFOSR grant FA9550-12-1-0048.

References and links

1. A. Couairon and A. Mysyrowicz, “Femtosecond filamentation in transparent media,” Phys. Rep. 441, 47–189 (2007). [CrossRef]

2. M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, and P. B. Corkum, “Theory of high harmonics generation by low-frequency laser fields,” Phys. Rev. A 49, 2117–2132 (1994). [CrossRef] [PubMed]

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). [CrossRef]

4. M. Gaarde, J. Tate, and K. Schafer, “Macroscopic aspects of attosecond pulse generation,” J. Phys. B: At. Mol. Opt. Phys. 41, 132001 (2008). [CrossRef]

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). [CrossRef]

6. E. A. Volkova, A. M. Popov, and O. V. Tikhonova, “Nonlinear polarization response of an atomic gas medium in the field of a high-intensity femtosecond laser pulse,” JETP Lett. 94, 519–524 (2011). [CrossRef]

7. I. Christov, “Propagation of ultrashort pulses in gaseous medium:breakdown of the quasistatic approximation,” Opt. Express 6, 34–39 (2000). [CrossRef] [PubMed]

8. E. Lorin, S. Chelkowski, and A. D. Bandrauk, “The WASP model: A Micro-Macro system of Wave-Schrödinger-Plasma 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). [CrossRef] [PubMed]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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). [CrossRef]

15. V. M. Villalba and L. A. González-Dí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). [CrossRef]

16. S. Geltman, “Ionisation dynamics of a model atom in an electrostatic field,” J. Phys. B: At. Mol. Phys. 11, 3323–3337 (1978). [CrossRef]

17. G. V. Dunne and C. S. Gauthier, “Simple soluble molecular ionization model,” Phys. Rev. A 69, 053409 (2004). [CrossRef]

18. Q. Su, B. P. Irving, C. W. Johnson, and J. H. Eberly, “Stabilization of a one-dimensional short-range model atom in intense laser fields,” J. Phys. B: At. Mol. Opt. 29, 5755–5764 (1996). [CrossRef]

19. T. Dziubak and J. Matulewski, “Stabilization of one-dimensional soft-core and singular model atoms,” Eur. Phys. J. D 59, 321–327 (2010). [CrossRef]

20. A. Sanpera, Q. Su, and L. Roso-Franco, “Ionization suppression in a very-short-range potential,” Phys. Rev. A 47, 2312–2318 (1993). [CrossRef] [PubMed]

21. S. Geltman, “Short-pulse model-atom studies of ionization in intense laser fields,” J. Phys. B - At. Mol. Opt. 27, 1497–1514 (1994). [CrossRef]

22. Z. X. Zhao, B. D. Esry, and C. D. Lin, “Boundary-free scaling calculation of the time-dependent Schrödinger equation for laser-atom interactions,” Phys. Rev. A 65, 023402 (2002). [CrossRef]

23. A. Teleki, E. M. Wright, and M. Kolesik, “Microscopic model for the higher-order nonlinearity in optical filaments,” Phys. Rev. A 82, 065801 (2010). [CrossRef]

24. J. M. Brown, E. M. Wright, J. V. Moloney, and M. Kolesik, “On the relative roles of higher-order nonlinearity and ionization in ultrafast light-matter interactions,” Opt. Lett. 37, 1604–1606 (2012). [CrossRef] [PubMed]

25. J. Brown, A. Lotti, A. Teleki, and M. Kolesik, “Exactly solvable model for non-linear light-matter interaction in an arbitrary time-dependent field,” Phys. Rev. A 84, 063424 (2011). [CrossRef]

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). [CrossRef] [PubMed]

27. M. Kolesik and J. V. Moloney, “Nonlinear optical pulse propagation simulation: From Maxwell’s to unidirectional equations,” Phys. Rev. E 70, 036604 (2004). [CrossRef]

28. A. Bideau-Mehu, Y. Guern, R. Abjean, and A. Johannin-Gilles, “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). [CrossRef]

29. P. B. Corkum, “Plasma perspective on strong field multiphoton ionization,” Phys. Rev. Lett 71, 1994–1997 (1993). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

31. F. Schapper, M. Holler, T. Auguste, A. Zaïr, M. Weger, P. Salières, L. Gallmann, and U. Keller, “Spatial finger-print of quantum path interferences in high order harmonic generation,” Opt. Express 18, 2987–2994 (2010). [CrossRef] [PubMed]

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 semi-infinite gas cell,” Laser Phys. Lett. 9, 207–211 (2012). [CrossRef]

33. A. L’Huillier, K.J. Schafer, and K.C. Kulander, “High-order harmonic generation in Xenon at 1064 nm: The role of phase matching,” Phys. Rev. Lett. 66, 2200–2203 (1991). [CrossRef]

34. P. Salières, A. L’Huillier, and M. Lewenstein, “Coherence control of high-order harmonics,” Phys. Rev. Lett. 743776–3779 (1995). [CrossRef] [PubMed]

35. E. Constant, D. Garzella, P. Breger, E. Mvel, 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). [CrossRef]

36. R. J. Jones, K. D. Moll, M. J. Thorpe, and J. Ye, “Phase-coherent frequency combs in the vacuum ultraviolet via high-harmonic generation inside a femtosecond enhancement cavity,” Phys. Rev. Lett. 94, 193201 (2005). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

38. D. R. Carlson, J. Lee, J. Mongelli, E. M. Wright, and R. J. Jones, “Intracavity ionization and pulse formation in femtosecond enhancement cavities,” Opt. Lett. 36, 2991–2993 (2011). [CrossRef] [PubMed]

Space-time resolved simulation of femtosecond nonlinear light-matter interactions using a holistic quantum atomic model : Application to near-threshold harmonics

Abstract

1. Introduction

2. Computational approach and model

Linear dispersive properties

Quantum atomic model for the nonlinear response

The nonlinear Kerr effect

3. Illustrations involving near-threshold high-harmonic generation

HHG in a Xenon gas jet

Tighter focus case

Weaker focus case

HHG in a femtosecond enhancement cavity

4. Conclusion

5. Appendix

Acknowledgments

References and links

Cited By

Figures (9)

Equations (11)

Optics Express