Optimal Basis Set for Electron Dynamics in Strong Laser Fields: The case of Molecular Ion H2.

A clear understanding of the mechanisms that control the electron dynamics in a strong laser field is still a challenge that requires interpretation by advanced theory. Development of accurate theoretical and computational methods, able to provide a precise treatment of the fundamental processes generated in the strong field regime, is therefore crucial. A central aspect is the choice of the basis for the wave function expansion. Accuracy in describing multiphoton processes is strictly related to the intrinsic properties of the basis, such as numerical convergence, computational cost, and representation of the continuum. By explicitly solving the 1D and 3D time-dependent Schrödinger equation for H2+ in the presence of an intense electric field, we explore the numerical performance of using a real-space grid, a B-spline basis, and a Gaussian basis (improved by optimal Gaussian functions for the continuum). We analyze the performance of the three bases for high-harmonic generation and above-threshold ionization for H2+. In particular, for high-harmonic generation, the capability of the basis to reproduce the two-center interference and the hyper-Raman phenomena is investigated.

Introduction

The optical response of a molecular system to an intense and ultrashort laser pulse is a subject of increasing interest since the advent of the attosecond laser pulses.[1] Recent advances in laser technology are continuously triggering the introduction of new time-resolved spectroscopies, offering the opportunity to investigate electron dynamics in molecules with unprecedented time resolution.[2] For example, electronic charge migrations have been traced in molecules using attosecond pulses,[3] electron correlation effects have been also observed in photoemission processes on the attosecond scale,[4,5] and above-threshold ionization (ATI) together with high-harmonic generation (HHG) spectra have been used to explain the attosecond dynamics of electronic wave packets in molecules.[6,7]

Despite these exciting experimental achievements, reaching a clear understanding of the mechanisms that control the electron dynamics under the action of a strong laser field is still a challenge that requires theoretical support.[6] It is crucial to develop accurate theoretical and computational methods capable to provide precise treatments of the fundamental processes generated by a strong laser field.[8−11]

Nowadays, the electron dynamics problem in strong fields is tackled by two main families of methods: time-dependent density-functional theory (TDDFT) and time-dependent wave function methods.[6,12−16] With these methods, developments have been focused on the accurate description of electron correlation. However, because of the complexity of nonlinear optical phenomena, such as HHG and ATI, another important aspect needs to be carefully addressed: the choice of the one-electron basis for representing the time-dependent wave function. In fact, a reliable description of the electron dynamics in strong laser fields depends on the accuracy in reproducing the bound states and, even more important, the continuum states of the molecular system considered. In addition, choosing a good basis can improve the numerical convergence of the results and reduce the computational cost of simulations.

Most of the proposed numerical methods in literature directly describe the system wave function on a real-space grid[17−20] or through a numerically defined grid-based basis set of functions, as in the case of the discrete-variable representation method,[21] the pseudospectral grid method, or the finite-element method.[22] Within these approaches, schemes have been proposed to compute ATI spectra in molecules[23] and to study the different molecular orbital contributions to HHG spectra.[24,25] Grid-based basis sets have demonstrated to be very accurate to describe nonlinear optical phenomena. However, the computational cost can be very high and strategies involving multilevel parallelization schemes have had to be developed.[26]

Another recurrent basis, in the context of ultrafast electron dynamics, is composed by B-splines, defined as piecewise polynomial functions with compact support.[27] They were first introduced in atomic calculations by Shore[28] and later extensively used to treat ionized and excited states.[29,30] B-splines have proved to be a very powerful tool to describe multiphoton ionization processes in atoms and molecules in the frameworks of TDDFT and wave function methods.[31−34] The success of B-splines is due to a remarkable feature: B-splines are able to reproduce accurately both bound and continuum states. This numerical property is directly related to their effective completeness.[35] Nowadays atomic packages based on B-splines are available[36−38] and recent studies show their ability to reproduce HHG and ATI spectra of molecules under the action of a strong laser field.[39] However, new algorithms have to be developed in order to increase the computational efficiency of complex calculations with B-splines.

More recently, Gaussian-type orbital functions (abbreviated as Gaussian functions in the following), in the framework of the time-dependent configuration-interaction (TDCI) method, have been used to calculate HHG spectra in atoms and molecules.[12,40−43] The importance of the cardinal number (related to the maximal angular momentum) of the basis set and the number of diffuse basis functions was investigated.[12,40] Two strategies to improve continuum states have been studied: multicentered basis functions[12,41] and, alternatively, Gaussian functions with exponents specially optimized to improve the continuum.[42,44] This latter strategy proved to be more efficient than using multicentered basis functions and it has also lower computational cost, however it remains to be tested on molecular systems. These works permitted us to identify the best basis sets to be used in order to capture the features of HHG spectra.

Finally, to overcome some of the limitations of the grid, B-spline, and Gaussian basis, hybrid approaches have been proposed in the last years. For example, Gaussian functions were used together with grid-based functions to reproduce electron dynamics in molecular systems,[45] and also Gaussian functions have been combined with B-splines for studying ionization in H and He atoms.[46,47]

The aim of the present work is to compare the performance of the three families of basis, briefly reviewed above, i.e., grid, B-splines, and Gaussians, for the calculation of HHG and ATI spectra of the molecular ion H2+. This system has been chosen because it has the advantage of having only one electron, which allows us not to bias our investigation with possible effects due to electron correlation. Indeed, with this simple case, we can focus on the effectiveness of the representation of the continuum states for the electron dynamics and the computational advantages of each basis. Moreover, the presence of two nuclei in H2+ offers the opportunity to observe intricate physical features, such as quantum interferences in the HHG process.[48−50]

This article is organized as follows. In Section we present the 1D theoretical model to solve the electronic time-dependent Schrödinger equation (TDSE) with grid, B-spline, and Gaussian bases. In Section we present and discuss the results for the 1D approach. In Section we present the 3D theoretical model to solve the electronic TDSE with grid and Gaussian basis. In Section we present and discuss the results for the 3D approach. We compare the bound and the continuum energy spectra of H2+, as well as HHG and ATI spectra for grid, B-spline, and Gaussian bases, emphasizing the advantages and disadvantages of each representation. In particular, for HHG spectra, we investigate the capability of the different basis to reproduce specific quantum features, such as the hyper-Raman[51] and the two-center interference phenomena.[48−50] Finally, Section contains our conclusions.

1D Theoretical Model of H2+

The electronic TDSE for a 1D model of H2+ is given by, in atomic units (au),where ψ(x, t) is the time-dependent electron wave function. Here, Ĥ0(x) is the field-free Hamiltonian,with a soft Coulomb electron–nuclei interaction given bywhere R is the interatomic distance and α is a parameter chosen to reproduce the exact ionization energy Ip (taken as −1.11 Ha for all the three bases employed here) of the real H2+ molecule at a given value of R (α = 1.44 at R = 2.0 au).[50]

The interaction between the electron and the laser electric field E(t) is taken into account by the time-dependent interaction potential, which is given in the length gauge bywhere E(t) is the laser electric field and x̂ is the electron position operator. The laser electric field is chosen as E(t) = E0f(t) sin(ω0t) where E0 is the maximum amplitude of the pulse, ω0 is the carrier frequency, and f(t) is a trapezoidal envelopewith T0 = 2π/ω0. The duration of the pulse is thus τ = 10T0 (i.e., 10 optical cycles).

HHG and ATI Spectra

A HHG spectrum, experimentally accessible by measuring the emission spectrum in the presence of an intense laser field, can be calculated as the acceleration power spectrum over the duration of the laser pulse τ[50]where −∇V̂ – E(t) is the electron acceleration operator, as defined by the Ehrenfest theorem, and W(t) is an apodization function that we chose to be of the sin-square window form. An alternative way to obtain the HHG spectrum is to calculate the dipole power spectrum as

It can be shown that the two forms are related,[12,52−54] ω4P(ω) ≈ Pa(ω), under reasonable conditions (see Appendix in ref (12)). The function W(t) is a sin-square window function chosen empirically to minimize the noise and especially to remove the artifacts arising from the discrete Fourier transform due to the fact that we integrate only over a limited time duration and not from −∞ to +∞.

An ATI spectrum, which is experimentally accessible by measuring the photoelectron spectrum of the molecule, can be calculated by spectrally analyzing the system wave function ψ(τ) at the time τ corresponding to the end of the laser pulse. Specifically, using the window operator method, one calculates the probability P(E, n, γ) to find the electron in the energy interval [E – γ, E + γ] as[55,56]where γ and n are parameters chosen to allow flexibility in the resolution and accuracy of the energy analysis. In our case we chose n = 2 and γ = 2 × 10–3 au.

Representation of the Time-Dependent Wave Function and Propagation

Real-Space Grid

The time-dependent wave function is discretized on a real-space grid of N points x separated by a constant step Δx = x – x, in the interval [x1 = −(N – 1)Δx/2, x = (N – 1)Δx/2]. It is thus represented by the vectorwhere x = (i – 1 – (N – 1)/2)Δx.

The Laplacian operator is computed with the second-order central difference formula which gives rise to a tridiagonal matrix representation of the Hamiltonian Ĥ0.[17] The TDSE (eq ) is solved by means of the Crank–Nicholson propagation algorithm.[57] The H2+ ground state, computed by inverse iteration,[58] is taken as the initial state for the propagation. In addition, to avoid unphysical reflections at the boundaries of the simulation grid, a mask-type absorber function[17] was implemented with a spatial extension of 50 au.

For ATI spectra, converged results were obtained with N = 200001 and Δx = 0.02 au, and with a time step Δt = 8.41 × 10–4 au. For HHG spectra, we obtained converged results with N = 160001, Δx = 0.01 au, and Δt = 1.35 × 10–2 au.

B-Spline Basis Set

The time-dependent wave function with the B-spline basis set is represented aswhere c(t) are time-dependent coefficients and {B(x)} are a set of B-spline functions of order k and dimension M. To completely define B-spline functions a sequence of knots must be given. Each function B(x) is defined on a supporting interval [t, t] which contains k + 1 consecutive knots, and the function B(x) vanishes outside this interval. We have chosen the first and the last knots to be k-fold degenerate, t1 = t2 = ··· = t = Rmin and t = t = ··· = t = Rmax, while the multiplicity of the other knots is unity. The width of an interval is t – t = Rmax/(M – k + 1).[32] In our calculations we used k = 8, M = 15008, Rmin = 0, and Rmax = 8000 au. The system was placed at the center of the box at x = 4000 au.

ATI and HHG spectra were obtained by solving the TDSE (eq ) within the Cranck–Nicholson propagation algorithm[57] using a time step of Δt = 1.35 × 10–2 au. The H2+ ground state was computed by inverse iteration[58] and taken as the initial state for the propagation. We did not need to use any absorber during the propagation because of the very large size of the simulation box.

Gaussian Basis Set

For the Gaussian basis set we followed the TDCI procedure developed in our previous work[12] and adapted it to the present 1D H2+ model. The time-dependent wave function is represented here aswhere ϕ(x) are the eigenstates of the field-free Hamiltonian Ĥ0, composed by the ground state (k = 0) and all the excited states (k > 0). The ϕ(x) are expanded on the Gaussian basis set. We use uncontracted Gaussians localized on each nucleus and two “angular momenta” (), corresponding to odd and even functions. The basis functions are thus of the form e–α(, where = 0 or 1. The Gaussian exponents α are of two different types. The first type of exponents are optimized to describe the bound part of the wave function. We used the uncontracted STO-3G basis set, i.e., three uncontracted Gaussians whose exponents are taken from the STO-3G basis set with Slater exponent ζ = 1. We take the same exponents for = 0 and = 1. The second type of exponents are optimized for the representation of the continuum.[12] They are computed with the procedure developed by Kaufmann[59] adapted to the 1D model, i.e., by optimizing the overlap between a 1D Slater type function N((ζ)xe–ζ| with ζ = 1 and a Gaussian function , where N( and are normalization factors. Note that, in this case, the exponents used for the = 0 shell and for the = 1 shell are different. In the following, we will denote these Gaussian functions optimized for the continuum as K functions. To sum up, we use three functions with STO-3G exponents and four K functions for each angular momentum, localized on each nucleus, which makes a total of (3 + 4) × 4 = 28 uncontracted Gaussian basis functions. However, when we orthonormalize this basis set, we find linear dependencies that needs to be removed. For this we define a cutoff ε = 10–8 under which the eigenvalues of the overlap matrix are considered to be zero, and their corresponding eigenvectors are removed from the space. We get an orthonormalized basis set of 24 basis functions. The basis-set exponents are collected in Table S1 of Supporting Information. To solve the TDSE (eq ) we used the split-operator propagator with Δt = 1.35 × 10–2 au.

In order to compensate for the unphysical absence of ionization, we used the double-d heuristic lifetime model proposed in ref (12). This model requires two parameters, d0 and d1, which represent different electron escape lengths after ionization. We have chosen these parameters on the basis of the rescattering model[60,61] where an electron is ionized by a strong laser field, accelerated in the continuum, and then brought back close to its parent ion where it can recombine or scatter. From this model, d0 is equal to the maximum electron excursion after ionization which is , while d1 < d0. In our calculations we always used d1 = 20 au. Moreover d0 affects all the continuum states below the cutoff energy Ecutoff = Ip + 3.17Up[60,61] (Up = E02/(4ω02) is the ponderomotive energy of the electron) while d1 handles the ionization for those continuum energy states above Ecutoff. This allows better retention of the contribution of continuum states for the recombination step of the HHG process. Table collects the values of d0 used in this work.

Table 1

d0 Values, Taken as xmax, Used in the Double-d Heuristic Lifetime Model for the Laser Intensities Employed in This Work

I (W/cm2)	d0 (au)
5 × 1013	23
1014	33
2 × 1014	46
3 × 1014	57
4 × 1014	66
5 × 1014	74
7 × 1014	87

There is a fundamental difference between this approach and the grid and B-spline ones. Indeed, the TDSE with the Gaussian basis set is solved in the energy space. This fact permits having a more direct and intuitive interpretation of the role of bound and continuum states in HHG and ATI spectroscopies. In addition, the use of Gaussians reduces considerably the computational time required in time propagation. This makes it a more promising tool for the modelization of larger molecules.

1D Results and Discussion

Spectrum of the Field-Free Hamiltonian

The spectrum of Ĥ0 should be strictly independent of the choice of the basis set in the limit of a complete basis set. However, because our basis sets are not complete, differences in the eigenstates and eigenvalues from grid, B-spline, and Gaussian basis sets can arise, especially at high-energy values. In order to investigate the behavior of the three basis sets, the spectrum of Ĥ0 is analyzed in this section.

In Figure the ground-state wave function is shown. The three basis sets reproduce exactly alike the ground state of the 1D H2+ model, at the equilibrium internuclear distance of R = 2.0 au. The panel (a) of Figure shows the eigenvalues given by each basis set up to the 30th energy state, and in panel (b) of Figure one finds the inverse of the density of continuum states which is defined as ρ(E) = 1/(E – E) where E is a positive eigenvalue. In order to compare the three bases, the density of the states has been normalized to the length of the simulation box in the case of the grid and B-splines and to a constant in the case of the Gaussians. This constant was chosen to force the first Gaussian continuum eigenvalue to match the first continuum eigenvalue of the grid and B-splines, which are identical. For all the three basis sets, the continuum part of the spectrum is represented as a finite number of eigenstates as, in numerical calculations, the basis set is always incomplete. However, the discreteness of the Gaussians is much larger than that of the grid and B-splines. The spectrum obtained with the Gaussians starts to diverge from the grid and B-spline ones already at around the 13th state. This issue is a direct consequence of the relatively small size of the Gaussian basis set compared to the number of grid points or B-spline functions used. Indeed, the STO-3G+4K basis contains only 24 Gaussian basis functions whereas we used 400001 grid points and 15000 B-splines. In principle, we could increase the number of Gaussians, but this will quickly lead to the linear dependency problem. This problem prevents us from using more than a few tens of optimized Gaussian functions. This fact, as we will see in the following sections, can have important consequences on the calculation of HHG and, in particular, of ATI spectra.

Figure 1

Ground-state wave function of H2+ (at the equilibrium internuclear distance of R = 2.0 au) calculated using grid, B-spline, and Gaussian basis.

Figure 2

(a) Eigenvalues of H2+ up to the 30th eigenstate. (b) Inverse of the normalized density of continuum states.

To investigate the accuracy of the grid, B-spline, and Gaussian bases in the description of continuum wave functions, we have chosen two different continuum energies, both representative of two different continuum energy regions: low energy (E = 0.06 Ha) and high energy (E = 1.97 Ha). For each of these energies, we reported in (Figure ) the corresponding wave functions φ(x). For the grid, the continuum wave functions were obtained by propagating the TDSE at the chosen positive energy E with a fourth-order Runge–Kutta algorithm,[58] and then normalized with the Strömgren procedure.[62,63] Instead, for B-splines and Gaussians, the wave functions were obtained from a direct diagonalization of Ĥ0. In this case, the resulting continuum states were renormalized using the procedure proposed by Macías et al.[64,65] We verified that the Strömgren and Macías procedures are equivalent.[66] The continuum wave functions computed with both grid and B-spline basis sets reproduce the same oscillations in the low- and high-energy regions of the continuum. On the other hand, Gaussians can reproduce just a few of the oscillations. We already observed this behavior in the case of the hydrogen atom in a 3D calculation[42] where the crucial role of the K functions was pointed out in order to obtain these oscillations (in that case a much larger basis set was employed). Here, we want to draw the attention to the fact that Gaussians can still be reasonable in the low-energy continuum but become unsuitable to reproduce oscillations for high-energy continuum states. The probability of propagating an electron in one of the two regions depends on the laser parameters used in the simulation. This fact can have important implications in the description of HHG and ATI spectra as we will see in the following sections.

Figure 3

(a) Spatial dependence of the even wave function φ(x) corresponding to E = 0.06 Ha. (b) Spatial dependence of the odd wave function φ(x) corresponding to E = 1.97 Ha.

HHG

HHG spectra have been calculated in the dipole and the acceleration forms for H2+ at different internuclear distances, R = 1.8, R = 2.0 (equilibrium distance), and R = 2.2 au, for a Ti:sapphire laser pulse with a carrier frequency ω0 = 0.057 Ha (1.55 eV, 800 nm) and different intensities, I = 5 × 1013, I = 1 × 1014, I = 2 × 1014, I = 5 × 1014, and I = 7 × 1014 W/cm2.

In Figure we show the dipole form of the HHG spectra at R = 2.0 au for three different laser intensities. All three basis sets reproduce the general expected features of an HHG spectrum: the intensity of the low-order harmonics decreases rapidly, then a plateau region follows where the intensity remains nearly constant, and at high frequencies the harmonic intensity decreases again. As H2+ has a center-of-inversion symmetry, only odd harmonics are presented in the spectrum. We estimated the cutoff energies by calculating Ecutoff = Ip + 3.17Up, as given in the semiclassical rescattering model.[60,61]

Figure 4

HHG spectra calculated from the electron dipole at the equilibrium internuclear distance R = 2.0 au with laser intensities (a) I = 1014 W/cm2, (b) I = 2 × 1014 W/cm2, and (c) I = 5 × 1014 W/cm2. Intensities I = 5 × 1013 and 7 × 1014 W/cm2 are reported in the Supporting Information. For each HHG spectrum, the dot-dashed lines indicate the cutoff energies, which are given by the rescattering model as Ecutoff = Ip + 3.17Up; see refs (60 and 61). (a) Ecutoff = 31.7ω0, (b) Ecutoff = 43.9ω0, and (c) Ecutoff = 80.5ω0. The arrow points to the expected position of the two-center interference minimum extracted from the recombination dipole.

We observe that the grid and B-spline HHG spectra are indistinguishable for all the laser intensities. This fact is consistent with the analysis reported above on the spectrum of Ĥ0 (see Section ). On the other hand, the agreement between the spectra obtained with the Gaussian basis and those obtained with the grid or B-splines deteriorates when the laser intensity increases. This is clearly observed for the plateau region for the intensity I = 5 × 1014 W/cm2 but also detected for the plateau and cutoff regions for the intensity I = 7 × 1014 W/cm2 (see Supporting Information). Most of these observations are also valid when using the acceleration form of the HHG spectrum. The only exception we found was with the Gaussian basis set and laser intensities I = 5 × 1014 W/cm2, as shown in Figure , and I = 7 × 1014 W/cm2 (see Supporting Information). For these largest intensities, the spectrum extracted from the acceleration seems to largely underestimate the position of the cutoff but to much better reproduce the harmonics of the plateau.

Figure 5

HHG spectra calculated from the electron dipole and the electron acceleration at the equilibrium internuclear distance of R = 2.0 au with a laser intensity of I = 5 × 1014 W/cm2 using Gaussian basis sets. The dot-dashed line is the cutoff energy Ecutoff = 80.5ω0 and the arrow points to the expected position of the two-center interference minimum, extracted from the recombination dipole which is identical to the one extracted from the recombination acceleration.

To analyze in more detail the fine structure of the HHG peaks, in Figure we show HHG spectra only up to the 15th harmonics. The B-spline and the grid spectra are almost identical except for some very small differences when the laser intensity is very high. Gaussian spectra reproduce the features of the B-spline and grid ones, but when the laser intensity increases the Gaussian spectrum becomes much more noisy.

Figure 6

HHG spectra calculated from the electron dipole at the equilibrium internuclear distance R = 2.0 au up to the 15th harmonic with laser intensities (a) I = 1014 W/cm2 and (b) I = 5 × 1014 W/cm2. The dashed lines indicate the position of the harmonics while the dotted lines indicate the hyper-Raman lines at position ω̃ ± 2kω0[67] where k is an integer and ω̃ = 6.69ω0 is the resonance with the first excited state.

From panel (a) of Figure it is also possible to identify another series of peaks besides those corresponding to the harmonics. These peaks correspond to hyper-Raman lines with position given by ω̃ ± 2kω0,[67] where k is an integer and ω̃ = 6.69ω0 is the resonance with the first excited state. We observe that the three basis sets describe with the same accuracy the hyper-Raman lines. Moreover, at sufficiently large laser intensity, the HHG process dominates, and the hyper-Raman lines are not observed anymore (panel (b) of Figure ).

The accuracy of the grid, B-spline, and Gaussian calculations was also investigated through their ability to reproduce the two-center interference in the HHG spectrum. This interference was predicted by Lein et al.[50] for diatomic molecules such as H2+. In this model, the electron that recombines with the ionic core can interact with either of the two nuclei. The two atomic centers can therefore be interpreted as coherent point sources, and the whole system can be seen as a microscopic analogue of Young’s two-slit experiment. The light emitted by each nucleus will interfere either constructively or destructively depending on its frequency, and the interference pattern will superimpose to the HHG spectrum. Since Lein’s model has been proposed, a great number of numerical analyses came forth pointing out the role of the internuclear distance, molecular orientation, recombination to excited states, and laser intensity.[11,48,68−75]

According to Lein’s model, the position of the minimum in the spectrum is independent from the laser intensity and can be extracted from the analysis of the recombination dipole drec(E) = ⟨φ0|x̂|φ⟩ where φ0 is the ground state and φ is a continuum state at energy E of Ĥ0. This quantity is plotted in panel (a) of Figure for R = 1.8 au and in panel (a) of Figure for R = 2.2 au. For R = 2.0 au, we report the recombination dipole in the Supporting Information. The minimum described in the two-center interference corresponds to the energy which makes the recombination dipole vanish. We found that the corresponding frequency is ω = 34.0ω0 for R = 1.8 au, ω = 26.4ω0 for R = 2.0 au, and ω = 20.8ω0 for R = 2.2 au. We note that the extraction of the minimum from the recombination dipole is straightforward for the grid and B-spline basis sets, while in the case of the Gaussian basis only a rough estimate can be given. Lein’s model predicts the position of the minimum at ω = π2/(2R2ω0) which gives ω = 26.7ω0 for R = 1.8 au, ω = 21.6ω0 for R = 2.0 au, and ω = 17.9ω0 for R = 2.2 au. The underestimation of the minimum position by Lein’s model has already been pointed out.[70] The main reasons must be searched in the different description of the ground state and the continuum between our 1D theoretical model and Lein’s model.

Figure 7

Two-center interference at R = 1.8 au: (a) recombination dipole and (b) HHG spectrum at I = 2 × 1014 W/cm2. The arrow points to the expected position of the two-center interference minimum extracted from the recombination dipole. The dot-dashed line is the cutoff energy Ecutoff = 43.8ω0. E0 is the ground-state energy.

Figure 8

Two-center interference at R = 2.2 au: (a) recombination dipole and (b) HHG spectrum at I = 2 × 1014 W/cm2. The arrow points to the expected position of the two-center interference minimum extracted from the recombination dipole. The dot-dashed line is the cutoff energy Ecutoff = 43.8ω0. E0 is the ground-state energy.

We report in panel (b) of Figure and in panel (b) of Figure the HHG spectra for R = 1.8 au and for R = 2.2 au with I = 2 × 1014 W/cm2, and we observe that all the basis sets reproduce the position of the minimum of the two-center interference. Also the minimum for R = 2.0 au is very well reproduced as can be seen in Figure . Another observation is that the sharpness of the minimum depends on the laser intensity and on the internuclear distance. We confirm the fact that the minimum is more visible for smaller internuclear distances.[76] We did the same investigation considering the recombination acceleration arec(E) = ⟨φ0|−∇V̂|φ⟩ and the HHG spectrum from the acceleration. We obtained the same results (see Supporting Information) explained before. From these studies we deduce that all the basis sets are capable of accurately reproducing the two-center interference.[50] However, in the case of the Gaussian basis, the acceleration seems to better reproduce the minimum for I = 5 × 1014 W/cm2 (panel (c) of Figure ) and I = 7 × 1014 W/cm2 (see Supporting Information).

From the detailed analysis of HHG spectra presented in this section, we conclude that for a good performance of the Gaussian basis the laser intensity cannot be “very large”. For example, for intensity lower than I = 5 × 1014 W/cm2 we obtain correct HHG spectra while for higher intensities only the harmonic peaks in the low-energy part of the plateau are correct. A strategy to improve the Gaussian basis set could be to modify the cutoff ε below which the eigenvalues of the overlap matrix are set to zero. This will change the number of kept eigenvectors. In Figure we compare an HHG spectrum for I = 5 × 1013 W/cm2 calculated with the grid and with the Gaussian basis while changing the linear-dependency threshold ε: ε = 10–4 (17 basis functions), ε = 10–8 (24 basis functions, which is the standard choice throughout the article), and ε = 10–10 (26 basis functions). This analysis shows that for a “low” intensity (I = 5 × 1013 W/cm2) the quality of the HHG spectrum in the plateau and cutoff regions is not affected by the specific choice of the threshold of eigenvalues.

Figure 9

HHG spectra from the dipole at the equilibrium internuclear distance R = 2.0 au with I = 5 × 1013 W/cm2 obtained with the grid and with the Gaussian basis sets with linear-dependency thresholds ε = 10–4, ε = 10–8, and ε = 10–10.

ATI

We calculated ATI spectra with intensities I = 5 × 1013, 1 × 1014, and 5 × 1014 W/cm2. In panel (a) of Figure we show the ATI spectrum with laser intensity I = 1014 W/cm2, while the spectra for intensities I = 5 × 1013 and 5 × 1014 W/cm2 are reported in the Supporting Information.

Figure 10

(a) ATI spectrum calculated at the equilibrium interatomic distance R = 2.0 au with intensity I = 1 × 1014 W/cm2. (b) Photoelectron spectrum calculated with the Gaussian basis at the equilibrium distance R = 2.0 au with intensity I = 1 × 1014 W/cm2 and three photon energies ω0 = 1.34 Ha (black), ω0 = 1.47 Ha (red), and ω0 = 1.61 Ha (blue). The ground-state energy (−1.11 Ha) and the continuum-state energies (0.06 Ha, 0.22 Ha, and 0.50 Ha) which correspond to transitions allowed by symmetry are displayed (green dots).

The ATI spectrum of Figure has positive energy peaks (bound-continuum transitions) corresponding to the electron density ionized during the propagation, i.e., the photoelectron spectrum, while the peaks in the negative region (bound–bound transitions) represent the electron density remaining in the ground state and that has been transferred to excited states. We remember that only the positive energy region of an ATI spectrum is experimentally measurable.

As already seen for the HHG spectra, the grid and B-spline basis sets describe with the same accuracy both bound–bound and bound–continuum transitions. Their ATI spectra coincide and correctly reproduce the expected features of an ATI spectrum: the distance between two consecutive ATI peaks (in the positive energy region) is constant and equal to the energy of a photon, i.e., 0.057 Ha.

The Gaussian basis is only able to reproduce bound–bound transitions. The negative energy part of the spectrum is quite close to the one obtained with the grid and B-splines, while bound-continuum transitions are out of reach for the Gaussian basis set. This limitation is due to the low density of states in the continuum. Indeed, with the basis-set parameters used here, only six continuum states are reproduced in the energy region between 0 and 1 Ha, as we can see in the bottom panel of Figure . This low density of states is far from reproducing the correct ATI energy distribution and explains why no more than six peaks are observed in the positive energy region of the spectrum. The energies of the six ATI peaks correspond to the energies of the six continuum states reported in Figure . To detail more on this feature, we plot in panel (b) of Figure the photoelectron spectrum, computed with the Gaussian basis, after absorption of one photon and for three different photon energies ω0 = 1.34 Ha, ω0 = 1.47 Ha, and ω0 = 1.61 Ha. Together, we also plot the energy position of the ground state and of the first continuum energies corresponding to symmetry-allowed transitions. One clearly sees that if the photon energy matches the energy of a transition from the ground state to one of the continuum states, then we get a photoelectron peak. However, if the photon energy does not match any transition then no ionization is observed. This crucial feature forbids the computation of a correct photoelectron or ATI spectrum with the Gaussian basis set used here. We believe that larger Gaussian basis sets can in principle describe ATI. Indeed, in 3D calculations,[12] one can easily produce tens of low-energy (<1 Ha) continuum states, leading to a possible improvement of the ATI spectrum.

3D Theoretical Model of H2+

The electronic TDSE for a 3D model of H2+ is given by, in atomic units (au),where ψ(r, t) is the time-dependent electron wave function. Here, Ĥ0(r) is the field-free Hamiltonian,with V̂(r) the Coulomb electron–nuclei interaction.

The interaction between the electron and the laser electric field E(t) is taken into account by the time-dependent interaction potential, which is given in the length gauge bywhere E(t) is the laser electric field polarized along the z axis, corresponding to the H2+ internuclear axis, and ẑ is the electron position operator along this axis. We have chosen the same type of laser as in the 1D model (see Section ) except that the duration of the pulse is τ = 6T0 (i.e., 6 optical cycles). We calculated HHG spectra from the dipole as in eq .

Representation of the Time-Dependent Wave Function and Propagation

Concerning the 3D calculations on a grid, we used the Octopus code which is a software package for TDDFT calculations.[26] For our calculations we have chosen the “independent particle” option which allows getting the numerically exact solution for the TDSE in the case of one electron. We have chosen as the simulation box a cylinder with radius 50 au and height 100 au with a grid space Δr = 0.435 au. The TDSE of eq is solved by means of the Crank–Nicholson propagation algorithm[57,58] with a time step Δt = 5 × 10–2 au. Also in this case to avoid unphysical reflections at the boundaries of the simulation box, a mask-type absorber function was used with a spatial extension of 22 au.

Gaussian Basis Set

In this case, we used the approach we developed and detailed in refs (12 and 40) which consists of solving the TDSE using the TDCI approach. For the Gaussian calculations, we used a development version of the MOLPRO software package[77] and the external code LIGHT[40] to perform the time propagation using also in this case a time step Δt = 5 × 10–2 au. As Gaussian basis set we used a 6-aug-cc-pVTZ with 5 K functions, which we denote as 6-aug-cc-pVTZ+5K, which is the largest basis without linear dependencies. The basis-set exponents and contraction coefficients are collected in Table S2 of Supporting Information. To treat ionization we used a double-d heuristic model where the parameters d1 and d0 have been chosen as in the 1D model. The value of Ip is in this case −1.10 Ha.

3D Results and Discussion

We calculated HHG spectra in the dipole form for H2+ at internuclear distance R = 2.0 au (equilibrium) for a Ti:sapphire laser with a carrier frequency ω0 = 0.057 Ha and intensities I = 5 × 1013, 1 × 1014, 2 × 1014, 3 × 1014, 4 × 1014, and 5 × 1014 W/cm2.

In Figure we show the HHG spectra for three laser intensities (the spectra for the other intensities are reported in the Supporting Information). Both the Gaussian and grid basis sets reproduce well the expected features of an HHG spectrum, regardless of the applied field intensity, as already pointed out for the 1D case. However, starting from intensity I = 3 × 1014 W/cm2, the quality of the spectrum obtained with the Gaussian basis set tends to diminish, especially in the cutoff region. For 3D calculations, obtaining a good HHG spectrum with optimized Gaussians seems to be more difficult than for 1D calculations, due to the computational complexity.

Figure 11

HHG spectra in the dipole form at the equilibrium internuclear distance R = 2.0 au with laser intensities (a) I = 5 × 1013 W/cm2, (b) I = 2 × 1014 W/cm2, and (c) I = 3 × 1014 W/cm2. For each HHG spectrum, the dot-dashed line gives the cutoff energy Ecutoff = Ip + 3.17Up given by the rescattering model[60,61] which is (a) Ecutoff = 25.4ω0, (b) Ecutoff = 43.7ω0, and (c) Ecutoff = 55.9ω0. The arrow points to the expected position of the two-center interference minimum extracted from the recombination dipole.

However, it is interesting to note that the low-energy harmonics are still well described when compared to the grid calculations. We show this behavior by analyzing the fine structures of the peaks as shown in Figure . Here, we plot the HHG spectra up to the 13th harmonic for different intensities. For the grid calculations (panel (a)) with I = 5 × 1013 W/cm2 only the first and the third harmonic peaks are clearly visible together with a strong and large peak at around 7.65ω0, due to the emission from the first excited state. Also in this case we observe hyper-Raman lines at position ω̃ ± 2kω0[67] where k is an integer and ω̃ = 7.65ω0 is the resonance with the first excited state. Observing the evolution of the harmonics and the resonant peaks as a function of the laser intensity (from I = 5 × 1013 W/cm2 to I = 5 × 1014 W/cm2), the harmonics become more and more intense while the hyper-Raman lines almost disappear. The same behavior was already observed in the 1D model. The spectra obtained with the Gaussian basis set show exactly the same trend as shown in panel (b) of Figure .

Figure 12

HHG spectra in the dipole form at the equilibrium internuclear distance of R = 2.0 au up to the 13th harmonic with laser intensities I = 5 × 1013 W/cm2, I = 2 × 1014 W/cm2, and I = 5 × 1014 W/cm2 for (a) grid and (b) Gaussian basis sets. For each HHG spectrum, the dashed line indicates the position of the harmonics, and the dotted line indicates the hyper-Raman lines at position ω̃ ± 2kω0[67] where k is an integer and ω̃ = 7.65ω0 is the resonance of the first excited state.

Conclusions

We explicitly solved the 1D and 3D TDSE for H2+ in the presence of an intense electric field, and we explored the numerical performance of using a real-space grid, a B-spline basis, or a Gaussian basis optimized for the continuum. We analyzed the performance of the three basis sets for calculating HHG and ATI spectra. In particular, for HHG, the capability of the basis set to reproduce the two-center interference and the hyper-Raman lines was investigated. We showed that the grid and B-spline representations of the time-dependent wave function give the same results for both HHG and ATI. On the contrary, the performance of the Gaussian basis is more mixed and depends on the intensity of the laser. It is possible to optimize Gaussian functions to describe the low-energy part of the continuum. However, this optimization is limited by the issue of linear dependencies among Gaussian functions. This implies that for HHG the Gaussian basis can perform well up to the laser intensity I = 5 × 1014 W/cm2 for 1D and up to I = 2 × 1014 W/cm2 for 3D. For higher intensities we have found that only low-energy harmonics are still correct. Moreover, for 3D calculations, obtaining a good HHG spectrum with optimized Gaussian functions seems to be more difficult than in 1D calculations. Despite their limitations, Gaussian basis sets can reproduce intricate features of the HHG spectrum at low energy. Instead, in the case of ATI, Gaussian basis sets make impossible the description of a correct spectrum.

In conclusion, from our investigation we noticed that the grid and B-spline basis sets have very similar behavior and computational cost. These basis sets are very accurate to describe the continuum and phenomena such as HHG and ATI. Gaussian basis sets are less efficient to describe the continuum. The effect on ATI and HHG spectra is however different: on one hand, ATI spectrum is not reproduced by Gaussian basis functions, but on the other hand the most important features and fine structures (minimum/resonances) at low energy of the HHG spectrum are correctly described. A clear advantage of Gaussian functions with respect the other basis sets is their computational cost which continues to make them interesting for many-electron systems.