Interaction of the hydrogen molecule with the environment: stability of the system

We study the stability of the hydrogen molecule interacting with the environment according to the balanced gain and loss energy scheme. We determined the properties of the molecule taking into account all electronic interactions, where the parameters of the Hamiltonian have been computed by using the variational method. The interaction of the hydrogen molecule with the environment was modeled parametrically ($\gamma$) with the help of the non-hermitian operator. We have shown that the hydrogen molecule is dynamically unstable. The dissociation time ($T_{D}$) decreases, if the $\gamma$ parameter increases (for $\gamma\rightarrow 0$, we get $T_{D}\rightarrow +\infty$). At the dynamic instability of the hydrogen molecule overlaps its static instability as the coupling constant $\gamma$ increases. We observed the decrease in the dissociation energy and the existence of the metastable state of the molecule ($\gamma_{MS}=0.659374$~Ry). The hydrogen molecule is statically unstable for $\gamma>\gamma_{D}=1.024638$~Ry. One can also observed the $\mathcal{PT}$ symmetry breaking effect for the electronic Hamiltonian ($\gamma_{\mathcal {PT}}=0.520873$~Ry). However, it does not affect the properties of the hydrogen molecule, such as: the electronic Hamiltonian parameters, the phonon and rotational energy, and the values of the electron-phonon coupling constants.

(2) is the abbreviation for plus the Hermitian conjugate and it means that an additional term being the Hermitian conjugate of the preceding term should be added. The Hamiltonian parameters are defined by the following integrals:  The meaning of above quantities is as follows: ε represents the energy of the molecular orbital, t is the electronic hopping integral, U denotes the on-site Coulomb repulsion, K is the energy of the inter-site Coulomb repulsion, J stands for the integral of the exchange, and V is called the correlated hopping. The integrals were calculated numerically, which is the complicated procedure that requires the use of the large computer resources. We notice that the contribution of the individual integrals to the energy eigenvalues is very diverse (see Table 1), nevertheless omitting any interaction would lead to the non-physical shortening of the distance between protons. We chose the Wannier's functions in the form of: where the coefficients ensuring normalization are expressed in the formulas: The atomic overlap (S) has the form: 3 1 2 , where 1 s Slater-type orbital can be written as: , α is the inverse size of the orbital. It should be noted that the second quantization method is completely equivalent to the Schrödinger analysis [38][39][40][41][42] .
The effective interaction of hydrogen molecule with the environment will be taken into account by supplementing the Hubbard Hamiltonian Ĥ H e with the balanced gain and loss operator 3,43 : 1 2

H H
The interpretation is that the operators î n 1 γ and γ −i n 2 enable the effective description of the inward and outward fluxes of the probability amplitude (the interaction with the environment) 3,44 . The addition or removal of an electron from the system would look more realistic, but this approach requires solving the master equations [45][46][47] , which is a quite complicated numerical task. Notice that the operator = + represents the non-Hermitian Hamiltonian, nevertheless it remains invariant due to the P PT T symmetry -at least to the characteristic γ P PT T value for which the symmetry is broken.

Results
Static stability of the system: the electron, the phonon and the electron-phonon properties.
In the Fig. 1, we plotted the dependence of the eigenvalues E j on γ. The analytical formulas for E j have been collected in the Appendix. Analyzing the obtained results, we found that for P PT T γ = . 0 520873 Ry there occurs the breaking of P PT T symmetry of the electronic Hamiltonian. This fact is manifested by the appearance of the complex values of E 5 and E 6 . Physically, this means that the P PT T symmetry breaking reduces the number of the available electronic states from six to four. Nevertheless, the considered effect has no physical significance due to the fact that the states E 5 and E 6 have the highest energy values. They cannot be thermally occupied -the k B T energy is of the order of 25 meV, while the difference between E 6 and E 4 is around 25.5 eV (E 4 is the ground state energy of the electronic subsystem). When discussing the results, it should be clearly emphasized that E 4 is always the real number.
Although the P PT T symmetry breaking is not manifested physically, the interaction of the hydrogen molecule with the environment can significantly change its physical state. This fact is connected with the dependence of the total energy on the γ parameter. In the Fig. 2, we presented the total energy values ( hydrogen molecule, and the influence of the γ parameter on the ground state energy E T (4) . One can see that an increase in the minimum energy value E R ( ) T (4) 0 is observed with the increase in γ, whereas the molecule is in the stable state. Above γ = .

MS
Ry the hydrogen molecule can exist only in the metastable state: Ry, where = .
) on the value of the γ parameter. The insert (b) shows the influence of the γ parameter on the equilibrium distance R 0 . Figure 3(a-c) trace the change of the distribution of electron charge for the stable case at γ = 0 (R 0 = 1.41968 a 0 ), for the metastable case (R MS 0 ( ) ), and at the dissociation point (R D 0 ( ) ). The density of electron charge was calculated according to the formula: The determination of the explicit function E T (R) for the given parameter γ allows to trace the influence of the environment on the vibrational energy. In the simplest approach (the harmonic approximation), the potential can be calculated as follows: . The energy of quantum oscillator has the form: The symbol n indexes the energy level: = . .. . The more advanced approach is based on the Morse potential: . The Morse energy is given by Table 2). The energy formula has the more complex form than for the harmonic case: Mo ω on the value of the γ parameter. There is a clear difference in the courses of the functions under consideration. It results from the applied method of approximation of the exact dependence of the total energy on the inter-proton distance (see Fig. 4(b)). It is worth noticing that the anharmonic approximation can be used only for the γ values smaller than γ MS ; the Morse curve incorrectly parameterizes the ground state energy function E R ( ) for higher values of γ. The rotational energy of the hydrogen molecule should be calculated from the expression: Having the explicit dependence of the ˆγ H H e parameters on R (see the appendix), we computed the electron-phonon coupling functions according to the formula: We plotted the obtained results in the Fig. 5. One can easily see that the absolute values of the considered functions at R 0 increase as the γ parameter increases. The couplings associated with the ε, t, U, and K parameters are of the greatest physical importance. The other two quantities g J and g V take very small values as compared to other electron-phonon coupling functions. Note the relatively high values of the g U and g K functions. The obtained result is caused by the fact that the electrons in the hydrogen molecule form the strongly correlated system. the dynamic instability of the hydrogen molecule. The basic observable of the electronic subsystem is the occupation number σ n j of the j th proton of the hydrogen molecule, where the symbol ... means the expectation value. In the Hermitian case (γ = 0) the dynamics of σ n j can be analyzed using the conventional Heisenberg equation:ˆˆĤ is the electron Hamiltonian in the mean-field approximation:     The new symbols have the following meanings:  , and ∆ = ↓ ↑ˆĉ c j j j . It should be emphasized that for ˆγ H H e the required operator calculations are not feasible due to their extensiveness. The mean-field approximation transforms the operator ˆγ e H H into the Hamiltonian, in which the energy of the molecular state and the hopping integral explicitly depend on the proton index j and the spin σ. In addition, the Hamiltonian have the part that models the reversal of the spin due to the exchange interaction J. It is also worth paying attention to the quantity of ∆ j , which has the formal structure of the Cooper pair annihilation operator in the real space. This analogy is not complete, because the Hamiltonian ˆγ e MF H H term containing ∆ j and ˆ † ∆ j does not correspond to BCS pairing operator [48][49][50] (the integral of the exchange J 0 has the positive value instead of negative -see Table 1).
After performing the required operator calculations, we get the set of sixteen first-order differential equations, which is explicitly written in the appendix. In the non-Hermitian case (γ ≠ 0), determining the time dependence of the electron observables is the more subtle issue 16,25,51 . First of all, one must define the operators: Tedious, but not difficult operator calculations lead to the complex system of the differential equations, which is presented in the appendix in the explicit form.
We plotted the time dependence of the observables ↑ n 1 , ↑ n 2 for the selected values of the γ parameter in the Fig. 6(a-e). As expected, the system in the Hermitian case is in the stable state, which manifests itself by the time invariance of the expectation value. In the non-Hermitian case, the weak interaction of the hydrogen molecule with the environment (γ = 0.1) causes oscillatory changes of the discussed quantities in time. However, this is not time-stable state of the system, because from the specific moment T D observables ↑ n 1 , ↑ n 2 accept complex values. From the physical point of view, the time T D should be interpreted as the moment in which the system dissociates. It is easy to show that as the γ parameter increases, the oscillations of the expectation values disappear and the value T D decreases very clearly (see Fig. 6(f)). The obtained results mean that any weak interaction of the hydrogen molecule with the environment modeled in the BGL scheme leads to the finite life time of the molecule.

Summary and Discussion of the Results
The obtained results show that the BGL type interaction of the hydrogen molecule with the environment leads to its dissociation. From the physical point of view, this means that the hydrogen molecule breaks down into two hydrogen atoms. Note that if the interaction of the hydrogen molecule with the environment would be modelled in the unbalanced gain and loss energy scheme, two other final states could be obtained: − H 2 or + H 2 . Our result is caused by the dynamic instability of the electronic subsystem. Note that the dynamic instability of the molecule is superimposed on the static instability for high values of γ parameter. We showed that the increase in the value of γ strongly reduces the dissociation energy of the molecule. Above γ = .
0 659374 MS Ry, the molecule is in the metastable state, decaying definitively for γ > .

D
Ry. An additional effect, that we observed for γ higher than γ = . 0 520873 P PT T Ry, is the P PT T symmetry breaking of the electronic Hamiltonian γ e H H . As a result, the two highest energies of the electron state assume complex values and the number of available electronic states of the molecule is reduced to four. This effect does not influence the stability of the considered system. Additionally, the P PT T symmetry breaking does not change either the values of the integrals of the electronic Hamiltonian, or the phonon or rotational properties of the hydrogen molecules, or the electron-phonon interaction constants. The dynamics of the electronic subsystem is also independent on the breaking of the P PT T symmetry of ˆγ H H e . It should be noted that all of the mentioned topics have more than just an academic value. Regarding the area of modern technology, particular attention should be paid to the nanoelectronic section linked to the molecular junctions [52][53][54][55][56][57][58][59][60][61] . Particularly interesting are the hydrogen molecule-bridged junctions of the type X/H 2 /X, where symbol X means metals like Pt 52,53,56-58 , Pd 55 , Au 58,59 , Cu 60 or Ni 61 . It is obvious that in nanojunctions there is no issue with the stability of a molecular bridge interacting with environment, because of the whole system being stabilized by electrodes. However, it does not mean that the issue of reducing the molecular levels of the bridge caused by the correspondingly strong interaction of the hydrogen molecule with the electrodes of the joint can be omitted (γ γ > ′ P PT T , whereas γ ′ P PT T means the value of γ parameter for which the PT symmetry breaking of the electronic bridge sub-system happens in the junction).
It should be emphasized that the formalism presented in the work enables the detailed analysis of the electronic structure of the hydrogen bridge in a nanojunction. For this purpose, the initial determination of physical parameters of the considered nanojunction, particularly of the equilibrium distance between the hydrogen atoms in the bridge ( ′ R 0 ), should be done according to the method based on the density-functional theory (DFT) 62 . The physical state of the bridge when there is no flux (γ = 0) corresponds to the minimum enthalpy value: , where the symbol F denotes the force exerted on the bridge by the junction anchors. The value of F, within the scheme presented in the work, should be selected so that the minimum enthalpy value is www.nature.com/scientificreports www.nature.com/scientificreports/ achieved for the ′ R 0 distance. Further calculations in order to characterise the electronic structure of the bridge for γ ≠ 0 should be performed according to the presented scheme, applying the generalised formula for the enthalpy: Noticeably, the dynamics of electronic observables of the molecular bridge interacting with electrodes should not be analyzed with the use of classical Heisenberg equation, but rather with a formalism of non-Hermitian quantum mechanics 16,25,51 .  www.nature.com/scientificreports www.nature.com/scientificreports/ where: . By using the operator (15) there was brought out the preliminary formulas for the eigenvalues, which has the form as follows:   4( 2 6 )), wherein: An attentive reader will notice that the energies ε, t, U, etc. are explicit functions of the inter-proton distance R and the parameter α. In the Fig. 7, we plotted the discussed values of the energies as the function of R and γ. Additionally, in Table 1 we give the equilibrium values of the ˆγ e H H parameters. The explicit dependence of the variational parameter α on the distance R is shown in the Fig. 8. the equilibrium values of phonon energy, rotational energy and the electron-phonon coupling function. In the Tables 2 and 3, we collected the equilibrium values of the phonon parameters for selected γ.
The Table 4 presents the equilibrium values of the electron-phonon coupling functions. The Table 5 presents the equilibrium distance R 0 , the equilibrium inverse size of the orbital α 0 , and the ground-state energy E T www.nature.com/scientificreports www.nature.com/scientificreports/ The set of the differential equations for electron observables (γ = 0). The system of differential equations has the form:ˆˆˆ †         The system of differential equations for electron observables (γ ≠ 0). The system of differential equations can be written in the form: