Molecular vibrational trapping revisited: a case study with D2+

The present theoretical study is concerned with the vibrational trapping or bond hardening, which is a well-known phenomenon predicted by a dressed state representation of small molecules like and in an intense laser field. This phenomenon is associated with a condition where the energy of the light induced, vibrational level coincides with one of the vibrational levels on the field-free potential curve, which at the same time maximizes the wave function overlap between these two levels. One-dimensional numerical simulations were performed to investigate this phenomenon in a more quantitative way than has been done previously by calculating the photodissociation probability of for a wide range of photon energy. The obtained results undoubtedly show that the nodal structure of the field-free vibrational wave functions plays a decisive role in the vibrational trapping, in addition to the current understanding of this phenomenon.

The present theoretical study is concerned with the vibrational trapping or bond hardening, which is a well-known phenomenon predicted by a dressed state representation of small molecules like H 2 + and D 2 + in an intense laser field. This phenomenon is associated with a condition where the energy of the light induced, vibrational level coincides with one of the vibrational levels on the field-free potential curve, which at the same time maximizes the wave function overlap between these two levels. Onedimensional numerical simulations were performed to investigate this phenomenon in a more quantitative way than has been done previously by calculating the photodissociation probability of D 2 + for a wide range of photon energy. The obtained results undoubtedly show that the nodal structure of the field-free vibrational wave functions plays a decisive role in the vibrational trapping, in addition to the current understanding of this phenomenon.
Understanding the behavior of atoms and molecules in a strong laser field is an intensively studied area and there are numerous experimental and theoretical investigations which have uncovered and explored several new phenomena of light-matter interactions including high harmonic generation, above threshold ionization or dissociation [1][2][3][4] . The bond softening and the bond hardening effects [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] are similar phenomena which are often visualized during the photodissociation or photofragmentation processes of diatomic or polyatomic molecules. The mechanism of bond softening is well known. It was first discovered experimentally in the dissociation spectra of the + H 2 and + D 2 ions by Bucksbaum and coworkers [11][12][13] and can easily be understood by the illustrative Floquet's picture 25,26 or dressed state representation, which is often used to explain various strong field physics phenomena. The Floquet Hamiltonian for the net absorption of one photon can be represented by a 2 × 2 matrix, which explicitly includes the light-matter interaction. The change of the nuclear dynamics, due to the laser field, can be visualized as arising from "light-induced" or "adiabatic" molecular potentials (LIP) 27 . The molecular potential deforms due to strong radiative coupling and results in a strongly enhanced dissociation rate. The adiabatic curves show that as the laser intensity increases, the dissociation barrier moves to lower vibrational levels, leading to a noticeable growth of the dissociation probability. Bond hardening, which is often called molecular stabilization or vibrational trapping is the opposite effect with the same origins. Despite numerous experimental and theoretical works that have been performed in order to gain a better qualitative understanding of the mechanism behind this phenomenon [5][6][7][8][9][10][11][12][13][14][15][28][29][30][31] , there still remained some unresolved issues.
Recent efforts have been invested in studying the nature of the light-induced conical intersections (LICIs) 30,32 which was first discussed for diatomics by applying Floquet representation. In a diatomic molecule, which has only one nuclear vibrational coordinate, it is not possible for two electronic states of the same symmetry to form a conical intersection (CI). However, this former statement is only true in field-free space. Conical intersections can be formed, even in diatomic molecules, if an additional degree of freedom e.g. rotation, is associated with the system exposed to strong laser fields. In this situation, the interaction of the transition dipole moment of the molecule with the external electric field leads to an effective torque toward the polarization direction of the light. By considering the rotational coordinate in the description as a dynamical variable, the change of nuclear dynamics due to the external light can be considered as arising from a LICI. The positions of these LICIs are determined by the laser frequency while the laser intensity controls the strength of the nonadiabatic coupling. Several theoretical and experimental works have demonstrated that the LICIs have strong impact on the spectroscopic and dynamical properties of molecules 31,33-48 . There are many topical studies investigating the dissociation process of the + D 2 molecule within the LICI framework [37][38][39][40][41] . In these works it was found that using relatively small intensities (1 × 10 11 to 1 × 10 12 W/cm 2 ) and propagating the dynamics of the + D 2 results in very small (practically zero) rates for the dissociation probabilities of the photofragments 39 for certain vibrational eigenstates and photon energies. The dissociation yield still remains relatively small even at the highest studied intensity (1 × 10 14 W/cm 2 ). The obtained results are almost independent of the origins of the initial nuclear wave packet, which can start either from a vibrational eigenstate or a superposition of eigenstates, the Franck-Condon (FC) distribution, of the + D 2 . Intensity dependence at fixed frequency regarding the vibrational trapping phenomenon in the + H 2 and + D 2 ions was intensively theoretically studied by Bandrauk and coworkers [5][6][7][8]27 . They highlighted that laser-induced avoided crossings provide a very important source of molecular stabilization and hence suppress dissociation at high intensities. They also found that when the diabatic (field-free) and the adiabatic (light-induced) vibrational levels coincide, the resulting resonances can minimize the photodissociation probabilities. Increasing the field intensities, trapping of the initial state into stable adiabatic states occurs due to the decreasing nonadiabatic couplings between the light-induced states. Their work resulted in a qualitatively correct picture of the vibrational trapping phenomenon mechanism.
In this article we investigate why, in certain situations, there is virtually no dissociation rate of the photofragments arising from the dissociation process of the + D 2 molecule. This paper expands previous investigations by trying to find a more accurate and quantitative connection between physical quantities including the coincidence of the diabatic and adiabatic vibrational levels, maximum overlaps between the diabatic and adiabatic wave functions and minima of the photodissociation probabilities. This is realized by calculating the difference between the diabatic and adiabatic vibrational energies, the overlap between the diabatic and adiabatic wave functions and the dissociation yield of the molecular photofragments over a wide range of photon energies (from ω 1 = 24.799 eV, λ 1 = 50 nm to ω 2 = 3.100 eV, λ 2 = 400 nm).
As was discussed in 39,40 , the bond hardening effect is highly dependent on the behavior of the system in low intensity field even if the actual applied intensity does not fall into this region. This is related to the fact that-not considering few cycle pulses-larger intensity pulses are built up gradually. Consequently, at the beginning of the pulse the role of the different adiabatic levels is determined by the low intensity behavior of the system. For this reason we use 1 × 10 11 W/cm 2 intensity in our simulations or even the zero intensity limit to determine the different features of the vibrational eigenstates of the adiabatic upper state. In the simulations we included only the two lowest lying electronic states of the + D 2 molecule. This decreases the validity of the results for the higher energy region-especially for the λ ≲ 100 nm-where the Rydberg states of the ion and even the doubly-ionized repulsive Coulomb states can play an essential role in the real processes. Fortunately, this inaccuracy in the description of the physical system does not invalidate our findings related to the bond hardening effect, which is widely studied within the two state approximation. One can consider that some of our results are valid only for a model system inspired by the + D 2 molecular ion.
Theory and Methods Figure 1 shows the potential energy curves for + D 2 in dressed state representation. The electronic ground (V 1 (R) = 1sσ g ) and the first excited (V 2 (R) = 2pσ u ) eigenstates are considered as diabatic potentials and together, with the kinetic energy, form the field free Hamiltonian where R is the vibrational coordinate and μ is the reduced mass. The ion is excited by a resonant laser pulse from the V 1 (R) ground state to the repulsive V 2 (R) state. An electronic transition occurs due to the nonvanishing transition dipole moment and the corresponding 2 × 2 field dressed Hamiltonian matrix reads The off-diagonal elements of eq. (2) represent the radiative couplings, where  0 is the laser field amplitude, ω is the laser frequency which couples the two electronic states, d(R) is the transition dipole and θ is the angle between the polarization direction of the light and the direction of the molecular axes. The potential energies of V 1 (R) and V 2 (R) and the transition dipole moment were sourced from refs 49,50. The exact, time-dependent (TD) form of eq. (2), which is used extensively in the time-dependent calculations is also given here where f(t) is the envelope function of the laser pulse. A convenient and approximate form for the adiabatic potential is assumed. The energy levels and the eigenfunctions of the upper adiabatic potential can easily be calculated for weak fields from a model potential (Fig. 2): where Θ is the Heaviside step function, λ is the wavelength and R CR is the crossing point between the ground V 1 (R) and the field dressed (V 2 (R) − ħω) excited states. This model potential in eq. (4) can be converted into the upper adiabatic potential at zero-field limit. It is now possible to calculate the vibrational energy levels (E ν ) and wave functions (ψ ν ) belonging to the different vibrational eigenstates of the diabatic potential in eq. (1) and the eigenvalues ( ν′  ) and eigenfunctions (ϕ ν′ ) of the upper adiabatic potential in eq. (4). A more detailed analysis requires new quantities derived from different combinations of the previously determined basic quantities. These are the difference of the adiabatic and diabatic eigenvalues , the adiabatic and diabatic eigenfunctions overlap, , and the total dissociation probability: diss 0 The ΔE ν,ν′ (λ) and S ν,ν′ (λ) were determined with a spacing of Δλ = 1.0 in the interval of (50 nm-400 nm), −iW is the complex absorbing potential (CAP) applied at the last 5 a.u. of the grid related to the vibrational degree of freedom (W = 0.00005 · (r 70) 3 , if r > 70 a.u. on the 1sσ g surface and W = 0.00236 · (r 75) 3 , if r > 75 a.u. on the 2pσ u surface.
The multi-configuration time-dependent Hartree (MCTDH) method 51-55 is used to solve the time-independent and time-dependent Schrödinger equations in conjunction with the eigenvalue and eigenfunction problems in the diabatic/adiabatic frameworks, (eqs (1 and 4)) and the dissociation dynamics (eq. (3)). Fast Fourier transformation-discrete variable representation (FFT-DVR) 56 is used to characterize the vibrational degree of freedom, with N R basis elements for the internuclear separation distributed within the range from 0.1 a.u. to 10.05 a.u. or 80 a.u. in the eigenstates or the dissociation yield calculations, respectively. These primitive basis sets (χ) are used to represent the wave function In the actual calculations, N R = 256 is used for the eigenstates and N R = 2048 for the dynamical calculations, respectively.

Results and Discussion
Our present studies intensively probe the effect of molecular stabilization using accurate numerical calculations, performed in a wide range of energies. Detailed analysis is performed using a one-dimensional (1D) nuclear dynamical simulation and the initial nuclear wave packet is assumed to be in one of the vibrational eigenstates (ν = 0, 1, 2, 3, 4, 5, 6,7,8,9) or Franck-Condon distribution of the vibrational states of the ion. The initial orientation of the molecules is assumed to be parallel to the external field in all cases. The dissociation rates for isotropic initial distribution can be approximated with large accuracy by dividing the results by 3 as this study is limited to low intensities. The energy interval (ω 1 = 24.799 eV, λ 1 = 50 nm to ω 2 = 3.100 eV, λ 2 = 400 nm) is chosen, so that it contains all the values belonging to the positions of the near zero dissociation probability of the studied vibrational levels. A linearly polarized Gaussian laser pulse, with intensity, I 0 = 10 11 W/cm 2 and a pulse duration in full width at half-maximum (FWHM) t pulse = 30 fs is used throughout the calculations.
Three physical properties, the difference of the diabatic and adiabatic eigenvalues (ΔE(λ)), the diabatic and adiabatic eigenfunctions' overlap (S(λ)) and the dissociation probability (P diss ) play an important role in the forthcoming analysis.
Dissociation probability, wave function overlap and energy difference. The dissociation probability, eq. (7), is calculated in the chosen energy interval by maintaining a fixed light intensity. Wavelengths which result in a virtually zero dissociation probability from a particular vibrational energy level are sought, like in previous investigations 39 . Calculations have been performed for all chosen vibrational eigenstate, ν = 0, 1, …, 9, and the previously determined energy interval ensured that for each eigenstates all the wavelengths which belong to a minimum value of the dissociation probability were assigned. For example, nine different values for the wavelength were found for the ν = 9 vibrational level. Figure 3 shows the results for the ν = 4 vibrational level for the case when the initial nuclear wave packet starts from (i) a vibrational eigenstate or from (ii) a superposition of eigenstates, the Franck-Condon, (FC) distribution, of the + D 2 . The dissociation rate has four minima at the λ D (4, 3) = 88.8 nm, λ D (4, 2) = 112 nm, λ D (4,1) = 139.5 nm and λ D (4, 0) = 177.1 nm wavelengths. Simulations in which the wave packets start from a vibrational eigenstate use pulses with a 30 fs duration centered around 0 fs and the total dissociation probabilities are calculated. In the superposition situation, pulses with a 30 fs duration are applied. These pulses are centered around 34 fs during which time approximately one and a half vibration cycle takes place on the diabatic lower surface V 1 38 . The yield of total dissociation probability does not provide any information about the dissociation rate resulting from a particular vibrational state because all of the vibrational eigenstates are included in the initial wave packet. Thus, attention turned to the kinetic energy (KER) release spectra and the dissociation rate at the E ν + ħω energy was studied.
Results for the dissociation yield obtained by assuming FC distribution of the nuclear wave packet are surprisingly similar to those obtained from initiating the dynamics from one of the vibrational eigenstates of the Hamiltonian. This similarity is related to the fact that, on one hand the applied laser pulse is long enough to resolve properly the vibrational energy eigenstates and on the other hand its intensity is in the linear respond regime. These two features allow us to conclude that the structure of the initial wave packet does not significantly affect the underlying physical mechanism in the bond hardening effect. (Not displayed results for other vibrational states clearly confirm this statement, as well). For the sake of simpler analysis initial wave packets assuming vibrational eigenstate form are used in the further calculations.  The next task was to obtain the diabatic and adiabatic eigenfunctions' overlap which was calculated using eq (6). For each diabatic vibrational level (ν) ν′ = 0, 1, … ν − 1, adiabatic eigenstates and correspondingly ν different overlap functions exist. Figure 4 shows these functions for the ν = 4 vibrational level. The 4 different panels belong to the corresponding values of ν′ from 0 to ν − 1, and the obtained curves are scaled on the left side of the panels.
The difference of the adiabatic and diabatic eigenvalues eq (5), is determined and ν energy difference curves are obtained for each value of ν. Figure 4 shows these curves which are scaled on the right side of the panels. It can be expected, based on previous theoretical studies 8,9,27 that, for a given pair of diabatic and adiabatic eigenstates with the feature of ν > ν′, both the maxima of the overlap between the diabatic and adiabatic vibrational wave functions at λ = λ O (ν, ν′) and the matching of the diabatic and adiabatic eigenvalues at λ = λ E (ν, ν′) are closely related to the minimal dissociation yield or bond hardening situation, which takes place at special photon energies λ = λ D (ν, ν′). These special values of the wavelength are presented in the first four columns of Table 1 for all the studied vibration levels (ν = 0, 1, …, 9). If the values of these three wavelengths (λ D , λ O and λ E ) are close to each other for a certain pair of vibrational frequency (ν, ν′), the values of the (λ O )/(λ D ) and (λ E )/(λ D ) ratios are also close to 1. However, Fig. 5 shows that this is not the case as the (λ O )/(λ D ) and (λ E )/(λ D ) wavelength rates systematically deviate from 1 within a range of 5-15% for all cases. The significantly large differences between the wavelengths suggest that the physics behind the trapping effect may not be described adequately by studying the overlap of the wave functions or energy differences and it is expected, that some other factors may play a role.
The role of the nodal structure. This section starts with a simple analysis of the adiabatic and diabatic wave functions overlap. It is obvious that the node positions of adiabatic wave function should be close to the node positions of the diabatic function in order to obtain a significantly large overlap. The adiabatic wave function fades to zero at such internuclear distances where the adiabatic potential is above the adiabatic energy level. Applying a photon energy which corresponds to the dissociation minimum results in the energies of the diabatic and adiabatic states being similar. In such a case, the disappearance of the adiabatic wave function can be expected to approximately take place where the adiabatic potential is equal to the energy of the diabatic eigenstate. Thus, it can be assumed that for a large overlap the matching of the adiabatic potential to the diabatic energy level should take place around one of the nodes of the diabatic eigenfunction. The wavelength λ ν ν′ ν ( , ) E is where the matching of the adiabatic potential to the energy of the diabatic eigenstate takes place at the (ν − ν′) th node (R n (ν, ν′)) of the diabatic wave function, i.e., after which there are additional ν′ nodes (Fig. 6): There is another possible way to demonstrate the importance of the nodes of the nuclear wave packet in the process of molecular stabilization. In the light-induced picture of dissociation, the nonadiabatic coupling between the two surfaces will be the largest around the crossing of the ground and the dressed excited state diabatic potentials. Therefore, a reasonable reduction in the dissociation yield can be expected if the initial diabatic wave packet has a node at, or close to this crossing. To check the validity of this approach, the wavelengths for which the crossing of the potentials takes place at one of the nodes of the diabatic wave function was determined. These values are denoted as λ × (ν, ν′) and the corresponding photon energies can be calculated as  Table 1 shows both of these newly introduced quantities and the positions of their dependant nodes. The accuracy of the above constructed approach was checked by calculating the ratios with λ D and the obtained (λ × )/(λ D ) and λ λ ν / E D rates are shown in Fig. 5a. For lower lying vibrational states, both these two new formulas, based on the position of the nodes, provide a better forecast of the wavelength of the minimal dissociation yield than the preliminary defined λ E or λ O values. λ × underestimates λ D approximately twice as much as λ ν E overestimates. This observation can be used to construct some additional semiempirical formulas which might balance these two opposite deviations in the λ × and λ ν E wavelengths (see in Fig. 5a). In these formulas the λ ν E wavelength always takes part with twice the weight in. Obtained combinations are: the weighted arithmetic, geometric and harmonic means (arithmetic mean of the photon energies) of the former two node based approximations for the position of the minimal dissociation rate: a E a E Figure 5. Ratios of the different wavelengths with respect to the minimum value of the dissociation probability λ D . Panel (a): λ O is the wavelength belonging to the maximum value of the wave functions' overlap, λ E is the wavelength where the adiabatic and diabatic eigenvalues are the same, λ ν E is the wavelength according to eq. (9) and λ × is the wavelength according to eq. (10); Panel (b): λ a , λ g and λ h are the weighted arithmetic, geometric and harmonic means of λ ν E and λ × as defined in eqs (11)(12)(13).
Scientific RepoRts | 6:31871 | DOI: 10.1038/srep31871 The ratios of the above constructed mixed quantities to the dissociation minimum (λ D ) are shown in Fig. 5b. All the three formulas lead to very similar results predicting the dissociation minima with less than a half percent error, except for a few cases. The harmonic mean, defined in eq. (13), is the most balanced combination as its values remain in the narrowest interval around λ D . The fact that these wavelengths, especially their special combination provide an almost perfect molecular stability, implies that the nodes of the nuclear vibrational wave packet play an essential role in the bond hardening effect.
The above numerical results support the theoretical expectations: initiating the nuclear dissociation dynamics from a given vibrational level, one can always find the appropriate photon energy that minimizes the dissociation yield. The latter holds for the total dissociation probability when the initial wave packet starts from one of the vibrational eigenstates and for the kinetic energy release spectra (KER) at a certain (E ν + ħω) energy when the initial wave packet starts from a FC distribution of ionic vibrational eigenstates. At the heart of this work are the nodes of the nuclear vibrational wave packet which can quantitatively be associated with the bond hardening effect.
Our showcase example throughout this study was the + D 2 molecule. One may assume that the present findings are independent from + D 2 and are valid for any other molecular systems, as well.

Conclusions
In this paper, a series of numerical calculations have been performed, over a wide range of photon energies in order to provide accurate results for the vibrational trapping phenomenon. It was demonstrated that the nodal structure of the nuclear vibrational wave packets have an important role in the molecular stabilization process. We have found an accurate quantitative connection between the nodes of the nuclear vibrational wave packet and the minimum of the dissociation rate. In contrast to the literature, the importance of the nodes of the vibrational wave packet is strongly stressed because an almost perfect molecular stabilization can be obtained using our recently proposed formula (eq. (13)) of the node based wavelength. However, we have yet to clearly identify the background mechanism behind the harmonic mean combination of the λ ν E and λ × wavelengths, which can provide the best description of the phenomena. Our studies will continue and will be extended in this direction and we expect to form a better understanding about the complex physical phenomena beyond the vibrational trapping effect. Nevertheless, we hope that the present findings will stimulate experimental works in the near future.  O (4, 1)) with the diabatic ν = 4 eigenstate. The zero line for the wave functions is placed to the energy level of the ν = 4 vibrational state. The V 1 diabatic and V 2 potential energies with five field-dressed excited states are also shown. The solid black line denotes the shift by a laser field with the wavelength λ D (4, 1) (ω D (4, 1) frequency), for which the dissociation of the + D 2 is minimal. Lines marked with empty squares and circles are shifted by the energies of ħω O and ħω E , respectively. The line marked with solid diamonds corresponds to a shift with ω E 4 , so that the energy of the dressed state at the third node (R n (4, 1)) of the ν 4 diabatic eigenfunction is equal to the E 4 vibrational eigen energy. The ħω × is the photon energy (λ × (4, 1) wavelength) by which shifting the V 2 state (marked with triangles) a crossing with the ground state potential is formed at the position of the R n (4, 1) node. The inset in the bottom right corner displays a zoomed version of the central area.