Strong-field double ionization dynamics of vibrating HeH$^+$ versus HeT$^+$

We study double ionization (DI) dynamics of vibrating HeH$^+$ versus its isotopic variant HeT$^+$ in strong laser fields numerically. Our simulations show that for both cases, these two electrons in DI prefer to release together along the H(T) side. At the same time, however, the single ionization (SI) is preferred when the first electron escapes along the He side. This potential mechanism is attributed to the interplay of the rescattering of the first electron and the Coulomb induced large ionization time lag. On the other hand, the nuclear motion increases the contributions of these two electrons releasing together along the He side. This effect differentiates DI of HeH$^+$ from HeT$^+$.

Present studies have revealed the importance of tunneling and rescattering in strong-field processes [15,16]. In comparison with other strong-field processes, double and multiple ionization which involves electron-electron correlation includes richer physical phenomena. As the double ionization (DI) from atoms and symmetric molecules have been studied widely [10], DI from polar molecules with a large permanent dipole [17,18] is less studied. Especially, when the nuclear motion is considered, the situation is more complex. It has been shown that for small polar molecules such as HeH + , the interaction of the strong laser field and the permanent dipole induces the rapid nuclear motion which has important influences on HHG and ATI of the asymmetric system [19][20][21][22][23]. In addition, the interplay of the Coulomb effect and the permanent-dipole effect also gives rise to a strong asymmetry in photoelectron momentum distributions (PMD) for single ionization (SI) of HeH + [24]. This asymmetry is closely related to the permanent dipole induced asymmetric ionization [25] and the Coulomb induced large ionization time delay [26]. Effects of these mechanisms on DI of the asymmetric system are not clear so far.
In this paper, we focus on DI of HeH + and its isotopic variant HeT + in strong linearly polarized laser fields beyond the Born-Oppenheimer (BO) approximation. The HeH + system, the simplest polar heteronuclear molecule, has served theoretically and experimentally as a fundamental benchmark system for understanding molecular formation and electron correlation [27]. Numerical solution of the time-dependent Schrödinger equation (TDSE) of the vibrating two-electron system in full dimensions is still not within reach and is easily limited by the existing computing capability. Thus, we use a simplified model where the motion of all particles is restricted to one dimension (1D). It has been shown that such a model can reproduce all important strong-field effects such as multiphoton ionization and HHG [28][29][30]. This simplified model can also describe qualitatively the correlation effects between electrons and interplay between the electronic and the nuclear motion [31][32][33][34].
The calculated PMDs of DI for HeH + or HeT + show a striking asymmetry with indicating that these two electrons in the DI process prefer to release together along the H(T) side. This phenomenon in DI differs remarkably from that in SI for HeH + or HeT + , which shows that the first electron in SI prefers to escape along the He side. This disagreement between DI and SI strongly implies that the rescattering of the first electron plays an important role in DI of the asymmetric system. This rescattering along with the Coulomb induced large ionization time delay remarkably increases the DI yields and results in preferred DI along the H(T) side. On the other hand, the rapid nuclear motion, which differs remarkably for HeH + and HeT + , increases the contributions of direct ionization (which is preferred along the He side) to DI.

II. NUMERICAL METHODS
The Hamiltonian of the asymmetric molecule HeH + studied here has the following form (in atomic untis of h = e = m e = 1): Here R is the internuclear separation and x j (j=1,2) is the electronic coordinate. µ N = M He M H /(M He +M H ) is the nuclear reduced mass and µ e = (M He +M H )/(M He + M H + 1) ≈ 1 is the electronic reduced mass. M He and M H are masses of He and H nuclei. The term V en denotes the interaction between the electron and nuclei and has the following form: where Z 1 = 2 and Z 2 = 1 are the charges for He and H centers, respectively. R 1 and R 2 are positions of He and H nuclei with R 1 = M H R/(M He + M H ) and R 2 = −M He R/(M He + M H ). ǫ = 0.59 is the smoothing parameter, which is adjusted such that the ground-state energy of the model HeH + molecule matches the real one of E 0 = −2.98 a.u.. The equilibrium separation of model HeH + studied here is R e = 2 a.u., which also holds for model HeT + and is somewhat larger than the real one of R e = 1.4 a.u.. Here, we have used the length-gauge form of the interaction Hamiltonian. The laser field used here is E(t) = E 0 f (t) sin(ω 0 t) with peak amplitude E 0 , envelope function f(t) and laser frequency ω 0 . In our simulations, we use a seven-cycle laser pulse which is linearly turned on and off for two optical cycles, and kept at a constant intensity for three additional cycles. We use Ψ(R, x 1 , x 2 , t) ≡ Ψ(t) on a three-dimensional grid to represent the wave function. The TDSE of iΨ(t) = H(t)Ψ(t) is solved numerically using the spectral method [35]. A grid size of L x1 × L x2 = 204.8 × 204.8 a.u. with the grid step of ∆x 1 = ∆x 2 = 0.4 a.u. for the electron, a range of R = 0.6...6.9 a.u. with the grid step of ∆R = 0.1 a.u. for the internuclear distance, and a time step of ∆t = 0.05 a.u. have proven sufficient convergence for describing the strong field dynamics. In order to avoid the reflection of the electron wave packet from the boundary and obtain the momentum space wave function, the coordinate space is split into the inner and the outer regions with Ψ(t) = Ψ in (t) + Ψ out (t), by multiplication using a mask function F (x 1 , x 2 , R) = F 1 (x 1 )F 2 (x 2 )F 3 (R). Here, F 1 (x 1 ) = cos 1/2 [π(|x 1 | − r 0 )/(L x1 − 2r 0 )] for |x 1 | ≥ r 0 and F 1 (x 1 ) = 1 for |x 1 | < r 0 . r 0 = 3/8L x1 is the critical boundary between the inner and the outer regions for one electron. The relevant electronic wave packet passing through this critical boundary will be absorbed by the mask function smoothly. The form of F 2 (x 2 ) is similar to F 1 (x 1 ). A similar absorbing procedure with the mask function F 3 (R) is also used for the upper boundary of R. In the inner region, the wave function Ψ in (t) is propagated with the complete Hamiltonian H(t). In the outer region, the time evolution of the wave function Ψ out (t) is carried out in momentum space with the Hamiltonian of the free electron in the laser field [36][37][38]. The mask function is applied at each time interval of 1 a.u. and the obtained new fractions of the outer wave function at the DI condition of |x 1 | ≥ r b and |x 2 | ≥ r b , denoted with Ψ d out (t), are added coherently or non-coherently to the corresponding momentumspace wave functionΨ d out (t). Finally, we obtain PMDs c d (p 1 , p 2 ) of DI fromΨ d out (t). Similarly, with coherently or non-coherently adding the obtained new fractions of the outer wave function at the SI condition of |x 1 | < r b or |x 2 | < r b , denoted with Ψ s out (t), to the corresponding momentum-space wave functionΨ s out (t), one can obtain PMDs c s (p 1(2) ) of SI fromΨ s out (t). Here, r b = 8 a.u., which defines the spacial region where the electron is considered to be located at bound states [34]. Accordingly, the integral of the loss at the DI (SI) grid boundaries over time gives the total DI (SI) probability P Here, we focus on the main characteristics of PMD from HeH + versus HeT + . For clarity, we present the non-coherent results, i.e., c d (p 1 , . This TDSE treatment for HeT + is similar to that for HeH + with replacing H by T in relevant expressions.

III. RESULTS AND DISCUSSIONS
A. Asymmetric PMDs of DI In Fig. 1, we show calculated PMDs of DI for HeH + and HeT + at different laser parameters. These distributions indicate momentum correlation between these  III  II  IV  I  IV  III  II  two electrons in DI. Firstly, the distributions in Fig. 1 show a strong asymmetry with the amplitudes in the third quadrant remarkably larger than those in the first quadrant. Secondly, when increasing laser intensities or wavelengthes, the contributions of the first quadrant increase and this asymmetry becomes smaller. Thirdly, for the same laser parameters, this asymmetry is smaller for HeH + than for HeT + . By contrast, the distributions in the second and the fourth quadrants are similar for both isotopic cases. In the following, we concentrate on the origin of this asymmetry, associated with PMDs of DI in the first and the third quadrants.

B. Mechanisms of SI
To explore the potential mechanism, in Fig. 2, we plot the time-dependent ionization probabilities in one laser cycle, which is approximately evaluated with P (t) = 1 − n=15 n=1 | n|Ψ(t) | 2 [39] and is relating to SI. Here, |n is the nth electronic bound eigenstate of the field- at the BO approximation. Excluding more bound-state components from |Ψ(t) , results are similar to P (t). We divide the one-cycle time region into four parts denoted with I-IV. First, for all cases in Fig. 2, the ionization is strong in the first half laser cycle and is weak in the second half laser cycle. This phenomenon  Probabilities Figure 3: Probabilities of SI and DI of HeH + and HeT + vs laser intensity at λ = 500 nm. Ratios of DI to SI and ionization probabilities for 1D HeH 2+ and HeT 2+ are also plotted here.
has been termed as asymmetric ionization and is identified as arising from the effect of the permanent dipole [25]. Secondly, in the first half laser cycle, the ionization mainly occurs in the region II after the laser field arrives at its peak. The reason has been attributed to the Coulomb induced large ionization time delay. Due to this delay, many electrons which tunnel out of the laser-Coulomb formed barrier near the peak of the laser field in region I are ionized finally in region II [26]. By comparison, in the second half laser cycle, the contributions in region IV after the time of peak intensity also dominate the ionization for the same reason of Coulomb induced delay as in the first half cycle, but the contributions of region III are also non-negligible. In Ref. [39], it has been shown that the contributions in region III arise from effects of excited states. Specifically, some electrons are pumped into the excited states from the ground state around the peak of the laser field in the first half laser cycle and survive the falling part of the laser field of region II. Then the excited-state electrons with lower ionization potentials are ionized mostly in the arising part of the laser field of region III. Thirdly, as increasing the laser intensity and wavelength, the contributions of region III decrease due to the decrease of the excited state effect, as discussed in [20]. Fourthly, when the ionization yields of HeH + with lighter nuclei are larger than HeT + , a careful analysis tells that the ionization asymmetry in the first and the second half laser cycle is somewhat more remarkable for HeT + than for HeH + . These SI characteristics will be used to analyze the potential DI mechanisms. A simple evaluation on the drift momentum of the electron with the classical expression p = −A(t) [16] tells that electrons born in regions I and IV (II and III) have minus (plus) momenta. Here, A(t) is the vector potential of E(t). The time-dependent asymmetric ionization results in Fig. 2 Figure 4: Sketches of the laser field E(t) (shaded area) in one laser cycle (a) and motions of electrons corresponding to possible DI routes for HeH + or HeT + (b-c). In (a), the onecycle time region is divided into four parts (I-IV), labeled by different colors. The insets in (a) plot the laser-dressed asymmetric Coulomb potential, the laser-dressed electronic states |0 ′ and |1 ′ corresponding to the free-field electronic ground state |0 and the first excited state |1 of HeH + or HeT + , when the laser polarization is antiparallel (regions I and II) or parallel (regions III and IV) to the permanent dipole which is directing from the He nucleus to the H(T) nucleus. For the antiparallel case, the electronic ground state is dressed up and the first excited state is dressed down, making the ionization easier to occur. This situation reverses for the parallel case. As a result, the ionization is strong (weak) in the first (second) half laser cycle for the present cases. These analyses are also applicable for HeH 2+ or HeT 2+ . In (b) and (c), possible sequential (L1 and L3) and non-sequential (L2 and L4) DI routes associated with the first electron born in regions II (b) and IV (c) are plotted.
with plus momenta have large amplitudes. This point has been shown in Ref. [24] for vibrating HeH + with two-dimensional single-electron dynamics. Here, we also show PMDs of SI for the vibrating two-electron system of HeH + (HeT + ) as the insets in Fig. 2. The insets clearly show the asymmetric PMD of SI where distributions for plus momenta show larger amplitudes. This asymmetry in PMD of SI is more remarkable for HeT + than for HeH + . Combing the PMD results in Fig. 1 and Fig. 2, we arrive at the conclusion that these two electrons in DI prefer to release together along the H side when the first electron enjoys escaping along the He side in SI. From the analyses of the ionization in Fig. 2, one can also conclude that if sequential DI (SDI) occurs, PMDs of DI for HeH + or HeT + will show large amplitudes in the first quadrant,  similar to SI. We therefore anticipate that non-sequential double ionization (NSDI) dominates in present cases. In fact, extended TDSE simulations for 1D vibrating HeH 2+ (HeT 2+ ) show that for the present laser parameters, the ionization probability of HeH 2+ (HeT 2+ ) from its ground state is remarkably lower than the ratio of DI to SI for HeH + (HeT + ), as shown in Fig. 3, suggesting that NDSI dominates here [34]. Next, we discuss possible routes of NSDI in detail.

C. Mechanisms of DI
In Fig. 4(a), we plot the electric field E(t) in one laser cycle, with dividing the time region into four parts as in Fig. 2. The laser-dressed Coulomb potential and the laser-dressed two lowest electronic states of HeH + (HeT + ) corresponding to the first and the second half laser cycle are also plotted here as the insets. When the SI mainly occurs in regions II and IV, as discussed in Fig. 2, we focus on possible DI routes associated with SI events in these two regions.
First, in Fig. 4(b), we plot DI routes associated with the birth of the first electron in region II. In this case, as the first electron ionizes, the second electron can be ionized directly by the external field in this region with contributing to the first quadrant in PMDs of DI. It should be noted that due to the Coulomb induced large ionization time lag, these two ionized electrons in region II can find their origins in region I. We denote these SDI routes contributing to the first quadrant "L1". The first electron born in region II can also return to and recollide with the second electron in region III (short NDSI route) and region IV (long one). For the short route, the Coulomb effect will also induce the delay of DI time, just as it does in region I, resulting in the emission of these two electrons in region IV and contributing to DI in the third quadrant. For the long one, both electrons will also contribute to DI in the third quadrant. We denote these NSDI routes contributing to the third quadrant "L2". For lower laser intensities, probabilities for direct ionization of the second electron by the laser field are small, and the route L2 dominates in DI for the cases in Fig.  4(b). For DI routes associated with the birth of the first electron in region IV, this situation reverses, as plotted in Fig. 4(c). In this case, NSDI routes associated with the rescattering of the first electron contribute to the first quadrant and SDI routes related to sequential ionization of these two electrons contribute to the third quadrant. We denote these SDI and NSDI routes "L3" and "L4" respectively. Due to that the SI amplitudes are smaller in region IV than those in region II and the second electron is bounded more deeply in region IV than in region II (see the insets in Fig. 2(a)), for a laser cycle, the main contributions to DI come from route L2, with PMDs of DI showing large amplitudes in the third quadrant. When the laser intensity increases, the ionization yields of HeH 2+ associated with direct ionization by the laser field increase and the contributions of route L1 increase. As a result, PMDs of DI for the asymmetric system become more symmetric. For increasing the laser wavelength, the asymmetric system stretches to a large distance at which direct ionization is usually easier to occur, resulting in somewhat similar results as increasing the laser intensity.
For HeT + with heavier nuclei than HeH + , at the same laser parameters, the laser-induced stretching for HeT + is smaller than for HeH + , as shown in Fig. 5. Generally, direct ionization prefers larger R at which the ionization potential of the vibrating system is lower. Accordingly, direct-ionization yields for HeT + are also smaller than for HeH + . As a result, PMDs of DI for HeT + usually  show a stronger asymmetry than for HeH + . To validate these above discussions, in Fig. 6, we plot the ratio of amplitudes of PMDs of DI in the first quadrant to those in the third quadrant. One can observe from Fig. 6, as increasing the laser intensity or wavelength, for both isotope cases, this ratio increases. In particular, for the same laser parameters, this ratio of HeH + is usually larger than that of HeT + . These results are in agreement with our above analyses.

D. R-resolved DI and SI
To provide further insight into electron-nucleus coupled DI dynamics of polar molecules, in Fig. 7, we plot R-dependent probabilities γ s (R) of SI and γ d (R) of DI for HeH + and HeT + , averaged by the corresponding total probabilities P s and P d , respectively, at different laser parameters.
First, as increasing laser intensities and wavelengthes, for SI, the structure of these R-dependent distributions in Fig. 7 changes from a relatively sharp hump to the plane and broad one, but the center of the hump does not change basically. For DI, however, the distributions extend to somewhat larger distances. These different responses of DI and SI on laser parameters revealed here agree with the experimental results in [21]. Secondly, in some cases such as for HeH + , the location of the hump both for SI and DI is remarkably larger than the equilibrium separation R e = 2 a.u.. By comparison, for the symmetric case of model H 2 [34], the position of the hump is nearer to the equilibrium distance and the structure of the hump is not sensitive to laser parameters. The results suggest that due to the permanent-dipole induced rapid stretching which remarkably diminishes the ion- ization potential of the target, both DI and SI of the asymmetric system prefer to occur at larger R when the laser intensity is not very high. Thirdly, for the same laser parameters, the position of the hump for SI or DI of HeH + is larger than that for HeT + , suggesting that the stretching of HeH + with lighter nuclei is stronger than for HeT + , in agreement with previous discussions [23]. Fourthly, on the whole, the position of the hump for DI of HeH + or HeT + is near to the corresponding one for SI at lower laser intensities and shorter wavelengthes and is somewhat larger than that at higher and longer ones. The results suggest that at relatively high laser intensities and long wavelengthes, SI and DI events prefer to occur at times with a larger relative time delay. Such events are expected to be mainly associated with SDI processes. For these processes, after the first electron ionizes, the asymmetric system stretches for a while with lowering its ionization potential remarkably, then the second electron is ionized by the laser field, resulting in a SDI event. These analyses also support our above discussions that for HeH + and HeT + , SDI events increase as increasing laser intensities and wavelengthes.
In Fig. 8, we also plot R-dependent PMDs β d (R) of DI for HeH + and HeT + at some distances R at which the function γ d (R) has larger amplitudes (see Fig. 7). We have chosen the laser-parameter cases in Figs. 1(e) and 1(f) where the distributions in the first quadrant have relatively large amplitudes. For both cases of HeH + and HeT + , one can observe from Fig. 8, the main contributions to the third (first) quadrant come from smaller (larger) R, at which the ionization potential of the system is higher (lower) and we expect that NSDI (SDI) dominates.

IV. CONCLUSION
In conclusion, we have studied the ionization dynamics of vibrating HeH + with comparing it to HeT + . The photoelectron momentum distributions for both DI and SI show an asymmetric structure but with the contrary trend. As the Coulomb induced large ionization time delay plays an important role in the asymmetry in SI, the rescattering of the first electron along with this delay contributes to the asymmetry in DI. The nuclear motion mainly influences the events of SDI. The contributions of SDI decreases the asymmetry in DI momentum distributions, and this decreases is more remarkable for HeH + with lighter nuclei and more rapid nuclear motion than for HeT + . This asymmetry in DI is expected to appear for other oriented polar molecules with a large permanent dipole. In addition, the proposed DI mechanism of rescattering followed by Coulomb induced ionization time delay holds for all of atoms and molecules including symmetric and asymmetric ones. For these symmetric cases without a permanent dipole, this asymmetry discussed in the paper does not appear in PMDs of DI. However, effects relating to this mechanism are possible to resolve with using two-dimensional laser fields which have shown the capability in probing sub-cycle strong-field electron dynamics.