Core-polarization effects in molecular high harmonic generation

By comparing three-dimensional multi-electron calculations with several simplifying models, we show that polarization of the molecular ion can crucially contribute to the high harmonic spectra of molecules. We solve the time-dependent Schrödinger equation for diatomic model molecules with up to 4 active electrons using the multi-configuration time-dependent Hartree–Fock method. Single electron models fail to reproduce the harmonic spectra as they do not account for polarization of the multi-electron ionic core by the laser field. We find the dominant mechanism in the modification of the recombination dipole matrix element, while ionization is only mildly affected.


Introduction
The radiative response of molecular gases to a strong laser pulse contains extensive information about the shape of the laser pulse, the structure of the gas molecules and the macroscopic propagation of the radiation. For the short wavelength part one can separate propagation effects from the microscopic response of the molecules and interpret the measured radiation in terms of the intriguing three-step model of high harmonic generation (HHG) [1]- [3]. In this picture, the harmonic response of a single molecule is split into three independent factors: initial release of an electron by tunnel-ionization, acceleration of a quasi-free electron in the laser field and recombination of the electron into the state from which it was originally released. Here, release and acceleration are usually considered to be largely independent of atomic or molecular species and depend only on the ionization potential and the laser field. According to that model, the gross features of the harmonic radiation are dominated by the laser pulse and by phase matching during propagation. Measurements of intensities, phases and time-structure of harmonic radiation generated in atomic gases were found to be in close agreement with predictions of the three-step model [4]- [7]. The details of the harmonic radiation are determined by structure and dynamics of the gas atoms or molecules. Effects of nuclear dynamics on the harmonic spectra were observed in a recent experiment [8]. To gain information about the system's electronic structure the technique of 'molecular orbital tomography' was proposed and first demonstrated in [9] for the N 2 molecule. Recently, the possibility of applying the method to the slightly more complex CO 2 molecule was demonstrated [10]. The idea behind the tomography experiments is that the main difference between the harmonic response of different gases arises from the recombination step of the three-step model, whose quantum amplitude is given by the bound-continuum dipole matrix element. By comparing the harmonic response from a molecular gas with unknown dipole matrix element to a known (atomic) reference, the matrix element and with it the valence electronic structure of the gas molecule may be reconstructed. The admissibility of this idea and possible corrections to it are discussed in 3 several recent publications [11]- [13]. In this whole picture, it is usually assumed that only a single electron is actually moving in the laser field and that all other electrons can be considered passive spectators in the process. It was pointed out in [14,15] that even for a static core the contribution of core electrons to the dipole matrix element leads to sizable corrections of the reconstructed electron wavefunction.
In all papers discussed so far, the electronic structure is considered to remain static without significant polarization by the laser field. In the present paper, we go beyond this assumption and investigate the effect of polarization of the ionic core by the driving laser field on the harmonic emission spectrum. We solve the time-dependent Schrödinger equation (TDSE) for diatomic model molecules with 2 and 4 active electrons using the multi-configuration time-dependent Hartree-Fock (MCTDHF) method [16,17] and compare results with several single-electron models of increasing complexity. We will show that core polarization profoundly changes the harmonic response thus excluding any single-electron based description of the details of the harmonics. Conversely, dynamical multi-electron effects must be included for the reconstruction of electronic structure from the harmonic response.

Multi-electron molecules in a laser field
A multi-electron system in a laser field is characterized by the time-dependent Hamiltonian Here, the number of electrons is denoted by f . The electron-field coupling is described in dipole approximation and in velocity gauge with the laser vector potential A(t). The nuclei are assumed to be fixed at their positions and the effective molecular potential V phenomenologically describes the interaction of the electrons with the nuclei and a frozen core of inner shell electrons. The repulsive Coulomb potential between the f 'active' electrons is denoted by V ee . The complex (in fact, purely imaginary) absorbing potential V cap is introduced for computational purposes in order to damp reflections at the boundaries of the finite simulation volume. Here and below we use atomic units e 2 =h = m e = 1 unless indicated otherwise.

MCTDHF
The time-dependent Schrödinger equation with the Hamiltonian (1) forms a (1 + 3 f )dimensional partial differential equation, whose problem size increases exponentially with f . Consequently, with present computing capabilities, its direct numerical solution is limited to a maximum of f = 2 active electrons [18]. To extend computations to larger electron numbers, we have developed the MCTDHF method, where the time-dependent multi-electron Schrödinger equation is solved in a restricted function space. The method and our implementation were reviewed in [16] and applications to systems with up to six active electrons were demonstrated in [17]. Here, we only give a short summary.

4
The MCTDHF method consists of making an ansatz for the f -electron wavefunction in the form Spatial and spin coordinates are denoted by q i = ( r i , m i ). Note that both the coefficients C j 1 ,..., j f and the orbitals φ j depend on time. The antisymmetrization of the wavefunction is ensured by the antisymmetry of the C j 1 ,..., j f in all their indices, which allows nonzero C j 1 ,..., j f 's only when all f subscripts j i differ. The summation over all permutations for a given set of { j 1 , . . . , j f } combines the orbitals into Slater determinants. Thus, the MCTDHF wavefunction consists of a linear combination of all Slater determinants which can be constructed from n f singleparticle spin orbitals φ j , j = 1, . . . , n.
For n = f we recover the (single-configuration) timedependent Hartree-Fock (TDHF) approximation. The time-evolution equations for the coefficients C j 1 ,..., j f (t) and the orbitals φ j (t) are obtained from the variational formulation of the time-dependent Schrödinger equation The resulting equations of motion form a system of coupled nonlinear partial differential equations in three spatial dimensions, which we solve numerically using a finite element discretization for the space coordinates and a self-adaptive variable-order Runge-Kutta method for time integration. It is an important feature of the MCTDHF method that the set of orbitals φ j (t) is not fixed, but dynamically adapted to provide an approximation to the correct wavefunction which is optimal in a variational sense and unique up to unitary transformations among the orbitals. Other than the spatial discretization of the orbitals, the only free parameter is the number of orbitals n. Because of the variational optimality of the orbitals, we can use the stabilization of a given observable with increasing n as a convergence criterion.
In the limit n → ∞ the exact time-dependent Schrödinger equation is recovered, while the limiting case n = f is the usual TDHF method. If we define the TDHF wavefunction as 'uncorrelated', MCTDHF allows to systematically include correlation into the approximation by using an increasing number of orbitals n. Although the ansatz (3) is best suited for weakly correlated systems, where the number of Slater determinants remains moderate, in principle the wavefunction converges to the exact solution in the limit n → ∞.

Length versus velocity form of the dipole operator
The harmonic response of a system of electrons to a laser field is proportional to the acceleration of the system's dipole, which can be written either in velocity form d dt or length form The equivalence of the two forms relies on the standard quantum mechanical commutation relations and on the fact that the wavefunctions are exact solutions of a standard Schrödinger equation with a hermitian Hamiltonian. For the computations and models used here, this condition is violated in several respects: (i) the discrete approximation of differential operators does not exactly reproduce the commutation relations, (ii) the complex absorbing potential V cap violates hermiticity, and (iii) the approximate wavefunctions of some of our simplifying models do not solve a standard Schrödinger equation consisting of a Laplacian representing kinetic energy plus potentials in the form of multiplication operators (see sections 4.1 and 4.3). It turns out that in our case absorption at the boundaries and discretization errors only lead to minor differences in spectra obtained in the length and velocity forms: when we use the MCTDHF wavefunctions, spectra calculated in either form-although differing by a factor of up to 10 for the lowest 20 harmonics-agree within 5% for all higher harmonics. The discrepancies at lower harmonics can be ascribed to violation of the hermiticity by the absorbing potential. However, for our model systems with non-standard Hamiltonians rather large discrepancies arise throughout the whole spectrum.
The choice of the dipole operator depends on which properties of the true solution are best reproduced by the approximate solution. The length form depends only on the electron density and emphasizes contributions from larger distances. The velocity form depends on the momentum distribution of the electrons irrespective of their location in space. This favors the length form of the operator for the strong field approximation (SFA), which becomes more accurate at larger distances, where the influence of the molecular potential is weaker. Similar arguments hold for the other models used in the present paper. For this reason and to allow a consistent comparison between models, here we always show the harmonic spectra calculated in the length form.
In [19] a third point of view was taken, namely that the acceleration operator ∇V should be considered as fundamental. The acceleration operator more strongly depends on regions of the electronic wavefunction near the singularities of the nuclear potential, which are dominated by multi-electron effects and where for HHG the dynamics of the electrons is important. While our MCTDHF calculations do include these effects, all models using only a single active electron by definition exclude them. For this reason again we prefer the length form over the acceleration form for the present study.

Model molecules
We choose homonuclear, diatomic molecules with an internuclear distance of | R| = 2.8 for our calculations. Only a few outer electrons are simulated explicitly, while inner shell electrons are accounted for by a screened nuclear Coulomb potential of the form with the potential for one nucleus being The charge Z on each nucleus is equal to half the number of active electrons. The screening parameters α and β have been adjusted such as to yield the ionization potential I p = 0.58 (the same as N 2 ). For homonuclear diatomic molecules, the electronic orbitals can be characterized by their symmetry properties with respect to inversion about the origin: the lowest orbital is of gerade, and the next-higher orbital is of ungerade type. We studied systems with f = 2 and 4 active electrons. For the 2-electron system, in the singlet spin state both electrons occupy a gerade orbital (we abbreviate this system by 2e-gg), while in the triplet spin state, gerade and ungerade orbital are occupied by one electron each (2e-gu system). For the 4-electron system, two electrons occupy the gerade orbital, and the other two the ungerade orbital (4e system). For comparison, we have also calculated harmonic spectra for a 2-electron atom with the same I p as for the molecules, the same general form of nuclear potential, and the same laser parameters.
In table 1, we list the Hartree-Fock (HF) orbital energies for our model systems. To the extent that the HF energies characterize the molecule's electronic structure it is worth noting that the energy spacing between the u and g levels of about E ≈ 0.2 is close to the value for CO 2 with E = 0.21, while being in between N 2 ( E = 0.06) and O 2 ( E = 0.28). Thus, while not meant to represent a specific molecular species, the characteristics of our systems are clearly in the typical range for small molecules.
The laser pulse has linear polarization parallel to the internuclear axis, a carrier wavelength of 800 nm (ω = 0.057), and a cos 4 -shaped envelope for the vector potential A(t) =R A(t) with a FWHM of 1 optical cycle: The rather short pulse duration and the pulse shape are dictated by computational requirements. The cos 4 envelope has finite support and a continuous third derivative everywhere. The latter is needed for calculating the dipole acceleration in the Lewenstein model. The electric field peaks at t = 0 (see figure 1). The laser intensity was I = 3 × 10 14 W cm −2 , which causes between 6 and 18% ionization (see table 2). In our simulations we use cylindrical coordinates (z, ρ) and neglect the dependence on the azimuthal angle. The dimensions of the simulation box ranged from 200 to 300 au in the z-direction (polarization direction), and from 40 to 60 au in ρ, with between 2.5 and 5 discretization points per au. Figure 2 shows the harmonic spectra for our model systems obtained with correlated MCTDHF wavefunctions using up to n = f + 4 orbitals. The overall shape of the spectra is similar in all cases. Due to the extremely short pulse duration there is a double plateau structure. The more intense plateau at lower harmonic numbers corresponds to ionization by the peak of the field with lower acceleration by the subsequent smaller field half-cycle (cf figure 1). The lower intensity plateau reaching to harmonic numbers of ∼ 45 originates from ionization one halfcycle before peak field strength.

Results
We verified the convergence of our calculations with respect to box size, discretization and number of orbitals n.  and 300 au has considerable effects only on the harmonic orders 20. The dependence of our observables on the number of orbitals n is of physical interest as it provides information about the degree of correlation in the molecule, i.e. the number of different Slater determinants needed for an accurate description of the multi-electron wavefunction.
In figure 2, we also included the spectra obtained from uncorrelated calculations (n = f , TDHF, blue) and weakly correlated calculations with n = f + 2 (green). The agreement between n = f + 2 and n = f + 4 for the 2e-gg molecule and 2e atom shows that the spectra can be considered as converged for all practical purposes. Convergence of the 4e spectrum has not been reached, but the main qualitative features at high harmonic orders are preserved.
The rather good agreement of harmonic spectra does not necessarily imply that the TDHF wavefunction reflects other important properties of the true, correlated wavefunction. Figure 3 shows the sum of ionization and excitation yields as a function of n. While for the 2e-gg molecule, correlation leads to a strong increase in ionization plus excitation, for the 2e atom and the 4e molecule ionization remains almost constant. This is consistent with the greater importance of correlation for the ionization of smaller molecular systems reported in [17].

Comparison with the Lewenstein model
The Lewenstein model of HHG [3] is based on the SFA [20]. It disregards excited states and therefore does not account for bound-state dynamics of the molecule in the laser field. In its original form it also neglects any influence of the binding potential on the scattering states. This amounts to representing the continuum by plane wave Volkov states. For a linearly polarized laser field with polarization direction z the expectation value of the dipole is given by where and Here E z (t) = −d/dt A(t) is the electric field strength, I p the ionization potential. The dipole transition matrix element d( k) is taken between the initial electronic ground state |0 and plane waves with the wave vector k. We perform integrations over the wave vectors k of the continuum electrons using the stationary phase approximation as in [3], but the time-integration is performed numerically. for our model systems in comparison with MCTDHF. Only the overall plateau structure and cutoffs, but neither harmonic intensities nor the spectral shape, are reproduced correctly for the molecular systems. By contrast, for the 2e atom good agreement in both aspects is achieved.
In a multi-particle system, the choice of the initial single-electron state |0 is not unique. Possible choices are the highest occupied molecular orbital of a HF wavefunction or the Dyson orbital (see section 4.3 for a discussion of the Dyson orbital). This choice will result in different dipole matrix elements d( k) in the Lewenstein formula. As was shown in [14], the difference between the dipole elements of HF-and Dyson-orbitals is only of second order in multi-particle perturbation theory. Although a Dyson orbital drawn from a correlated calculation changes the Lewenstein spectrum of the 2e-gg molecule near the cutoff, it in no way leads to better agreement with the MCTDHF result. For the 2e atom and 4e molecule, there is no difference between Dyson and HF orbitals. We therefore omit the HF curve in figure 4.

Accounting for the molecular potential
One of the known shortcomings of the Lewenstein model is the absence of the molecular binding potential: neither dynamics of excited states nor scattering from the molecular potential are included. This affects both the ionization and the recombination step in HHG. A variety of corrections for both effects has been proposed [11,21,22]. In order to assess the importance of these effects for our systems, we construct a single-electron Hamiltonian whose initial state coincides with the Dyson orbital, but which includes excited states and scattering from the molecular potential. It is based on the assumptions that a single active electron moves in the potential of a frozen ionic core and the Stark shift of the ground state is negligible. While the second of these assumptions seems quite reasonable, our results below will show that the ionic core cannot be considered as untouched by the laser field.  We start from a field-free single-electron Hamiltonian of the form which contains a molecular potential with screened point charges. We then choose the desired single-electron initial state orbital |0 and divide the Hilbert space into |0 and its orthogonal complement. The effective Hamiltonian H eff acts like H m on the orthogonal complement of |0 and gives the desired ground state energy E 0 for |0 : The interaction with the field has been included here in length gauge, as the field-free molecular Hamiltonian H m should not contain time-dependent terms. As the ground state energy we choose −I p of the multi-electron system. For |0 we use the Dyson orbital, but very similar results are obtained with the highest occupied HF orbital. The one-electron TDSE can be easily solved for this model, which we denote by 1e-TDSE. Figure 4 also includes the 1e-TDSE harmonic spectra (red) obtained with the same orbital |0 as in the Lewenstein model (blue). We immediately see that the dynamics of bound and/or continuum states does play an important role in molecules. Magnitude, overall decrease with harmonic order, and the positions of the minima and cutoff in the spectrum all differ between the two models, although the variations are comparatively smooth. Again, however, like the Lewenstein model, the 1e-TDSE model works well for the high harmonics of the 2e-atomic system.
The 1e-TDSE model describes the ionization process reasonably well. Table 2 compares the sum of ionization plus excitation of TDHF, MCTDHF and 1e-TDSE for the 2 and 4 electron system. The agreement of the model with the TDHF calculation is better than 1%. The use of a Dyson orbital from correlated ground state wavefunctions instead of from HF wavefunctions for |0 in (14) has little impact on the figures. However, in the multi-electron calculation of the 2e-gg molecule, correlation leads to a significantly enhanced ionization and excitation probability, which is not reproduced by 1e-TDSE. Comparing with the MCTDHF spectra of figure 4 we see that in spite of the inclusion of the molecular potential and a good match with the outer electron's orbital, the 1e-TDSE model completely fails to reproduce the molecular high harmonic spectra in intensity as well as spectral shape.
One possible source for the discrepancy between MCTDHF and 1e-TDSE is ionization. We have seen that 1e-TDSE gives overall ionization yields only to within factors 3 of the MCTDHF result, when correlation is included. However, even much larger correction factors to the overall ionization yield could not eliminate the difference in the spectral shapes. According to the three-step model, one important factor determining the spectral shape is the time-structure of the re-collision current. In [13] we have introduced a method to represent the current of recolliding electrons as a function of time. For that the wavefunction is multiplied by a Gaussian mask function centered at some distance z = z 0 from the molecule and the result is analyzed in terms of intensity and momenta. The method is discussed in detail in [13]. Figure 5 shows the momentum distribution of the electron currents at a distance z 0 = 15 with a FWHM of the Gaussian of 8.24 for different times during the pulse as obtained from the MCTDHF and 1e-TDSE wavefunctions. After multiplication by overall normalization constants between 0.40 and 0.90, which quite well reflect the errors in total ionization for the various systems, electron momentum distributions at different times agree within 10% or better. We conclude that the time-distribution of ionization is well described by the 1e-TDSE, though the magnitude is not. When we reason within the three-step model, this leaves only the recombination step as the source of error in the 1e-TDSE model.

Multi-electron corrections to the SFA
One important modification of the SFA has been introduced in [14,15], where multi-electron effects were included in the recombination matrix element. In this case, the field-free bound state |0 is the multi-electron ground state | 0 itself, and the continuum states are written as products of the (field free) ionic ground state | + and single-electron plane waves | k . The total wavefunction for a spin singlet becomes where j=1 P j f ) is the antisymmetrization operator, f is the number of electrons and ↑, ↓ denote the ±1/2 spin states. The dipole matrix element is given by which can be evaluated further as [15] where φ D denotes the Dyson orbital.  Figure 6 shows the importance of these corrections for our molecules. Dipole matrix elements were calculated with the respective neutral and ionic MCTDHF wavefunctions. We see that the exchange correction can be very substantial, amounting to more than one order of magnitude in the case of the 4e system.
The exception is the 2e-gg system, for the simple reason that the exchange term d 2 nearly vanishes. It vanishes exactly in the HF case since both the HF orbital φ and the ionic single-electron orbital φ + have gerade symmetry (z 1 refers to the z-coordinate of the first electronic orbital). The extra multi-configuration corrections remain small. The same remarks apply to the 2e atom. Unfortunately, it is evident from comparing figures 4 and 6 that also the exchange-corrected Lewenstein model fails to reproduce the multi-electron solution, and it still does so by orders of magnitude. A failure with respect to overall harmonic intensity alone might have been ascribed to one of the known shortcomings of the Lewenstein model, namely incorrect ionization rates. What is more disturbing and can certainly not be explained in this way is the strong deviation in the spectral shape at high harmonic energies, as molecular orbital tomography [9] depends on the unambiguous understanding of the harmonic power spectrum (and phases).

Factoring the multi-electron wavefunction
In view of the failure of the previous models, we now attempt to directly extract a timedependent single-electron function from the known multi-electron function such that the harmonic spectra agree for the two functions.
In the single-electron description one assumes that the movement of the electron that becomes ionized is on a different spatial and velocity scale from the electrons remaining closer to the nucleus. It may therefore be possible to write the total wavefunction as a product 14 of a core (x 1 , . . . , x f −1 ; t) of electrons times a single-particle orbital φ(x f ; t) with proper antisymmetrization: Note that, since and φ are, in general, not orthogonal, we have to normalize by a timedependent factor.
For the determination of the constituents and φ two requirements should be met: (i) the product function should well approximate the true solution and (ii) the inactive electron part should contribute little to our observable (the harmonic spectrum). To meet the first condition, φ is determined such that, given a particular , provides the best approximation to the true . Here, we define the best approximation as the one with the largest overlap | . Apart from this, the specific flavor of a particular model is introduced by . Its choice depends on the physical process under consideration. For HHG a suitable choice seems to be the cationic wavefunction. In this case, φ is the Dyson orbital.
Technically, the Dyson orbital [23] is defined by the integral over the ionic coordinates for which we use the short notation The norm of the Dyson orbital φ D (the 'pole strength'), is not necessarily equal to 1. Rather, the deviation from unity is a measure of how much the orbitals of the ion core change upon ionization. If it is impossible to satisfy requirement (i) with the ansatz (22), then the system is significantly correlated. We consider two possible choices of core electron factor : 1. The ion remains in its initial state (static), i.e. is not significantly deformed by the laser pulse, 2. The ion is considerably polarized by the laser field and we need to include its correct time evolution: where + (t) is the time-dependent ionic wavefunction. Figure 7 (left column) shows spectra obtained with the static and polarized ionic core. We see that polarization of the ionic core strongly affects the overall shape and intensity of the harmonic spectrum. As we have shown above, multi-electron corrections strongly affect the harmonic spectrum. It is therefore to be expected that even minor modifications of the inner shell electrons may strongly change the overall spectrum. Note that this change is not due to harmonic emission by the ion itself, which is by about two orders of magnitude lower than emission from the neutral (figure 7, right column). The differences between static and polarized core factorizations are on a similar scale as the differences between the 1e-TDSE models and the MCTDHF result. It is intriguing that for all our models the high harmonic part of MCTDHF spectra is well approximated by the polarized core factorization. Interestingly,  Figure 7. Comparison of harmonic spectra obtained by MCTDHF and from product wavefunctions as (22). The product wavefunctions with the static or dynamic ion core are denoted by stat and dyn , respectively. In the left column (a, c, e), stat (green) and dyn (magenta) are compared to the MCTDHF (black) solution. In the right column (b, d, f), we compare the Dyson orbital φ dyn D (cyan) with the product wavefunction dyn . For reference, the spectrum of the ion is also shown (gray). for the 2e atom the factorized wavefunction leads to no improvement over the single-electron models (Lewenstein and 1e-TDSE), which actually yield better agreement in terms of harmonic intensity. Figure 7 (right column) shows, that the 'dynamical Dyson orbital' φ dyn D dominates the harmonic emission: in all cases, including the exchange corrections between φ dyn D and + (t) only slightly changes the harmonic spectrum. The spectra for the Dyson orbital are obtained from the dipole expectation value φ dyn D |z|φ dyn D .

Accuracy and comparability of the calculations
As remarked in the discussion of the MCTDHF results (section 3) the convergence of spectra with respect to the number of orbitals (figure 2) is not fully reached numerically, with changes remaining in the order of a factor 2 in the worst case (4e molecule). Even so, this is considerably less than the observed deviations between MCTDHF and simpler models, which reach up to three orders of magnitude.
The importance of core-polarization effects derives from the in general strong contribution of the cross-matrix elements between the 'active' orbital and the core orbitals. This observation was originally made for the multi-electron corrected SFA [14,15] and is corroborated here for the case when a frozen core is factorized out from the full MCTDHF wavefunction. The same procedure with the frozen core replaced with the polarized core systematically gives largely correct HHG spectra for all our model systems. In that case contributions from cross-matrix elements between the active single-electron factor and the core factor are small.
The success of the polarized core factorization, together with the failure of the static core factorization, also provides indirect evidence that both the ionic and the neutral wavefunctions are well described by MCTDHF for the purpose of the calculation of high harmonic spectra. All other features like frozen inner shell electrons, cylindrical symmetry, or fixed nuclei are common to all models. Therefore the comparison with MCTDHF allowed us to disentangle specific shortcomings of the single-electron models from the true multi-electron effects in HHG. The magnitude of these effects is linked to our specific molecular models. However, with key model parameters chosen in realistic ranges and calculations made in three-dimensions (3D), we expect effects on the same scale for real molecules.

Conclusions
We found strong indications that HHG is a genuinely dynamic multi-electron effect. Neither incorrect ionization rates nor the plane wave approximation, which could be remedied by more refined single-electron models, are responsible for the failure of SFA for our multi-electron model systems. We also conclude that in the present case the multi-electron effects cannot be described in the spirit of [14,15] as the interaction of some frozen core with an active electron.
For these realistic 3D multi-electron molecular models, both the standard SFA and an advanced single-electron model that includes a phenomenological molecular potential completely fail to reproduce shape and intensity of the harmonic spectra. We could show that for our system parameters the ionization process is described with satisfactory fidelity by the single-electron model, but recombination is grossly incorrect. The discrepancy is due to the large size of the multi-electron contributions to the recombination matrix element [15] and their sensitivity to the polarization of the ionic core by the laser pulse. A largely correct description of the higher part of the harmonic spectrum can be obtained, when the full wavefunction is factorized into the polarized ion wavefunction times a single electron factor. In the case of such a factorization, the cross-terms between the polarized ion and the 'active' electron contribute little and the single-electron factor already produces a roughly correct spectrum. In this very general sense, one might argue that HHG from molecules is still a single electron process. However, the time-evolution of this factor wavefunction implicitly includes interactions with the polarized core-electron wavefunction. In particular, any dynamical description of its time-evolution must include terms accounting for the time-dependent exchange and correlation.
At present, due to computational limitations, we cannot conduct a similar calculation ab initio for systems most frequently used in experiments like N 2 or CO 2 . Even within our model systems confidence in the correctness of the results is based on an analysis in terms of various models rather than systematic convergence studies, as these would require larger computer resources than presently available to us. Still, we believe that the rather realistic nature of our MCTDHF simulations raises serious questions about the possibility of tomographic imaging for general systems that must be decided before the method can be put to practical use.