Dissociation and ionization of HeH + in sub-cycle-controlled intense two-color fields

Using quantum-mechanical, one-dimensional, non-Born–Oppenheimer simulations we study the control over the strong-field dynamics of the helium hydride molecular ion HeH + due to interaction driven by short and strong two-color laser pulses. We calculate yields of two competing fragmentation channels: electron removal and dissociation. We find that by changing the relative phase of the two colors, we can select the dominating channel. Nuclear motion is decisive for explaining ionization in this target. Ionization yields are vastly underestimated when nuclear motion is excluded and they are substantially reduced in the heavier isotopologue HeD + . Coupling of the two lowest electronic states is crucial even for the ground-state dissociation process.


Introduction
Understanding and controlling the interaction of molecules with short and strong laser pulses is a major aspect of coherent control and ultrafast dynamics [1][2][3][4]. A central goal is to be able to steer the outcome of chemical reactions. The simplest realization of this idea is the control over different fragmentation channels in laser-induced bond breaking. To this end, tailored femtosecond laser fields provide wide-ranging opportunities to manipulate molecules on their own time scale. Pulse-shaping techniques such as a spatial light modulator 5 Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. are traditionally used in femtosecond coherent control experiments [5,6]. They can generate complex pulse shapes with many control parameters. In high-intensity and attosecond physics, one usually resorts to conceptually simpler waveforms with less control parameters. In particular, by tailoring a laser field on the sub-cycle, i.e. attosecond time scale, we can impose a directionality on the laser field in the sense that adjacent half cycles are not equivalent. Two well known examples are laser pulses with controlled carrier-envelope phase (CEP) [7,8] and two-color fields, which are in the focus of this work.
The helium hydride molecular ion HeH + is one of the simplest diatomic molecules with an asymmetry in the electron configuration. As such it serves as a benchmark system. Although theoretical studies of HeH + exist [20][21][22][23][24][25], only few groups have carried out experiments with HeH + under the influence of light [26,27]. In a recent combined experimental and theoretical study at two different wavelengths, it was demonstrated that vibrational excitation significantly influences the ionization dynamics of HeH + [27]. To our knowledge, there are no previous investigations of HeH + or similar systems where dissociation and ionization are studied simultaneously.
Two-color laser irradiation is a simple and effective means of sub-cycle control. Already many years ago, it was shown that varying the relative phase between a laser pulse and a small portion of its second harmonic can control the breakup of atoms and molecules in strong laser fields [28][29][30]. Twocolor fields are also a frequently used tool for control and analysis of electron trajectories [31][32][33][34][35][36][37].
Here, we report on the results of a theoretical study of HeH + in two-color fields in order to identify laser parameters allowing for efficient control over two strong-field driven fragmentation channels: dissociation and ionization. We employ different models to study the roles of (i) nuclear motion and (ii) electronic excitation. We have chosen 1380 nm as the fundamental wavelength. This choice represents an intermediate value between small wavelengths where dissociation is negligible and large wavelengths where it becomes increasingly difficult to reach intensities high enough to observe ionization. Whereas previous studies of HeH + intentionally avoided laser parameters where dissociation and ionization compete, we explore this region in an attempt to shine light on the connection of the two. The chosen wavelength is easily experimentally accessible, e.g. using common optical parametric amplifiers pumped by Ti:Sa-based amplifiers [38]. This paper is organized as follows. In section 2, we describe the different models used in the simulations. In section 3, results are presented with and without intensity averaging. The pathways for different fragmentation channels are discussed in section 4. Section 5 gives a short conclusion of our work.

Model and methods
We solve the time-dependent Schrödinger equation (TDSE) for several models which are described in the following. They all have in common that the HeH + molecular ion is oriented along the laser polarization axis and that all movement (electronic and nuclear degrees of freedom) is confined to this axis as well. This is a reasonable approximation as discussed in [39].

Single-active-electron model with or without nuclear motion
As the most realistic model in this study, we use a singleactive-electron model system that includes both, one electronic and one nuclear degree of freedom. This model reproduces 3D Born-Oppenheimer potential-energy curves and has been used successfully to explain measured kinetic-energy release (KER) spectra [27].
The electron (coordinate x relative to the nuclear center of mass) and the nuclei (internuclear distance R) can move along the molecular axis. The time-dependent Schrödinger equation (TDSE) in dipole approximation and velocity gauge reads where p and P are the momentum operators of the electron and the internuclear separation, respectively, with κ = (m n + Z 1 + Z 2 )/(m n + 1), λ = (Z 1 m 2 − Z 2 m 1 )/m n . m 1 , m 2 and Z 1 , Z 2 are the nuclear masses and the effective charges on the proton and helium side respectively. The total nuclear mass is m n = m 1 + m 2 . The inactive electron is assumed to be fixed at the helium core, i.e. Z 1 = Z 2 = 1 and m 2 includes one electron mass. μ e = m e m n /(m n + 1), μ = m 1 m 2 /m n are the reduced masses. The electron-nuclear interaction is given by and V ion is the Born-Oppenheimer potential of HeH 2+ [40]. The soft-core parameters α 1 (R), α 2 (R) are tuned for each internuclear distance to reproduce the two lowest Born-Oppenheimer potential-energy curves of HeH + [41]. V ion and the two lowest Born-Oppenheimer potential-energy curves are shown in figure 1. The interaction with the laser field enters via the vector potential The wave function is initialized as a vibrational eigenstate with quantum number v of the electronic ground state. It is then numerically propagated on a grid using a split-operator scheme and Fourier transforms [42]. At the grid boundaries, absorbers are implemented using complex absorbing potentials [43]. Absorption at large |x| corresponds to electron removal, which we call ionization in the following (despite the fact that HeH + is already an ion). Absorption at large R corresponds to fragmentation without electron removal (i.e. the nuclei depart and the electron is bound to one of the nuclei), which we call dissociation in contrast to fragmentation by ionization. Usually, the grid is large enough to keep the amount of absorption at large |R| negligible. The ionized parts of the wave function from all time steps are coherently collected and propagated in the potential V ion and the laser field, neglecting the interaction term V ne .
Additional calculations are carried out for fixed nuclei. In this case, T R in (1) is removed, the wave function is prepared as the normalized electronic ground state at each internuclear distance, and propagation is carried out for all internuclear distances at the same time.
From the wave function at the end of the time evolution, the following observables are calculated: the total ionization yield (probability) is given by the norm of the collected absorption at large |x|; the total dissociation yield is calculated by projecting out bound states from the final wave function. The dissociation has many possible final states: He(n) + H + or He + + H(n) for any atomic eigenstate n. These are separated by projecting the final wave function onto selected electronic eigenstates for every internuclear distance. In the following, dissociation refers to the lowest dissociation channel, i.e. dissociation into a proton and the ground state of helium. This is the by far dominating dissociation channel. For both ionization and dissociation, the KER spectrum can be calculated by Fourier transforming the corresponding wave function along the R coordinate and integrating the modulus square over x.
The numerical grids have step sizes of 0.2 a.u. for the electron coordinate (4096 grid points) and 0.05 a.u. for the internuclear distance (2048 grid points). The time step size is 0.02 a.u. The absorbing potentials are given by monomials of degree 4 in the respective coordinate. For the electron, it covers an interval of 102.4 a.u. at each end of the grid; for the nuclear motion it covers an interval of 6.4 a.u. at the end of the grid. Eigenstates are calculated by repeatedly applying the shifted field-free short-time propagator exp(−iH 0 Δt) + 1 + i, where H 0 is the field-free Hamilton operator. After normalization, this converges to a state which is stationary in real time propagation and which corresponds to the smallest eigenvalue of H 0 . These states are numerically slightly different from eigenstates of H 0 (e.g. by diagonalization) or from imaginary time evolution.

Born-Oppenheimer calculations with or without electronic coupling
As the electronic ground state of HeH + is well separated from the first excited state (more than 20 eV, see figure 1), it is tempting to treat the dissociation process in the Born-Oppenheimer approximation. Neglecting the coupling between electronic states, the TDSE for the nuclear wave function ψ k (R) on the kth electronic potential curve V k reads Both, the (usually small) diagonal adiabatic correction to the Born-Oppenheimer approximation H ac,k and the coupling function λ k are calculated using the kth electronic eigenstate φ k (x; R) as a function of internuclear distance R, The coupling of the lowest electronic states is included in two-level calculations. Here, the TDSE reads i ∂ ∂t with H 1 , H 2 given by equation (5) and d 12 (R) = − φ 1 |x|φ 2 (x) . Equation (8) is solved by applying the split-operator scheme combined with the matrix exponential for the coupling term. Ionization is not included in this model. These calculations are carried out on one-dimensional numerical grids large enough such that absorbing boundaries are not necessary. The grid parameters are the same as for the calculations with electron motion. The dissociation yield is calculated from ψ 1 at the end of the time evolution (after the laser pulse) by projecting out the bound vibrational eigenstates of H 1 . The dissociation yield in the electronically excited state is directly provided by ψ 2 .

Electric field
The HeH + is exposed to two-color fields of the form E(t) = E ω (t) + E 2ω (t; φ) for 0 t T, where both E ω (t) and E 2ω (t) are pulses with a sin 2 envelope, Figure 3. Ionization (solid lines) and dissociation yields (broken lines) for HeH + (thick lines, left vertical axis) and HeD + (thin lines, right vertical axis) in a two-color field with 50 fs and 60 fs pulse duration for the ω and 2ω fields respectively and relative phase φ. The fundamental wavelength is 1380 nm. The peak intensity of the fundamental field is averaged from 1 × 10 14 Wcm −2 to 1 × 10 15 Wcm −2 assuming a flat target in a Gaussian beam profile (P avg. = P(I)dI/I) [44]. The intensity ratio of the two colors is 5:1. Simulations start in the v = 0 state of HeH + or HeD + respectively. Curves are calculated in steps of π/16 (33 data points per line).
φ is the carrier-envelope phase (CEP) of the 2ω field, T = N2π/ω is the total pulse duration. The full width at half maximum of the intensity envelope (FWHM) is 0.36 T. For φ = 0 the maxima of the electric fields coincide and the highest peak of the total field points from the proton towards the helium. The two-color phase φ plays the role of a delay between the two colors. The intensity ratio is fixed at 5:1. Throughout this work, values for intensities refer to the peak intensity of the fundamental field. An example of an electric field used in our calculations is shown in figure 2. Figure 3 shows the dependence of ionization and dissociation probabilities on the relative two-color phase φ after intensity averaging for HeH + and HeD + . We discuss the HeH + case (thick lines) first. Ionization is strongest when φ = 0 or φ = 2 π. In this case, the electric field maxima deform the molecular potential such that the electron tends to leave the molecule on the hydrogen side. This is the preferred ionization direction (in agreement with Dehghanian et al [24]). Dissociation, on the other hand, is favored at φ = π, i.e. when the preferred direction of the electric field is reversed and points from helium towards the hydrogen. By simply changing the two-color phase, the ratio of ionization versus dissociation signal can be changed from 6.3 to 0.74. Figure 4 shows similar results for a single intensity (i.e. without intensity averaging) and shorter pulses. The results for HeD + (thin lines in figures 3 and 4) are qualitatively similar, but they show some remarkable differences to the HeH + case: (i) overall, the yields are much lower,  Ionization (solid lines) and ground-state dissociation yields (broken lines) for HeH + (v = 0) in a two-color field with relative phase φ = 0 (purple asterisks) or φ = π (green circles) and variable pulse duration. The fundamental wavelength is 1380 nm and the peak intensity is 5 × 10 14 Wcm −2 . The intensity ratio of the two colors is 5:1.

Simulation results
(ii) the ratio of dissociation to ionization is weaker, and (iii) the suppression of ionization at φ = π is stronger than for HeH + .
For HeH + the dominant fragmentation process changes from ionization at φ = 0 to dissociation at φ = π. We investigate these two extreme cases in more detail. The dependence on the pulse duration is shown in figure 5. Data points from figure 4 can be found at the very left end of figure 5. The simulations show that for all considered pulse durations the yields of ionization at φ = 0 and dissociation at φ = π are of similar magnitude and scale approximately linearly with the pulse duration. When ionization is suppressed (φ = π), its yield grows approximately as T 3 . The dissociation yield for φ = 0 decreases for large pulse durations, probably because in long laser pulses, dissociating wave packets can still be ionized later in the pulse at large R. In this case, only dissociation that is initiated late in the pulse contributes to the dissociation yield.
In order to investigate the importance of nuclear motion for the ionization process, we carry out simulations with fixed nuclei. The internuclear distance-dependent ionization yields are plotted in figure 6(c) and show the region of enhanced ionization for R > 2. They are weighted with the  The results show that fixed-nuclei calculations dramatically underestimate the ionization yield, especially for the lighter isotopologue HeH + (where the fixed-nuclei yields have been increased by a factor of 10 for better visibility). In the case of HeD + , the agreement is better, which is not surprising as the nuclear motion is slower and thus less important during the laser pulse duration.
The fixed-nuclei curves in figures 6(a) and (b) have minima which are shifted away from φ = π. Although the enhanced ionization at increased R is much weaker for φ = π than for other values of φ (especially between R = 3 a.u. and R = 4 a.u., cf figure 6(c)), fixed-nuclei ionization is dominated by the region where the vibrational state is actually localized. In this region, the ionization probability is not minimal for φ = π, but for values near φ = 5π/8 and φ = 11π/8.

Discussion
The phase dependence of ionization dynamics can be understood considering the molecular orientation and nuclear motion. In general, asymmetric diatomic molecules have a preferred field direction for ionization. In the HeH + ground state, the electron is located at the helium site. The static ionization rate is increased when the field is directed such that the electron exits on the hydrogen side, compared to the situation where the electron exits on the helium side of the molecule. This is the main reason why ionization is much larger for φ = 0 than for φ = π in our asymmetric laser pulse.
Additionally, the internuclear distance plays a crucial role as resonances between the bound states can occur, leading to enhanced ionization at certain internuclear distances [24,45]. The KER spectrum for the fragments after ionization is displayed in figure 7(b). For every φ it is divided by the corresponding ionization yield (solid curve in panel (a)), i.e. the region around φ = π is amplified, revealing additional structures that are hardly visible for other delays. After the removal of the active electron, the nuclear motion of the remaining system evolves in the potential V ion (R) of the unstable HeH 2+ -basically a 1/R curve except at small R-gaining kinetic energy 1/R ion , where R ion is the internuclear distance at the time of ionization. Thus, low KER for the ionization is a sign of nuclear motion before ionization and subsequent ionization at increased internuclear distance [27].
The vibrational ground state has an average internuclear distance of R ≈ 1.5 a.u., i.e. direct ionization from the ground state should result in KERs around 0.66 a.u. Instead, we observe much lower KERs, indicating ionization from the region of enhanced ionization after vibrational excitation or even during the dissociation process.
This analysis also explains the discrepancy between the calculations with and without nuclear motion in figure 6. Nuclear motion is crucial for the ionization process. The strong connection between ionization and dissociation is supported by their similar behaviour as a function of the pulse duration, see figure 5. In both cases, the main part of the yield depends on how much time is available for the initial stretching of the molecule.
The KER spectrum for the dissociation is shown in figure 7(c). Clear multi-photon peaks can be seen, in agreement with calculations by Ursrey et al [22,23]. To identify the relevant pathways of the dissociation process, we compare results from 2D calculations (i.e. including electron dynamics) with 1D calculations (nuclear motion on the electronic ground state only or with coupling to the first excited state). These results are displayed in figure 8. The 1D model without coupling ('one level') can explain most of the phase-independent dissociation signal, but not the large modulation. Already the coupling to the first excited state ('two levels') is enough to increase the dissociation yield by more than a factor of 2 and to give a significant modulation. As ionization is stronger for φ = 0 than for φ = π, one can expect that including ionization in the two-level calculation would give quite good agreement with the 2D dissociation results.
A remarkable result is that the population of the excited state in the two-level calculation is much lower than the difference between one-level and two-level dissociation yields. This indicates that already the deformation of the groundstate potential curve due to the coupling is enough to alter the dissociation yield, or transient population in the excited state may play a role.
The instantaneous Born-Oppenheimer potential-energy curves at maximum field strength are shown in figure 9 for φ = 0 and φ = π. Even without coupling of the two states (dashed curve) the vibrational wave packet in the electronic ground state is temporarily completely unbound when the field points from He to H. This is the case in figure 9(b). The R-dependent coupling additionally lowers the barrier so that dissociation is increased compared to the case without coupling.
From the deformed potential curves for φ = 0 (figure 9(a)) we can infer the mechanism of the excitation pathway. In this case, there is a potential-curve crossing when the coupling is excluded. This arises because the permanent dipoles of the two states have opposite signs at large R where the ground state corresponds to He + H + while the excited state corresponds to He + + H. When the coupling is included, the curve Figure 9. Field-dressed instantaneous Born-Oppenheimer potential energy curves for the two-level system with (solid curves) and without (dashed curves) coupling between the two lowest states, calculated by diagonalizing the Hamiltonian operator in equation (8). The field strength is 0.14 a.u., corresponding to the maximum field in a fundamental intensity of 5 × 10 14 Wcm −2 . Field direction corresponding to (a) φ = 0, (b) φ = π. The black arrow in (a) indicates the pathway that leads to population of the excited state.
crossing turns into an avoided crossing with a barrier in the lower curve. A fraction of the already dissociating nuclear wave packet on the ground-state curve will pass the barrier in the avoided crossing region as indicated by the black arrow. This leads to population of the excited state after the end of the laser pulse.

Conclusion
Using TDSE simulations for HeH + and HeD + including both nuclear and electronic degrees of freedom, we have shown that the yields of fragmentation channels can be substantially controlled simply by changing the relative phase or delay of a two-color field on a sub-cycle time scale. Key ingredient is the asymmetric ionization rate for oriented asymmetric molecules and the coupling of electronic states, even for ground-state dissociation. Intensity-averaged calculations indicate that these effects should be observable in an experiment. As a striking result, we have found that there is a difference of π in the two-color delays that maximize ionization and dissociation.
The great importance of nuclear motion for the ionization yields has been shown by comparing predictions from fixednuclei calculations with those from simulations that include both degrees of freedom. The distribution of internuclear distances in the initial vibrational state is not enough to explain the ionization yields. When the nuclei are allowed to move during the laser pulse, they can reach higher distances and the ionization channel is dramatically enhanced.
HeH + is an extreme example of a diatomic polar molecule. The effects we have described can be expected to play a role-at least to some extent-in more complex polar molecules as well, possibly at other wavelengths. Therefore, applying tailored two-color laser fields to oriented molecules may be used to control fragmentation channels. The importance of nuclear motion will clearly depend on the reduced mass as well as the laser pulse duration. Therefore, for heavier molecules and short pulses, the nuclear motion will be less important for the prediction of the ionization yield.