Influence of large permanent dipoles on molecular orbital tomography

The influence of large permanent dipoles on molecular orbital tomography via high-order harmonic generation (HHG) is investigated in this work. It is found that, owing to the modification of the angle-dependent ionization rate resulting from the Stark shift, the one-side-recollision condition for the tomographic imaging can not be satisfied even with the few-cycle driving pulses. To overcome this problem, we employ a tailored driving pulse by adding a weak low-frequency pulse to the few-cycle laser pulse to control the HHG process and the recollision of the continuum electrons are effectively restricted to only one side of the core. Then we carried out the orbital reconstruction in both the length and velocity forms. The results show that, the orbital structure can only be successfully reproduced by using the dipole matrix elements projected perpendicular to the permanent dipole in both forms.


Introduction
The recent development in strong-field physics has opened wide perspectives for observing the structure and ultrafast electron dynamics in atoms and molecules withÅngström and attosecond resolutions [1][2][3][4][5]. One of the most fascinating topics at the frontier of this area is the molecular orbital tomography (MOT) based on high-order harmonic generation (HHG), which was first proposed by J. Itatani et al in [6]. By applying this method to a series of aligned molecules, they accomplished a tomographic reconstruction of the highest occupied molecular orbital (HOMO) of N 2 . Since then, the MOT has attracted a great deal of attention for its great significance in uncovering the basic properties of molecules and chemical reactions [7][8][9][10][11][12][13]. In most of the works, the target molecular orbitals are symmetric with either gerade or ungerade symmetry, while the MOT of nonsymmetric orbitals is more difficult. It has been discussed by Van der Zwan et al [14] that the MOT of orbitals with arbitrary symmetry requires one to restrict the recollision of the continuum electrons to only one side of the core in the HHG process [15,16]. Therefore, they suggested using extremely short few-cycle laser pulses to ensure the oneside-recollision and the molecular orbital of HeH 2+ was successfully reconstructed in their simulation. Besides, one can also consider using an ω + 2ω two-color laser pulse for the MOT of nonsymmetric orbitals [7,17].
In addition to the arbitrary symmetry, another important characteristic of the nonsymmetric orbital is the permanent dipole. As a large permanent dipole would lead to an unnegligible Stark shift of the orbital energy, it would have influence on the MOT in two aspects. Firstly, it will cause additional laser-induce Stark phase acquired in the HHG of the nonsymmetric orbital. As the phases of field-free dipole matrix elements are needed in order to reconstruct molecular orbitals, one has to find a way to subtract the Stark phase from experimental measurements. The method is found in [18], in which the authors derived an analytical expression of the firstorder Stark phase only depending on the measurable quantities. Secondly, previous works have shown that the Stark shift of the orbital energy can efficiently modify the angular dependence of the ionization rate and the distribution of the ionized electrons [19][20][21]. This effect will make it more complicated to control the continuum electron for one-side-recollision and to reconstruct the orbital. This influence has however not been discussed before to the best of our knowledge and requires more concerns. Moreover, the influence of the nonsymmetric electron distribution on the reconstructed result in MOT has not been extensively discussed and calls for more studies too.
In this paper, the influence of large permanent dipoles on the MOT of nonsymmetric molecular orbital is investigated. It is found that, the one-side-recollision condition in HHG can not be satisfied even with the few-cycle laser pulses, owing to the modification of the angle-dependent ionization rate induced by the Stark shift. We employ a tailored laser pulse to overcome this problem and the orbital reconstruction is carried out in both the length and velocity forms. The reconstructed results are compared to the ab initio orbital and the influence of the nonsymmetric electron distribution on the reconstructed result in MOT is discussed at the end.

Theoretical model
This section includes two independent parts: "Simulating the HHG" and "Reconstruction method". In the first part, we describe the modified Lewenstein model in the length form, which is employed in this work to simulate the HHG. In the second part, we introduce the orbital reconstruction method in both the length and velocity forms. Note that in the second part the inputting high-order harmonic spectra for the reconstruction process could be obtained from any ways: from calculations of different theories and forms or from experimental measurements.

Simulating the HHG
According to the Lewenstein model [16] and using the approximate saddle points [22,23], the time-dependent dipole moment is calculated by (atomic units are used throughout) where a ion , a prop and a rec are the amplitudes of tunneling ionization, propagation after tunneling and recombination. t d is the approximate saddle time determined by solving A f (t) is the vector potential corresponding to the electric field F(t), and p st is the stationary momentum calculated by For atoms the amplitudes in Eq. (1) are given by a rec (k) = < Ψ g |r|e ik·r > .
I p is the ionization energy and Ψ g is the ground state in the form of hydrogenlike atom [16].
is the momentum of the return electron at the instant of recombination. For molecules, the ionization amplitude a ion is calculated with the MO-ADK theory [24], which takes into account the structure features of molecular orbitals. Ψ g is the HOMO of the molecule obtained from Gaussian 03 ab initio code [25]. For polar molecules whose HOMO has an unnegligible permanent dipole µ h , the ionization energy becomes time-dependent in response to the external electric field [20,21,26]: where I p0 is the field free ionization energy. The permanent dipole µ h is calculated by ρ H (r) is the electron density of the HOMO calculated by ρ H (r) = |Ψ g (r)| 2 . For the Lewenstein model in the length form, a modification should be applied in Eq. (6) that the momentum k is replaced by the effective momentum [27,28] This modification implies that the recombination occurs within the potential well and the acceleration effect of the binding potential is considered in this modified Lewenstein model [29,30]. In Eq. (9), the factor γ(k) is given by .
This k-dependent factor avoids the unreasonable fact that all the return electrons have momenta larger than 2I p . The introduction of γ(k) also suppresses the strong background noises resulting from the jump of k e f f over the gap from − 2I p to 2I p . As γ = 1 for k 2 ≥ I p , this term only affects the HHG of very low orders. Finally, the harmonic spectrum is obtained by Fourier transforming the time-dependent dipole acceleration a(t): where a(t) =D(t), ω L is the frequency of the driving pulse and q corresponds to the harmonic order.

Reconstruction method
As shown in Fig. 1, (x,y,z) is the molecular frame and (x',y',z') is the laboratory frame. The driving laser field is linearly polarized along the x' axis and propagates along the z' axis. To reconstruct the molecular orbital, one must calculate or measure the high-order harmonic spectra (including the intensity and phase) from the oriented target molecule at various orientation angles θ . In addition, one should also calculated or measure the high-order harmonic spectrum from the reference atom with the ionization energy similar to that of the target molecule. Then the recombination dipole moment d L and the dipole velocity d V for the molecular orbital in both x' and y' components are obtained by [7] where E, A, ϕ denote the complex electric field, amplitude and phase of the high-order harmonic radiation respectively. η is the scaling factor in molecules depending on the ionization angle θ i calculated by the MO-ADK theory. ω = qω L is the harmonic frequency of the q-th order. The relation between ω and k is determined by the energy conversation ω = k 2 /2. To subtract the additional acquired Stark phase from E mol , the first-order Stark phase Φ (13) and (14), which could be approximately evaluated with the analytical expression derived in [18]. In the present paper, to focus our investigation on the influence of the permanent dipole in the aspect of modifying the ionization rate, the MOT is performed with the Stark phase perfectly subtracted in the simulation. In addition, the result of MOT with the Stark phase subtracted following [18] as well as the influence of the permanent dipole in phase are briefly discussed at the end of Section 3.
The dipole matrix elements of the "known" reference atom can be obtained by whereΨ re f (k) is the Fourier transform of the ground state of the reference atom Ψ re f (r). By using the dipole moment or the dipole velocity, the reconstruction of the molecular orbital can be performed in the length form or in the velocity form.
In the length form, the orbital can be reconstructed as: where Ψ L x (x, y) is the reconstructed orbital by using the dipole moments projected on the x axis of the molecular frame and Ψ L y (x, y) is the reconstructed orbital by using the dipole moments projected on the y axis. For symmetric molecules, the same orbital can in principle be reconstructed from both the x and y components of the dipole moments (d L x and d L y ). However, due to the limited discrete sampling in Fourier space, they will most likely not give the same result [7]. Generally, the final reconstructed orbital is defined as to hopefully average out the distortions.
Similarly, one can also reconstruct the orbital in the velocity form using the two projected components of the dipole velocity respectively: with

Result and discussion
This section consists of two parts. In the first part, we analyze the HHG of a polar molecule to demonstrate the influence of the large permanent dipole on MOT and suggest a method to overcome this unfavorable influence. In the second part, we compare and discuss the reconstructed results in different forms and by using different projected components of the dipole matrix elements.

Simulating and controlling the HHG
In this work, we apply CO as the target molecule. The HOMO of CO has a large permanent dipole of 1.72 a.u. pointing from C to O nucleus. The ionization energy of CO is 0.5152 a.u., and the reference atom is Kr with the ionization energy of 0.5145 a.u. We apply a linearly polarized three-cycle sin 2 driving pulse where T 0 is the optical cycle and ϕ 0 = 1.25π. This form of laser pulse (with an opposite sign) has been used to accomplish the MOT of the nonsymmetric molecule HeH 2+ [14]. The intensity and wavelength of the laser pulse we use are 1 × 10 14 W/cm 2 and 1.5 µm respectively. Employing the 1.5 µm mid-IR laser source allows us to extend the harmonic spectrum considerably with a low laser intensity below the barrier suppression intensity [12,31], which ensures the validity of the used MO-ADK model and strong-field approximation (SFA) in our theoretical description. Another benefit of keeping the intensity low is that it will reduce the effect of the Stark shift.
The electric field of the driving pulse is depicted by the dashed blue curve in Fig. 2(a). The HHG process is also illustrated by the black arrows in Fig. 2(a): under this three-cycle laser pulse, the HHG is dominantly contributed by the electrons ionized by the two peaks P1 and P2. For the electrons ionized by P1, they first drift towards the positive direction and finally return to the core at times around T1 with negative momenta, emitting the high-order harmonic photons. On the contrary, the electrons ionized by P2 finally return to the core at times around T2 with positive momenta.
When the effect of the permanent dipole of CO is ignored, the ionization energy is constant. The ionization probabilities in the +x' and -x' directions are dominantly determined by the waveform of the electric field, due to the exponential dependence of the tunneling ionization rate on the amplitude of electric field. Therefore, as the maximum amplitude of P2 is bigger than that of P1, electrons are more easily to be ionized by the P2 peak and the high-order harmonic radiation will be stronger at T2 than at T1. This is confirmed by our calculation shown in the first row in Fig. (3). In Fig. 3(a), the red and green solid curves show the x' and y' components of the high-order harmonic spectrum from CO oriented at θ = 10 • without taking into account the Stark shift. The harmonic spectrum from the reference atom is also shown in Fig. 3(a) by the dashed curve. Figures 3(b) and (c) show the corresponding time-frequency analyses [26,32] of the molecular HHG in the parallel (x') and perpendicular (y') directions respectively. It is found that the atomic spectrum exhibits double-plateau structure, while only the plateaus on the lower energy side are observed in the molecular spectra. This is because the plateaus on the higher energy side are so weak as to be below the background noise. The result implies that the plateaus on the lower energy side are at least 3-4 orders of magnitude higher than those on the higher energy side in both x' and y' components in the spectrum from CO. According to the time-frequency analyses, the plateaus on the lower energy side are mainly contributed by the Q2 peaks, which are generated by the recollision of electrons with positive momenta at T2. While the plateaus on the higher energy side are contributed by the Q1 peaks, which are generated by the recollision of electrons with negative momenta at T1. From both the observations in the spectra ( Fig. 3(a)) that the plateaus on the lower energy side are much higher than those on the higher energy side and in the time-frequency distributions (Figs. 3(b) and (c)) that Q2 are much stronger than Q1, it can be concluded that the high-order harmonics are generated by the recollision of electrons from only one side.
However, when the stark shift is considered, the one-side-recollision condition is worse satisfied. As shown in Figs. 3(d), the plateaus on the higher energy side appear and are in the same level with the plateaus on the lower energy side in the harmonic spectra of CO in both components. The time-frequency analyses also show that the peaks Q1 and Q2 have comparable intensities. It has been discussed that Q1 is generated by the electrons with negative return momenta and Q2 is generated by the electrons with positive momenta. Therefore, the results indicate that the electrons with both negative and positive momenta have comparable contributions to the HHG. To explain why the Stark shift makes the one-side-recollision condition worse satisfied, the time-dependent orbital energy E g = −I p in response to the external electric field F(t) is plotted in Fig 2(b). The electron is more likely to be ionized with a high E g and is harder to be ionized when E g is low. Therefore, the ionization rate of CO is decreased when µ h is parallel with the electric filed F and is increased in the antiparallel geometry. In the case of 10 • orientation, as shown by the red solid curve in Fig. 2(b), the ionization by P1 is enhanced by the Stark shift and the ionization by P2 is suppressed, leading to the comparable ionization probabilities from both peaks. And these comparable ionization probabilities finally lead to the comparable HHG contributions by electrons returning from opposite sides.
As the MOT requires extracting structure information (the slices) of the target orbital from different directions, the one-side-recollision condition should be satisfied at various orientation angles. To evaluate the ratio between the probabilities that an electron returns to the core with negative and positive momenta at different angles, a semi-classical method introduced in [33] is used. First of all the electrons are tunneling ionized with the probability calculated by the MO-ADK theory. After ionization, the freed electron follows a classical trajectory, and we collect the electrons which finally return to the core. Every classical return receives a weight based on the ionization probability and a factor of τ 3 , where τ is the time the electron spends in the continuum from ionization to recombination and reflects the effect of wave-packet spreading. After the collection, the return probability of the electron with negative momentum P − and the return probability of electron with positive momentum P + are obtained. The ratio of the probabilities P − /P + for orientation angles from 0 • to 180 • is plotted in Fig. 4, where the blue Fig. 4. The ratio between the probabilities that an electron returns to the core with negative and positive momentum for various orientation angles in the cases of: (1) using only the three-cycle driving pulse without taking into account the Stark shift (SS) (blue circles), (2) using only the three-cycle driving pulse with the SS (green diamonds), and (3) using the tailored pulse with the SS (red pentagrams). The inset shows the ratios on a logarithmic scale.
circles and the green diamonds present the calculated results without and with the Stark effect taken into account respectively. Comparing the green diamonds to the blue circles, the Stark shift leads to more recollision with negative momentum especially for the orientation at small angles. The ratio of P − /P + is bigger than 0.1 for θ < 60 • and goes up to bigger than 1 for θ = 0 and 10 • .
To overcome this problem, one may first think of using an even lower intensity for the driving pulse to further reduce the effect of the Stark shift. However, too low intensity will bring the problem of a bad signal to noise ratio. Instead, in previous experimental and theoretical works higher intensities than 1×10 14 W/cm 2 were preferred, e.g. Ref. [12,26]. Therefore, in this work we suggest a method to solve this problem by using a tailored laser pulse with the relatively high intensity. The tailored laser pulse is synthesized by adding a weak low-frequency laser field to the few-cycle laser field. The low-frequency field can be provided by different kinds of light sources such as the CO 2 laser [34, 35] and the THz radiation [36,37]. The basic idea of adding this assistant field is to enlarge the difference between the amplitudes of P1 and P2. In the present work, we apply the low-frequency pulse from CO 2 laser with the wavelength of 10.6 µm. The intensity of this pulse is 1% of that of the few-cycle driving pulse and the electric field is F = F 0 cos(ω L t + ϕ 0 ) with ϕ 0 = 1.7π. The low-frequency field and the tailored field are presented by the dash-dotted green and solid red curves in Fig. 2(a) respectively.
The ratio P − /P + under this tailored field is presented by the red pentagrams in Fig. 4. Comparing to the green diamonds, the recollision with negative momentum is suppressed: the ratio is two orders of magnitude decreased and is below 0.1 for all the orientation angles (for 180 • < θ < 360 • the ratio at θ equals that at 2π − θ ). The result shows that the tailored laser pulse can efficiently restrict the recollision of the continuum electrons to only one side of the Fig. 5. The calculated high-order harmonic spectra (the first column) and the corresponding time-frequency distributions for the x' and y' components of the molecular harmonic spectra (the second and third columns), by using the tailored pulse taking into account the Stark shift. The first, second and third rows correspond to θ = 10 • , 90 • and 170 • respectively. The values in the color bar correspond to 10log 10 (Intensity). core.
The high-order harmonic spectra at three typical orientation angles 10 • , 90 • and 170 • calculated with the quantum model described in Section 2.1 are shown in the first column in Fig. 5. The second and third columns present the time-frequency analyses of the molecular harmonic spectra in the x' and y' components respectively. These results are consistent with those from the semi-classical method. The peaks Q2 are much stronger than Q1 in all the time-frequency distributions. Regarding the molecular harmonic spectra, the plateaus on the lower energy side (predominantly contributed by electrons with positive momenta) are 3-4 orders higher than those on the higher energy side (contributed by electrons with negative momenta) at θ = 10 • as shown in Fig. 5(a). And the plateaus on the lower energy side are at least 4 orders of magnitude higher than those on the higher energy side for θ = 90 and 170 • as shown in Figs. 5(d) and (g). Therefore, the one-side-recollision condition is met and the HHG by this tailored pulse is satisfying for the MOT of the nonsymmetric target orbital.

Tomographic reconstruction
The MOT result in the length form is shown in Fig. (6). To reconstruct the orbital, we use the odd harmonics from the 23rd to the 79th orders (corresponds to 19-65 eV) generated from CO at θ from 0 • to 360 • with the angular step ∆θ = 10 • . Figures 6(a) and (b) display the reconstructed orbitals Ψ L x and Ψ L y using the dipole moments projected on the x and y axes in the molecular frame respectively. The ab initio orbital is also shown in Fig. 7(c). Compared with Fig. 7(c), Fig. 6(b) successfully reproduces the main structure of the nonsymmetric target orbital, while Fig. 6(a) fails to reproduce the structure. There is an artificial structure close to x=0 in Fig. 6(a) with a sharp jump from the negative value to the positive value. This is due to the nonsymmetric distribution of the orbital.
To illustrate this, in Figs. 6(c) and (d) we present the intermediate functions f x and f y , which are defined as [6] The relation between Ψ L and f is Both f x (x, y) and f y (x, y) are symmetrically distributed along the y direction and are nonsymmetrically distributed along the x axis, corresponding to the nonsymmetric electric density distribution of the target orbital along the permanent dipole. Therefore, one will have no problem in dividing f y by y to obtain Ψ L y , because the nodal plane of f y locates right at y = 0. But the nodal plane is bent and not coincident with x = 0 for f x . As a result, numerical problem arises when dividing f x by x with x→0, which leads to the artificial structure and sharp jump in Fig.  6(a). The reconstruction is also performed in the velocity form as shown in Figs. 7(a) and (b). The scan of θ starts from 5 • and increases in the step of 10 • to avoid cos(θ ) = 0 or sin(θ ) = 0 in the denominators in Eqs. (22) and (23). Ψ V x and Ψ V y are the reconstructed orbitals using the dipole velocities projected on the x and y axes respectively. Again, Ψ V x and Ψ V y give different results. Ψ V y successfully reproduces the main structure of the target orbital, while Ψ V x is incorrect. To further examine the influence of the nonsymmetric distribution of the orbital on the reconstructed result, we have also tried the MOT of another kind of orbital which is nonsymmetrically distributed along the y axis. On the contrary to the results shown in Figs. 6 and 7, Ψ L x and Ψ V x successfully reproduce the orbital this time, while Ψ L y and Ψ V y are incorrect (not shown). All the results indicate that, the MOT can only be accomplished by using the dipole matrix elements projected perpendicular to the permanent dipole in both the length and velocity forms for the nonsymmetric orbitals.
Further comparing the reconstructed orbitals Ψ L y and Ψ V y with the ab initio orbital, it is found that there is still detailed structure missed in the reconstructed results near the left core. This is due to the limited spectral range "detected" for the MOT [17]. In Fig. 7(d), we compare the slices along the x axis of Ψ L y , Ψ V y and the ab initio orbital. The comparison shows that, the result obtained in the length form agrees with the ab initio orbital well for the positions and amplitudes of the left and middle peaks, but the position of the right peak deviates far away from the ab initio orbital. The result obtained in the velocity form matches the ab initio orbital better for the positions of all the three peaks, but the values at x=-3 and 0 a.u. deviate a lot from those of the ab initio orbital. Both forms have limitation in reproducing the exact shape of the target orbital. Furthermore, the reconstruction of the orbital in the y direction is not satisfying. This error is also observed in the previous works especially for the σ orbital. Besides the limited spectrum range, the imaging errors may also result from several other reasons, such as the discrete sampling in the frequency domain with the interval of 2ω L , and the SFA and singleactive-electron approximation used in the theoretical model [7][8][9][10]12]. In the MOT process, when the Stark phase is approximately subtracted by using the analytical expression derived in [18], the main structure of the target orbital is also successfully reproduced in both the length and velocity forms. The reconstructed orbitals have similar shapes to those shown in Figs. 6 and 7, except that the whole structures are shifted towards the right in both forms. These shifts may therefore originate from the approximation in evaluating the Stark phase. On the other hand, this result implies that the influence of the permanent dipole in phase would lead to the shift of the reconstructed orbital, if the Stark phase is not exactly subtracted.

Conclusion
The influence of large permanent dipoles on the MOT is investigated in this paper. It is found that, owing to the modification of the angle-dependent ionization rate induced by the Stark shift of the orbital energy, the one-side-recollision condition can not be satisfied even with the few-cycle driving pulse. To overcome this problem, we employ a tailored laser pulse by adding a weak low-frequency pulse to the few-cycle pulse to control the HHG process and the requirement for the one-side-recollision is met. Then the reconstruction of the target orbital is performed in both the length and velocity forms. The results show that, the orbital structure can only be successfully reproduced by using the dipole matrix elements projected perpendicular to the permanent dipole in both forms.