Interaction of the hydrogen molecule with the environment: stability of the system and the ${\mathscr{PT}}$ symmetry breaking

Abstract

The stability of the hydrogen molecule interacting with the environment according to the balanced gain and loss energy scheme was studied. We determined the properties of the molecule taking into account all electronic interactions, the parameters of the Hamiltonian being computed by the variational method. The interaction of the hydrogen molecule with the environment was modeled parametrically (γ) by means of the non-Hermitian, ${\mathscr{PT}}$-symmetric Hamiltonian. We showed that the hydrogen molecule is dynamically unstable. Its dissociation time (T_D) decreases if the γ parameter increases (for γ → 0 we got T_D → + ∞). The dynamic instability of the hydrogen molecule is superimposed on the decrease in its static stability as γ increases. Then we can observe the decrease in the dissociation energy value and the existence of the metastable state of the molecule as γ_MS reaches 0.659374 Ry. The hydrogen molecule is statically unstable when γ > γ_D = 1.024638 Ry. Moreover, we can also observe the ${\mathscr{PT}}$ symmetry breaking effect for the electronic Hamiltonian when ${\gamma }_{{\mathscr{PT}}}$ = 0.520873 Ry. This effect does not affect such properties of the hydrogen molecule as: the electronic Hamiltonian parameters, the phonon and the rotational energies, and the values of the electron-phonon coupling constants neither it disturbs the dynamics of the electronic subsystem. However, the number of available quantum states goes down to four.

Introduction

Research on the impact of the environment (external quantum system) on the state of the quantum system is an interesting, but very difficult issue^1,2. This is due to two reasons: (i) usually in the case of the realistic quantum system it is impossible to accurately determine its internal state due to the mathematical complexity of the problem, and (ii) the interaction between the quantum system and the environment can be so complicated that it is impossible to obtain unambiguous results.

In the paper, we studied the physical properties of the hydrogen molecule, which interacts with the environment according to the Balanced Gain and Loss (BGL) energy scheme³. The hydrogen molecule is an interesting case because it represents the non-trivial quantum system, and its state can be described accurately using the variational method^4,5,6,7. On the other hand, the BGL scheme describes the interaction between the molecule and the environment in the realistic and simple way. From the mathematical point of view, the BGL type interaction is modeled by the non-Hermitian Hamiltonian^8,9 which is invariant due to the ${\mathscr{PT}}$ symmetry (the symmetry of reflection in space (${\mathscr{P}}$) and in time (${\mathscr{T}}$))^{10,11,12,13,14}. It should be emphasized that if the Hamiltonian is the non-Hermitian one, but has the unbroken ${\mathscr{PT}}$ symmetry, the energy spectrum of the system is real – at least to the characteristic value of the parameter controlling the interaction with the environment.

The interest in the non-Hermitian Hamiltonians, in the context of the description of the open systems, appeared in many areas of physics. The papers^15,16, in which the open Bose-Hubbard dimer was analyzed are worth mentioning here. Such system can be implemented experimentally in the form of trapped bosons, where the coupling constant between the studied system and the environment reflects the value of barrier potential¹⁷. The ${\mathscr{PT}}$ symmetry breaking was also analyzed in^18,19,20 in the context of Bose-Einstein condensates. It is worth noting that the existence of the ${\mathscr{PT}}$ symmetry breaking effects were also confirmed in the field of quantum optics^21,22,23,24. Additionally, the fact of the appearance of complex energy values can be applied in explanation of the dynamics of the physical systems, the probability of disintegration, or the transport mechanism^{25,26,27,28,29,30,31}. It should be noticed that although most of the existing models were introduced heuristically, nevertheless it was done on the basis of relatively satisfactory mathematical justification¹⁷.

On the basis of the discussed issues, we intended to analyze the hydrogen molecule interacting with the environment and to examine its stability. We assumed that the required calculations would be carried out in the extremely accurate manner (at the level required by the the quantum chemistry standards), so that the obtained results could be verified experimentally.

In our opinion, the results presented in the paper shall be useful for people interested in the statistical and dynamical stability of small quantum systems interacting with environment and for those who research into the impact of ${\mathscr{PT}}$ symmetry breaking on the physical properties of a system. The presented results are significant insomuch that they concern the real physical system, and calculations were carried out for the model which does not contain any free parameters.

The achieved results can also be important from the technical point of view, because they can be related to the properties of the electronic devices in the single-molecule scale. Indeed, we already know that single particles can act similarly to the crucial elements of contemporary microelectronics, in particular they can serve as the rectifiers³², the electronic mixers³³, and the switchers^34,35,36. Therefore one can reasonably hope that the molecular electronics will replace the current technologies over time. However, the real development in this branch of knowledge depend on the full understanding of the transport mechanisms in single-molecule junctions. The presented work, by the example of the hydrogen molecule interacting with the environment in the balanced gain and loss energy scheme, draws attention to the fact that the interaction of the molecular bridge with anchors of the nanojunction can lead to changes in the bridge energy levels and to the reduction of their number. This is a substantial effect, because the electronic structure of a single molecule controls the electrical properties of the junction, in which it is used as a building block³⁷. The way of description of physical properties of the molecular bridge in the nanojunction, applied within the formalism presented in the work, is discussed in detail in the concluding part of the work.

Formalism

The total energy (E_T) of hydrogen molecule is defined as:

$${E}_{T}={E}_{p}+{E}_{e\gamma },$$

(1)

where: E_p = 2/R represents the energy of proton repulsion, with R = |R| as the distance between protons, and E_eγ means the energy of the lowest electronic state in the presence of the loss and gain effect (γ represents the coupling between the molecule and the environment). For γ = 0, the energy ${E}_{e}={E}_{e(\gamma \mathrm{=0})}$ should be determined using the Hubbard Hamiltonian, which takes into account all electronic interactions. In the second quantization formalism, we have⁶:

$$\begin{array}{rcl}{\hat{\mathscr{H}}}_{e} & = & \varepsilon ({\hat{n}}_{1}+{\hat{n}}_{2})+t\sum _{\sigma }\,({\hat{c}}_{1\sigma }^{\dagger }{\hat{c}}_{2\sigma }+{\hat{c}}_{2\sigma }^{\dagger }{\hat{c}}_{1\sigma })+U({\hat{n}}_{1\uparrow }{\hat{n}}_{1\downarrow }+{\hat{n}}_{2\uparrow }{\hat{n}}_{2\downarrow })\\ & & +\,(K-\frac{J}{2}){\hat{n}}_{1}{\hat{n}}_{2}-2J{\hat{{\bf{S}}}}_{1}{\hat{{\bf{S}}}}_{2}+J({\hat{c}}_{1\uparrow }^{\dagger }{\hat{c}}_{1\downarrow }^{\dagger }{\hat{c}}_{2\downarrow }{\hat{c}}_{2\uparrow }+h.\,c\mathrm{}.)\\ & & +\,V\sum _{\sigma }\,[({\hat{n}}_{1-\sigma }+{\hat{n}}_{2-\sigma })({\hat{c}}_{1\sigma }^{\dagger }{\hat{c}}_{2\sigma }+{\hat{c}}_{2\sigma }^{\dagger }{\hat{c}}_{1\sigma })],\end{array}$$

(2)

where ${\hat{n}}_{j}$ is given by: ${\hat{n}}_{j}={\sum }_{\sigma }{\hat{n}}_{j\sigma }={\sum }_{\sigma }{\hat{c}}_{j\sigma }^{\dagger }{\hat{c}}_{j\sigma }$, and ${\hat{c}}_{j\sigma }^{\dagger }$ (${\hat{c}}_{j\sigma }$) is the electron creation (annihilation) operator, which refers to the j-th hydrogen atom, σ represents the electronic spin: $\sigma \in \{\uparrow ,\downarrow \}$. In the last part of the Hamiltonian, ${\hat{\mathscr{H}}}_{e}$ the symbol −σ (in the subscript) denotes the spin direction opposite to the direction marked with σ. The product of spin operators ${\hat{{\bf{S}}}}_{i}{\hat{{\bf{S}}}}_{j}$ takes the form of: $\frac{1}{2}({\hat{S}}_{i}^{+}{\hat{S}}_{j}^{-}+{\hat{S}}_{i}^{-}{\hat{S}}_{j}^{+})+{\hat{S}}_{i}^{z}{\hat{S}}_{j}^{z}$, where ${\hat{S}}_{j}^{+}={\hat{c}}_{j\uparrow }^{\dagger }{\hat{c}}_{j\downarrow }$, ${\hat{S}}_{j}^{-}={\hat{c}}_{j\downarrow }^{\dagger }{\hat{c}}_{j\uparrow }$, and ${\hat{S}}_{j}^{z}=\frac{1}{2}({\hat{n}}_{j\uparrow }-{\hat{n}}_{j\downarrow })$. The symbol +h.c. in Eq. (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:

$$\begin{array}{rcl}\varepsilon & = & \int \,{d}^{3}{\bf{r}}{\Phi }_{1}({\bf{r}})[-{\nabla }^{2}-\frac{2}{|{\bf{r}}-{\bf{R}}|}]{\Phi }_{1}({\bf{r}}),\\ t & = & \int \,{d}^{3}{\bf{r}}{\Phi }_{1}({\bf{r}})[-{\nabla }^{2}-\frac{2}{|{\bf{r}}-{\bf{R}}|}]{\Phi }_{2}({\bf{r}}),\\ U & = & \int \int \,{d}^{3}{{\bf{r}}}_{1}{d}^{3}{{\bf{r}}}_{2}{\Phi }_{1}^{2}({{\bf{r}}}_{1})\frac{2}{|{{\bf{r}}}_{1}-{{\bf{r}}}_{2}|}{\Phi }_{1}^{2}({{\bf{r}}}_{2}),\\ K & = & \int \int \,{d}^{3}{{\bf{r}}}_{1}{d}^{3}{{\bf{r}}}_{2}{\Phi }_{1}^{2}({{\bf{r}}}_{1})\frac{2}{|{{\bf{r}}}_{1}-{{\bf{r}}}_{2}|}{\Phi }_{2}^{2}({{\bf{r}}}_{2}),\\ J & = & \int \int \,{d}^{3}{{\bf{r}}}_{1}{d}^{3}{{\bf{r}}}_{2}{\Phi }_{1}({{\bf{r}}}_{1}){\Phi }_{2}({{\bf{r}}}_{1})\frac{2}{|{{\bf{r}}}_{1}-{{\bf{r}}}_{2}|}{\Phi }_{1}({{\bf{r}}}_{2}){\Phi }_{2}({{\bf{r}}}_{2}),\\ V & = & \int \int \,{d}^{3}{{\bf{r}}}_{1}{d}^{3}{{\bf{r}}}_{2}{\Phi }_{1}^{2}({{\bf{r}}}_{1})\frac{2}{|{{\bf{r}}}_{1}-{{\bf{r}}}_{2}|}{\Phi }_{1}({{\bf{r}}}_{1}){\Phi }_{2}({{\bf{r}}}_{2})\mathrm{}.\end{array}$$

(3)

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:

$${\Phi }_{j}(r)=a[{\phi }_{j}({\bf{r}})-b{\phi }_{l}({\bf{r}})],$$

(4)

where the coefficients ensuring normalization are expressed in the formulas:

$$a=\frac{1}{\sqrt{2}}\sqrt{\frac{1+\sqrt{1-{S}^{2}}}{1-{S}^{2}}},\,b=\frac{S}{1+\sqrt{1-{S}^{2}}}\mathrm{}.$$

(5)

Table 1 The values of the Hubbard Hamiltonian integrals calculated for the equilibrium distance of the hydrogen molecule. The selected values γ have been taken into account.

The atomic overlap (S) has the form: $S=\int \,{d}^{3}r{\phi }_{1}({\bf{r}}){\phi }_{2}({\bf{r}})$, where 1 s Slater-type orbital can be written as: ${\phi }_{j}({\bf{r}})=\sqrt{{\alpha }^{3}/\pi }\,\exp [-\,\alpha |{\bf{r}}-{{\bf{R}}}_{j}|]$, α 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 ${\hat{\mathscr{H}}}_{e}$ with the balanced gain and loss operator^3,43:

$${\hat{\mathscr{H}}}_{\gamma }=i\gamma ({\hat{n}}_{1}-{\hat{n}}_{2})\mathrm{}.$$

(6)

The interpretation is that the operators $i\gamma {\hat{n}}_{1}$ and $-i\gamma {\hat{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 ${\hat{\mathscr{H}}}_{e\gamma }={\hat{\mathscr{H}}}_{e}+{\hat{\mathscr{H}}}_{\gamma }$ represents the non-Hermitian Hamiltonian, nevertheless it remains invariant due to the ${\mathscr{PT}}$ symmetry – at least to the characteristic ${\gamma }_{{\mathscr{PT}}}$ 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 ${\gamma }_{{\mathscr{PT}}}=0.520873$ Ry there occurs the breaking of ${\mathscr{PT}}$ symmetry of the electronic Hamiltonian. This fact is manifested by the appearance of the complex values of E₅ and E₆. Physically, this means that the ${\mathscr{PT}}$ 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}\rangle$ and $|{E}_{6}\rangle$ have the highest energy values. They cannot be thermally occupied - the k_BT energy is of the order of 25 meV, while the difference between E₆ and E₄ is around 25.5 eV (E₄ is the ground state energy of the electronic subsystem). When discussing the results, it should be clearly emphasized that E₄ is always the real number.

Figure 1

The real and the imaginary part of the eigenvalues of the Hamiltonian ${\hat{\mathscr{H}}}_{e\gamma }$. The equilibrium distances between protons (R₀) were assumed. The hatched areas correspond to the $\gamma$ values for which the operator ${\hat{\mathscr{H}}}_{e\gamma }$ ceases to be invariant due to the ${\mathscr{PT}}$ symmetry.

Although the ${\mathscr{PT}}$ 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 (${E}_{T}^{(j)}={E}_{p}+{E}_{j}$) of the isolated 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}_{T}^{(4)}({R}_{0})$ is observed with the increase in γ, whereas the molecule is in the stable state. Above ${\gamma }_{MS}=0.659374$ Ry the hydrogen molecule can exist only in the metastable state: ${E}_{T}^{(4)}({R}_{0}^{(MS)}) > {E}_{T}^{(4)}(R\to +\infty )=2$ Ry, where ${R}_{0}^{(MS)}=1.244701$ a₀. After exceeding γ_D = 1.024638 Ry, which corresponds to ${R}_{0}^{(D)}=1.196587$ a₀, the molecule breaks down (see the appendix). The insert (a) presents the dependence of the hydrogen dissociation energy (${E}_{{\rm{D}}}=2{\rm{R}}y-{E}_{T}^{(4)}$) on the value of the γ parameter. The insert (b) shows the influence of the γ parameter on the equilibrium distance R₀. Figure 3(a–c) trace the change of the distribution of electron charge for the stable case at γ = 0 (R0 = 1.41968 a₀), for the metastable case (${R}_{0}^{(MS)}$), and at the dissociation point (${R}_{0}^{(D)}$). The density of electron charge was calculated according to the formula: $\rho ({\bf{r}})={\sum }_{j}|{\Phi }_{j}({\bf{r}}){|}^{2}$.

Figure 2

The dependence of the total energy ${E}_{T}^{(j)}$ on the distance between protons. Additionaly, we took into account the influence of γ on the ground state energy ${E}_{T}^{(4)}$. Insert (a) the dissociation energy E_D versus γ parameter. Insert (b) equilibrium distance R₀ versus γ parameter.

Figure 3

The distribution of electronic charge in the hydrogen molecule: (a) the stable case; (b) metastable case; (c) just before dissociation.

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: ${V}_{{\rm{H}}}(R)={E}_{T}^{(4)}({R}_{0})+\frac{1}{2}{k}_{{\rm{H}}}{(R-{R}_{0})}^{2}$, where: ${k}_{{\rm{H}}}={[{d}^{2}{E}_{T}^{(4)}(R)/d{R}^{2}]}_{R={R}_{0}}$. The energy of quantum oscillator has the form:

$${E}_{o}^{{\rm{H}}}={\omega }_{0}^{{\rm{H}}}(n+\mathrm{1/2})\mathrm{}.$$

(7)

The symbol n indexes the energy level: $n=\mathrm{0,}\,\mathrm{1,}\,\mathrm{2,}\,\mathrm{..}.$. Additionally, ${\omega }_{0}^{{\rm{H}}}=\sqrt{{k}_{{\rm{H}}}/m{\prime} }$, where m′ is the reduced mass of the protons: $m{\prime} ={m}_{p}/2=918.076336$ (m_p is the proton mass). The more advanced approach is based on the Morse potential: ${V}_{{\rm{M}}o}(R)={E}_{T}^{(4)}({R}_{0})+{E}_{{\rm{D}}}{[1-\exp (-{\alpha }_{{\rm{M}}o}(R-{R}_{0}))]}^{2}$, where α_Mo means measure of the curvature of the potential about its minimum. The force constants, k_Mo should be calculated according to the formula: ${k}_{{\rm{Mo}}}={[{d}^{2}{V}_{{\rm{Mo}}}(R)/d{R}^{2}]}_{R={R}_{0}}$. The Morse energy is given by: ${\omega }_{0}^{{\rm{Mo}}}=\sqrt{{k}_{{\rm{M}}o}/m{\prime} }$ (see Table 2). The energy formula has the more complex form than for the harmonic case:

$${E}_{o}^{{\rm{Mo}}}={\omega }_{0}^{{\rm{Mo}}}(n+\mathrm{1/2})+({({\omega }_{0}^{{\rm{Mo}}})}^{2}\mathrm{/4}{E}_{D}){(n+\mathrm{1/2})}^{2}\mathrm{}.$$

(8)

Table 2 The Morse potential parameters for different values of γ.

Figure 4(a) depicts the dependence of the energies ${\omega }_{0}^{{\rm{H}}}$ and ${\omega }_{0}^{{\rm{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}_{T}^{(4)}(R)$ for higher values of γ.

Figure 4

(a) The influence of the γ parameter on the energy values ${\omega }_{0}^{{\rm{H}}}$ and ${\omega }_{0}^{{\rm{Mo}}}$. (b) The exemplary parameterization of the total energy curve in the harmonic and the anharmonic Morse case. (c) The rotational constant B₀ versus γ parameter.

The rotational energy of the hydrogen molecule should be calculated from the expression:

$${E}_{r}={B}_{0}l(l+1),$$

(9)

where: ${B}_{0}=\mathrm{1/}m{\prime} {R}_{0}^{2}$ and $l=0,\,1,\,2,\,\mathrm{..}.$. The influence of the γ parameter on the rotational energy value has been presented in the Fig. 4(c). From the physical point of view, the increase of the energy B₀ results from the decrease of the equilibrium distance R₀, which we observe when the γ parameter grows (see the Fig. 2 - insert (b)).

Having the explicit dependence of the ${\hat{\mathscr{H}}}_{e\gamma }$ parameters on R (see the appendix), we computed the electron-phonon coupling functions according to the formula: g_x = dx/dR, where $x\in \{\varepsilon ,t,U,K,J,V\}$. We plotted the obtained results in the Fig. 5. One can easily see that the absolute values of the considered functions at R₀ 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.

Figure 5

The electron-phonon coupling as a function of inter-proton distance for selected values of the γ parameter. The symbols placed on the curves point to the equilibrium value of the inter-proton distance.

The dynamic instability of the hydrogen molecule

The basic observable of the electronic subsystem is the occupation number $\langle {\hat{n}}_{j\sigma }\rangle$ of the j^th proton of the hydrogen molecule, where the symbol $\langle \mathrm{..}.\rangle$ means the expectation value. In the Hermitian case (γ = 0) the dynamics of $\langle {\hat{n}}_{j\sigma }\rangle$ can be analyzed using the conventional Heisenberg equation:

$$i\frac{d\langle {\hat{n}}_{j\sigma }\rangle }{dT}=\langle {[{\hat{n}}_{j\sigma },{\hat{\mathscr{H}}}_{e\gamma }^{MF}]}_{-}\rangle ,$$

(10)

where ${\hat{\mathscr{H}}}_{e\gamma }^{MF}$ is the electron Hamiltonian in the mean-field approximation:

$${\hat{\mathscr{H}}}_{e\gamma }^{MF}=\sum _{j\sigma }\,{\varepsilon }_{j\sigma }(\gamma ){\hat{n}}_{j\sigma }+\sum _{j\sigma }\,{t}_{j\sigma }{\hat{n}}_{j\overline{j}\sigma }+\sum _{j\sigma }\,{J}_{j\sigma }{\hat{n}}_{j\sigma -\sigma }+\sum _{j}({P}_{j}{\hat{\Delta }}_{j}^{\dagger }+{P}_{j}{\ast}{\hat{\Delta }}_{j})\mathrm{}.$$

(11)

The Hamiltonian parameters were defined by the expressions:

$$\begin{array}{rcl}{\varepsilon }_{j\sigma }(\gamma ) & = & \varepsilon +U\langle {\hat{n}}_{j-\sigma }\rangle +K\sum _{\sigma {\prime} }\,\langle {\hat{n}}_{\overline{j}\sigma {\prime} }\rangle -J\langle {\hat{n}}_{\overline{j}\sigma }\rangle \\ & & +\,V[\langle {\hat{n}}_{j\overline{j}-\sigma }\rangle +\langle {\hat{n}}_{\overline{j}j-\sigma }\rangle ]-{(-1)}^{j}i\gamma ,\\ {t}_{j\sigma } & = & t+V[\langle {\hat{n}}_{j-\sigma }\rangle +\langle {\hat{n}}_{\overline{j}-\sigma }\rangle ],\\ {J}_{j\sigma } & = & -J\langle {\hat{n}}_{\overline{j}-\sigma \sigma }\rangle ,\\ {P}_{j} & = & J\langle {\hat{\Delta }}_{\overline{j}}\rangle \mathrm{}.\end{array}$$

(12)

The new symbols have the following meanings:

$$\bar{j}=\{\begin{array}{lll}1 & {\rm{for}} & j=2\\ 2 & {\rm{for}} & j=\mathrm{1,}\end{array}$$

(13)

${\hat{n}}_{j\overline{j}\sigma }={\hat{c}}_{j\sigma }^{\dagger }{\hat{c}}_{\overline{j}\sigma }$, ${\hat{n}}_{j\sigma -\sigma }={\hat{c}}_{j\sigma }^{\dagger }{\hat{c}}_{j-\sigma }$, and ${\hat{\Delta }}_{j}={\hat{c}}_{j\downarrow }{\hat{c}}_{j\uparrow }$.

It should be emphasized that for ${\hat{\mathscr{H}}}_{e\gamma }$ the required operator calculations are not feasible due to their extensiveness. The mean-field approximation transforms the operator ${\hat{\mathscr{H}}}_{e\gamma }$ 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 ${\hat{\Delta }}_{j}$, which has the formal structure of the Cooper pair annihilation operator in the real space. This analogy is not complete, because the Hamiltonian ${\hat{\mathscr{H}}}_{e\gamma }^{MF}$ term containing ${\hat{\Delta }}_{j}$ and ${\hat{\Delta }}_{j}^{\dagger }$ does not correspond to BCS pairing operator^48,49,50 (the integral of the exchange J₀ 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: ${\hat{\mathscr{H}}}_{e\gamma \pm }^{MF}=\frac{1}{2}({\hat{\mathscr{H}}}_{e\gamma }^{MF}\pm {\hat{\mathscr{H}}}_{e\gamma }^{MF\dagger })$, where ${\hat{\mathscr{H}}}_{e\gamma \pm }^{MF}=\pm \,{\hat{\mathscr{H}}}_{e\gamma \pm }^{MF\dagger }$. Then we use the generalized form of the Heisenberg equation:

$$i\frac{d\langle {\hat{n}}_{j\sigma }\rangle }{dT}=\langle {[{\hat{n}}_{j\sigma },{\hat{\mathscr{H}}}_{e\gamma +}^{MF}]}_{-}\rangle +\langle {[{\hat{n}}_{j\sigma },{\hat{\mathscr{H}}}_{e\gamma -}^{MF}]}_{+}\rangle -2\langle {\hat{n}}_{j\sigma }\rangle \langle {\hat{\mathscr{H}}}_{e\gamma -}^{MF}\rangle \mathrm{}.$$

(14)

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 $\langle {n}_{1\uparrow }\rangle$, $\langle {n}_{2\uparrow }\rangle$ 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 $\langle {n}_{1\uparrow }\rangle$, $\langle {n}_{2\uparrow }\rangle$ 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.

Figure 6

(a–e) The time evolution of $\langle {n}_{1\uparrow }\rangle$ and $\langle {n}_{2\uparrow }\rangle$, respectively, for γ equal to: 0, 0.1 Ry, ${\gamma }_{{\mathscr{PT}}}$, ${\gamma }_{{\mathscr{MS}}}$, and 1 Ry. The hatched areas correspond to the value greater than T_D. T_D is a time, after which an imaginary part of the observables exceeds the value of ±0.005. Figure (f) shows the form of T_D(γ) function. The dashed curve was obtained from the formula: ${T}_{D}=\frac{a}{{\gamma }^{b}}\exp (c\gamma )$, where a = 0.095689, b = 2.38294, and c = 2.80013.

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 ${\gamma }_{MS}=0.659374$ Ry, the molecule is in the metastable state, decaying definitively for ${\gamma }_{D} > 1.024638$ Ry.

An additional effect, that we observed for γ higher than ${\gamma }_{{\mathscr{PT}}}=0.520873$ Ry, is the ${\mathscr{PT}}$ symmetry breaking of the electronic Hamiltonian ${\hat{\mathscr{H}}}_{e\gamma }$. 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 ${\mathscr{PT}}$ 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 ${\mathscr{PT}}$ symmetry of ${\hat{\mathscr{H}}}_{e\gamma }$.

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₂/X, where symbol X means metals like Pt^{52,53,56,57,58}, Pd⁵⁵, Au^58,59, Cu⁶⁰ or Ni⁶¹. 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 ($\gamma > {\gamma }_{{\mathscr{PT}}}^{{\prime} }$, whereas ${\gamma }_{{\mathscr{PT}}}^{{\prime} }$ means the value of γ parameter for which the ${\mathscr{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{\prime} }_{0}$), should be done according to the method based on the density-functional theory (DFT)⁶². The physical state of the bridge when there is no flux (γ = 0) corresponds to the minimum enthalpy value: $H={E}_{p}+{E}_{e}+FR{\prime}$, 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 achieved for the ${R{\prime} }_{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: $H={E}_{p}+{E}_{e\gamma }+FR{\prime}$.

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.

The eigenvalues of the electronic Hamiltonian of the hydrogen molecule interacting with the environment

The Hamiltonian ${\hat{\mathscr{H}}}_{e\gamma }$ should be written in the matrix form:

$$(\begin{array}{cccccc}{h}_{1}^{+} & 0 & t+V & 0 & J & t+V\\ 0 & {h}_{2} & 0 & 0 & 0 & 0\\ t+V & 0 & 2\varepsilon +K & 0 & t+V & -\,J\\ 0 & 0 & 0 & {h}_{2} & 0 & 0\\ J & 0 & t+V & 0 & {h}_{1}^{-} & t+V\\ t+V & 0 & -\,J & 0 & t+V & 2\varepsilon +K\end{array}).$$

(15)

where: ${h}_{1}^{\pm }=2\varepsilon +U\pm 2i\gamma$, and ${h}_{2}=2\varepsilon +K-J$. By using the operator (15) there was brought out the preliminary formulas for the eigenvalues, which has the form as follows:

$${E}_{1}=-\,J+K+2\varepsilon ,$$

(16)

$${E}_{2}=-\,J+K+2\varepsilon ,$$

(17)

$${E}_{3}=J+K+2\varepsilon ,$$

(18)

$$\begin{array}{rcl}{E}_{4} & = & \frac{1}{3}(\,-\,J+K+2U+[\,-\,4{J}^{2}+2J(K-U)\\ & & -\,{(K-U)}^{2}+12(\,-\,{(t+V)}^{2}+{\gamma }^{2})]/[\,-\,8{J}^{3}+6{J}^{2}(K-U)\\ & & +\,3J({(K-U)}^{2}+12(\,-{(t+V)}^{2}+{\gamma }^{2}))-(K-U)({(K-U)}^{2}\\ & & +\,\mathrm{18((}t+V{)}^{2}+2{\gamma }^{2}))+\frac{1}{2}\sqrt{A-B}{]}^{\frac{1}{3}}\\ & & -\,[\,-\,8{J}^{3}+6{J}^{2}(K-U)+3J((K-U{)}^{2}\\ & & +\,\mathrm{12(}\,-\,{(t+V)}^{2}+{\gamma }^{2}))-(K-U)((K-U{)}^{2}\\ & & +\,\mathrm{18((}t+V{)}^{2}+2{\gamma }^{2}))+\frac{1}{2}\sqrt{A-B}{]}^{\frac{1}{3}}+6\varepsilon ),\end{array}$$

(19)

$$\begin{array}{rcl}{E}_{5} & = & \frac{1}{12}\mathrm{([2(1}+i\sqrt{3}\mathrm{)(4}{J}^{2}+{(K-U)}^{2}+2J(\,-\,K+U)\\ & & +\,\mathrm{12(}t+V-\gamma )(t+V+\gamma ))]/[\,-\,8{J}^{3}+6{J}^{2}(K-U)\\ & & +\,3J((K-U{)}^{2}+\mathrm{12(}\,-\,{(t+V)}^{2}+{\gamma }^{2}))-(K-U)((K-U{)}^{2}\\ & & +\,\mathrm{18((}t+V{)}^{2}+2{\gamma }^{2}))+\frac{1}{2}\sqrt{A-B}{]}^{\frac{1}{3}}+\mathrm{2(1}-i\sqrt{3})[\,-\,8{J}^{3}\\ & & +\,6{J}^{2}(K-U)+3J((K-U{)}^{2}+12(\,-\,{(t+V)}^{2}+{\gamma }^{2}))\\ & & -\,(K-U)((K-U{)}^{2}+\mathrm{18((}t+V{)}^{2}+2{\gamma }^{2}))+\frac{1}{2}\sqrt{A-B}{]}^{\frac{1}{3}}\\ & & +\,\mathrm{4(}\,-\,J+K+2U+6\varepsilon )),\end{array}$$

(20)

$$\begin{array}{rcl}{E}_{6} & = & \frac{1}{12}\mathrm{([2(1}-i\sqrt{3}\mathrm{)(4}{J}^{2}+{(K-U)}^{2}+2J(\,-\,K+U)\\ & & +\,\mathrm{12(}t+V-\gamma )(t+V+\gamma ))]/[-\,8{J}^{3}+6{J}^{2}(K-U)\\ & & +\,3J((K-U{)}^{2}+\mathrm{12(}\,-\,{(t+V)}^{2}+{\gamma }^{2}))-(K-U)((K-U{)}^{2}\\ & & +\,\mathrm{18((}t+V{)}^{2}+2{\gamma }^{2}))+\frac{1}{2}\sqrt{A-B}{]}^{\frac{1}{3}}+\mathrm{2(1}+i\sqrt{3})[\,-\,8{J}^{3}\\ & & +\,6{J}^{2}(K-U)+3J((K-U{)}^{2}+\mathrm{12(}\,-\,{(t+V)}^{2}+{\gamma }^{2}))\\ & & -\,(K-U)((K-U{)}^{2}+\mathrm{18((}t+V{)}^{2}+2{\gamma }^{2}))+\frac{1}{2}\sqrt{A-B}{]}^{\frac{1}{3}}\\ & & +\,\mathrm{4(}\,-\,J+K+2U+6\varepsilon )),\end{array}$$

(21)

wherein:

$$A=4{[(2J+K-U)[(J-K+U)(4J-K+U)+18{(t+V)}^{2}]-36(J-K+U){\gamma }^{2}]}^{2},$$

(22)

and

$$B=4{[4{J}^{2}+{(K-U)}^{2}+2J(-K+U)+12(t+V-\gamma )(t+V+\gamma )]}^{3}.$$

(23)

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 ${\hat{\mathscr{H}}}_{e\gamma }$ parameters. The explicit dependence of the variational parameter α on the distance R is shown in the Fig. 8.

Figure 7

The integral of the Hamiltonian ${\hat{\mathscr{H}}}_{e\gamma }$ as a function of the inter-proton distance for the selected values of the parameter modeling the interaction of the molecule with the environment. The balls placed on the curves point to the equilibrium value of the inter-proton distance R₀.

Figure 8

The variation parameter α as a function of the proton distance for selected γ values.

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₀, the equilibrium inverse size of the orbital α₀, and the ground-state energy ET(4)(R₀) for different values of γ.

Table 3 The harmonic potential parameter k_H and the quantum energy for different values of γ.

Table 4 The values of the electron-ion coupling constants at the hydrogen-molecule equilibrium for different values of γ.

Table 5 The equilibrium distance R₀, the equilibrium inverse size of the orbital α₀, and the ground-state energy ${E}_{T}^{(4)}({R}_{0})$ for different values of γ.

The set of the differential equations for electron observables (γ = 0)

The system of differential equations has the form:

$$\begin{array}{ccc}i\frac{d\langle {\hat{n}}_{1\uparrow }\rangle }{dT} & = & {t}_{1\uparrow }\langle {\hat{n}}_{12\uparrow }\rangle -{t}_{2\uparrow }\langle {\hat{n}}_{21\uparrow }\rangle +{J}_{1\uparrow }\langle {\hat{n}}_{1\uparrow \downarrow }\rangle \\ & & -\,{J}_{1\downarrow }\langle {\hat{n}}_{1\downarrow \uparrow }\rangle +{P}_{1}\langle {\hat{\Delta }}_{1}^{\dagger }\rangle -{P}_{1}^{\ast }\langle {\hat{\Delta }}_{1}\rangle ,\end{array}$$

(24)

$$\begin{array}{ccc}i\frac{d\langle {\hat{n}}_{1\downarrow }\rangle }{dT} & = & {t}_{1\downarrow }\langle {\hat{n}}_{12\downarrow }\rangle -{t}_{2\downarrow }\langle {\hat{n}}_{21\downarrow }\rangle +{J}_{1\downarrow }\langle {\hat{n}}_{1\downarrow \uparrow }\rangle \\ & & -\,{J}_{1\uparrow }\langle {\hat{n}}_{1\uparrow \downarrow }\rangle +{P}_{1}\langle {\hat{\Delta }}_{1}^{\dagger }\rangle -{P}_{1}^{\ast }\langle {\hat{\Delta }}_{1}\rangle ,\end{array}$$

(25)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{2\uparrow }\rangle }{dT} & = & -{t}_{1\uparrow }\langle {\hat{n}}_{12\uparrow }\rangle +{t}_{2\uparrow }\langle {\hat{n}}_{21\uparrow }\rangle +{J}_{2\uparrow }\langle {\hat{n}}_{2\uparrow \downarrow }\rangle -{J}_{2\downarrow }\langle {\hat{n}}_{2\downarrow \uparrow }\rangle \\ & & +\,{P}_{2}\langle {\hat{\Delta }}_{2}^{\dagger }\rangle -{P}_{2}^{\ast }\langle {\hat{\Delta }}_{2}\rangle ,\end{array}$$

(26)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{2\downarrow }\rangle }{dT} & = & -{t}_{1\downarrow }\langle {\hat{n}}_{12\downarrow }\rangle +{t}_{2\downarrow }\langle {\hat{n}}_{21\downarrow }\rangle +{J}_{2\downarrow }\langle {\hat{n}}_{2\downarrow \uparrow }\rangle \\ & & -\,{J}_{2\uparrow }\langle {\hat{n}}_{2\uparrow \downarrow }\rangle +{P}_{2}\langle {\hat{\Delta }}_{2}^{\dagger }\rangle -{P}_{2}^{\ast }\langle {\hat{\Delta }}_{2}\rangle ,\end{array}$$

(27)

$$i\frac{d\langle {\hat{n}}_{12\uparrow }\rangle }{dT}=-\,{\varepsilon }_{1\uparrow }\langle {\hat{n}}_{12\uparrow }\rangle +{\varepsilon }_{2\uparrow }\langle {\hat{n}}_{12\uparrow }\rangle +{t}_{2\uparrow }\langle {\hat{n}}_{1\uparrow }\rangle -{t}_{2\uparrow }\langle {\hat{n}}_{2\uparrow }\rangle ,$$

(28)

$$i\frac{d\langle {\hat{n}}_{12\downarrow }\rangle }{dT}=-\,{\varepsilon }_{1\downarrow }\langle {\hat{n}}_{12\downarrow }\rangle +{\varepsilon }_{2\downarrow }\langle {\hat{n}}_{12\downarrow }\rangle +{t}_{2\downarrow }\langle {\hat{n}}_{1\downarrow }\rangle -{t}_{2\downarrow }\langle {\hat{n}}_{2\downarrow }\rangle ,$$

(29)

$$i\frac{d\langle {\hat{n}}_{21\uparrow }\rangle }{dT}={\varepsilon }_{1\uparrow }\langle {\hat{n}}_{21\uparrow }\rangle -{\varepsilon }_{2\uparrow }\langle {\hat{n}}_{21\uparrow }\rangle +{t}_{1\uparrow }\langle {\hat{n}}_{2\uparrow }\rangle -{t}_{1\uparrow }\langle {\hat{n}}_{1\uparrow }\rangle ,$$

(30)

$$i\frac{d\langle {\hat{n}}_{21\downarrow }\rangle }{dT}={\varepsilon }_{1\downarrow }\langle {\hat{n}}_{21\downarrow }\rangle -{\varepsilon }_{2\downarrow }\langle {\hat{n}}_{21\downarrow }\rangle +{t}_{1\downarrow }\langle {\hat{n}}_{2\downarrow }\rangle -{t}_{1\downarrow }\langle {\hat{n}}_{1\downarrow }\rangle ,$$

(31)

$$i\frac{d\langle {\hat{n}}_{1\uparrow \downarrow }\rangle }{dT}=-\,{\varepsilon }_{1\uparrow }\langle {\hat{n}}_{1\uparrow \downarrow }\rangle +{\varepsilon }_{1\downarrow }\langle {\hat{n}}_{1\uparrow \downarrow }\rangle +{J}_{1\downarrow }\langle {\hat{n}}_{1\uparrow }\rangle -{J}_{1\downarrow }\langle {\hat{n}}_{1\downarrow }\rangle ,$$

(32)

$$i\frac{d\langle {\hat{n}}_{1\downarrow \uparrow }\rangle }{dT}={\varepsilon }_{1\uparrow }\langle {\hat{n}}_{1\downarrow \uparrow }\rangle -{\varepsilon }_{1\downarrow }\langle {\hat{n}}_{1\downarrow \uparrow }\rangle +{J}_{1\uparrow }\langle {\hat{n}}_{1\downarrow }\rangle -{J}_{1\uparrow }\langle {\hat{n}}_{1\uparrow }\rangle ,$$

(33)

$$i\frac{d\langle {\hat{n}}_{2\uparrow \downarrow }\rangle }{dT}=-\,{\varepsilon }_{2\uparrow }\langle {\hat{n}}_{2\uparrow \downarrow }\rangle +{\varepsilon }_{2\downarrow }\langle {\hat{n}}_{2\uparrow \downarrow }\rangle +{J}_{2\downarrow }\langle {\hat{n}}_{2\uparrow }\rangle -{J}_{2\downarrow }\langle {\hat{n}}_{2\downarrow }\rangle ,$$

(34)

$$i\frac{d\langle {\hat{n}}_{2\downarrow \uparrow }\rangle }{dT}={\varepsilon }_{2\uparrow }\langle {\hat{n}}_{2\downarrow \uparrow }\rangle -{\varepsilon }_{2\downarrow }\langle {\hat{n}}_{2\downarrow \uparrow }\rangle +{J}_{2\uparrow }\langle {\hat{n}}_{2\downarrow }\rangle -{J}_{2\uparrow }\langle {\hat{n}}_{2\uparrow }\rangle ,$$

(35)

$$i\frac{d\langle {\hat{\Delta }}_{1}^{\dagger }\rangle }{dT}=-\,{\varepsilon }_{1\uparrow }\langle {\hat{\Delta }}_{1}^{\dagger }\rangle -{\varepsilon }_{1\downarrow }\langle {\hat{\Delta }}_{1}^{\dagger }\rangle +{P}_{1}^{\ast }\langle {\hat{n}}_{1\uparrow }\rangle +{P}_{1}^{\ast }\langle {\hat{n}}_{1\downarrow }\rangle -{P}_{1}^{\ast },$$

(36)

$$i\frac{d\langle {\hat{\Delta }}_{1}\rangle }{dT}={\varepsilon }_{1\uparrow }\langle {\hat{\Delta }}_{1}\rangle +{\varepsilon }_{1\downarrow }\langle {\hat{\Delta }}_{1}\rangle -{P}_{1}\langle {\hat{n}}_{1\downarrow }\rangle -{P}_{1}\langle {\hat{n}}_{1\uparrow }\rangle +{P}_{1},$$

(37)

$$i\frac{d\langle {\hat{\Delta }}_{2}^{\dagger }\rangle }{dT}=-\,{\varepsilon }_{2\uparrow }\langle {\hat{\Delta }}_{2}^{\dagger }\rangle -{\varepsilon }_{2\downarrow }\langle {\hat{\Delta }}_{2}^{\dagger }\rangle +{P}_{2}^{\ast }\langle {\hat{n}}_{2\uparrow }\rangle +{P}_{2}^{\ast }\langle {\hat{n}}_{2\downarrow }\rangle -{P}_{2}^{\ast },$$

(38)

$$i\frac{d\langle {\hat{\Delta }}_{2}\rangle }{dT}={\varepsilon }_{2\uparrow }\langle {\hat{\Delta }}_{2}\rangle +{\varepsilon }_{2\downarrow }\langle {\hat{\Delta }}_{2}\rangle -{P}_{2}\langle {\hat{n}}_{2\downarrow }\rangle -{P}_{2}\langle {\hat{n}}_{2\uparrow }\rangle +{P}_{2}.$$

(39)

The system of differential equations for electron observables (γ ≠ 0)

The system of differential equations can be written in the form:

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{1\uparrow }\rangle }{dT} & = & {t}_{1\uparrow }\langle {\hat{n}}_{12\uparrow }\rangle -{t}_{2\uparrow }\langle {\hat{n}}_{21\uparrow }\rangle +{J}_{1\uparrow }\langle {\hat{n}}_{1\uparrow \downarrow }\rangle -{J}_{1\downarrow }\langle {\hat{n}}_{1\downarrow \uparrow }\rangle \\ & & +\,{P}_{1}\langle {\hat{\Delta }}_{1}^{\dagger }\rangle -{P}_{1}^{\ast }\langle {\hat{\Delta }}_{1}\rangle +i\gamma \langle {\hat{n}}_{1\uparrow }\rangle -2i\gamma |{\Delta }_{1}{|}^{2}\\ & & -\,2i\gamma < {\hat{n}}_{1\uparrow } > \sum _{\sigma }(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(40)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{1\downarrow }\rangle }{dT} & = & {t}_{1\downarrow }\langle {\hat{n}}_{12\downarrow }\rangle -{t}_{2\downarrow }\langle {\hat{n}}_{21\downarrow }\rangle +{J}_{1\downarrow }\langle {\hat{n}}_{1\downarrow \uparrow }\rangle -{J}_{1\uparrow }\langle {\hat{n}}_{1\uparrow \downarrow }\rangle \\ & & +\,i\gamma \langle {\hat{n}}_{1\downarrow }\rangle -2i\gamma |{\Delta }_{1}{|}^{2}+{P}_{1}\langle {\hat{\Delta }}_{1}^{\dagger }\rangle -{P}_{1}^{\ast }\langle {\hat{\Delta }}_{1}\rangle \\ & & -\,2i\gamma \langle {\hat{n}}_{1\downarrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(41)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{2\uparrow }\rangle }{dT} & = & -{t}_{1\uparrow }\langle {\hat{n}}_{12\uparrow }\rangle +{t}_{2\uparrow }\langle {\hat{n}}_{21\uparrow }\rangle +{J}_{2\uparrow }\langle {\hat{n}}_{2\uparrow \downarrow }\rangle -{J}_{2\downarrow }\langle {\hat{n}}_{2\downarrow \uparrow }\rangle \\ & & +\,{P}_{2}\langle {\hat{\Delta }}_{2}^{\dagger }\rangle -{P}_{2}^{\ast }\langle {\hat{\Delta }}_{2}\rangle -i\gamma \langle {\hat{n}}_{2\uparrow }\rangle +2i\gamma |{\Delta }_{2}{|}^{2}\\ & & -\,2i\gamma \langle {\hat{n}}_{2\uparrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(42)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{2\downarrow }\rangle }{dT} & = & -{t}_{1\downarrow }\langle {\hat{n}}_{12\downarrow }\rangle +{t}_{2\downarrow }\langle {\hat{n}}_{21\downarrow }\rangle +{J}_{2\downarrow }\langle {\hat{n}}_{2\downarrow \uparrow }\rangle -{J}_{2\uparrow }\langle {\hat{n}}_{2\uparrow \downarrow }\rangle \\ & & +\,{P}_{2}\langle {\hat{\Delta }}_{2}^{\dagger }\rangle -{P}_{2}^{\ast }\langle {\hat{\Delta }}_{2}\rangle -i\gamma \langle {\hat{n}}_{2\downarrow }\rangle +2i\gamma |{\Delta }_{2}{|}^{2}\\ & & -\,2i\gamma \langle {\hat{n}}_{2\downarrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(43)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{12\uparrow }\rangle }{dT} & = & -{\varepsilon }_{1\uparrow }\langle {\hat{n}}_{12\uparrow }\rangle +{\varepsilon }_{2\uparrow }\langle {\hat{n}}_{12\uparrow }\rangle +{t}_{2\uparrow }\langle {\hat{n}}_{1\uparrow }\rangle -{t}_{2\uparrow }\langle {\hat{n}}_{2\uparrow }\rangle \\ & & +\,i\gamma \langle {\hat{n}}_{12\uparrow }\rangle \\ & & -\,2i\gamma \langle {\hat{n}}_{12\uparrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(44)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{12\downarrow }\rangle }{dT} & = & -{\varepsilon }_{1\downarrow }\langle {\hat{n}}_{12\downarrow }\rangle +{\varepsilon }_{2\downarrow }\langle {\hat{n}}_{12\downarrow }\rangle +{t}_{2\downarrow }\langle {\hat{n}}_{1\downarrow }\rangle -{t}_{2\downarrow }\langle {\hat{n}}_{2\downarrow }\rangle \\ & & +\,i\gamma \langle {\hat{n}}_{12\downarrow }\rangle \\ & & -\,2i\gamma \langle {\hat{n}}_{12\downarrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(45)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{21\uparrow }\rangle }{dT} & = & {\varepsilon }_{1\uparrow }\langle {\hat{n}}_{21\uparrow }\rangle -{\varepsilon }_{2\uparrow }\langle {\hat{n}}_{21\uparrow }\rangle +{t}_{1\uparrow }\langle {\hat{n}}_{2\uparrow }\rangle -{t}_{1\uparrow }\langle {\hat{n}}_{1\uparrow }\rangle \\ & & -\,i\gamma \langle {\hat{n}}_{21\uparrow }\rangle \\ & & -\,2i\gamma \langle {\hat{n}}_{21\uparrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(46)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{21\downarrow }\rangle }{dT} & = & {\varepsilon }_{1\downarrow }\langle {\hat{n}}_{21\downarrow }\rangle -{\varepsilon }_{2\downarrow }\langle {\hat{n}}_{21\downarrow }\rangle +{t}_{1\downarrow }\langle {\hat{n}}_{2\downarrow }\rangle -{t}_{1\downarrow }\langle {\hat{n}}_{1\downarrow }\rangle \\ & & -\,i\gamma \langle {\hat{n}}_{21\downarrow }\rangle \\ & & -\,2i\gamma \langle {\hat{n}}_{21\downarrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(47)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{1\uparrow \downarrow }\rangle }{dT} & = & -{\varepsilon }_{1\uparrow }\langle {\hat{n}}_{1\uparrow \downarrow }\rangle +{\varepsilon }_{1\downarrow }\langle {\hat{n}}_{1\uparrow \downarrow }\rangle +{J}_{1\downarrow }\langle {\hat{n}}_{1\uparrow }\rangle -{J}_{1\downarrow }\langle {\hat{n}}_{1\downarrow }\rangle \\ & & +\,i\gamma \langle {\hat{n}}_{1\uparrow \downarrow }\rangle \\ & & -\,2i\gamma \langle {\hat{n}}_{1\uparrow \downarrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(48)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{1\downarrow \uparrow }\rangle }{dT} & = & {\varepsilon }_{1\uparrow }\langle {\hat{n}}_{1\downarrow \uparrow }\rangle -{\varepsilon }_{1\downarrow }\langle {\hat{n}}_{1\downarrow \uparrow }\rangle +{J}_{1\uparrow }\langle {\hat{n}}_{1\downarrow }\rangle -{J}_{1\uparrow }\langle {\hat{n}}_{1\uparrow }\rangle \\ & & +\,i\gamma \langle {\hat{n}}_{1\downarrow \uparrow }\rangle \\ & & -\,2i\gamma \langle {\hat{n}}_{1\downarrow \uparrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(49)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{2\uparrow \downarrow }\rangle }{dT} & = & -{\varepsilon }_{2\uparrow }\langle {\hat{n}}_{2\uparrow \downarrow }\rangle +{\varepsilon }_{2\downarrow }\langle {\hat{n}}_{2\uparrow \downarrow }\rangle +{J}_{2\downarrow }\langle {\hat{n}}_{2\uparrow }\rangle -{J}_{2\downarrow }\langle {\hat{n}}_{2\downarrow }\rangle \\ & & -\,i\gamma \langle {\hat{n}}_{2\uparrow \downarrow }\rangle \\ & & -\,2i\gamma \langle {\hat{n}}_{2\uparrow \downarrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(50)

$$\begin{array}{rcl}i\frac{d\langle {\hat{n}}_{2\downarrow \uparrow }\rangle }{dT} & = & {\varepsilon }_{2\uparrow }\langle {\hat{n}}_{2\downarrow \uparrow }\rangle -{\varepsilon }_{2\downarrow }\langle {\hat{n}}_{2\downarrow \uparrow }\rangle +{J}_{2\uparrow }\langle {\hat{n}}_{2\downarrow }\rangle -{J}_{2\uparrow }\langle {\hat{n}}_{2\uparrow }\rangle \\ & & -\,i\gamma \langle {\hat{n}}_{2\downarrow \uparrow }\rangle \\ & & -\,2i\gamma \langle {\hat{n}}_{2\downarrow \uparrow }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(51)

$$\begin{array}{rcl}i\frac{d\langle {\hat{\Delta }}_{1}^{\dagger }\rangle }{dT} & = & -{\varepsilon }_{1\uparrow }\langle {\hat{\Delta }}_{1}^{\dagger }\rangle -{\varepsilon }_{1\downarrow }\langle {\hat{\Delta }}_{1}^{\dagger }\rangle +{P}_{1}^{\ast }\langle {\hat{n}}_{1\uparrow }\rangle +{P}_{1}^{\ast }\langle {\hat{n}}_{1\downarrow }\rangle -{P}_{1}^{\ast }\\ & & +\,2i\gamma \langle {\hat{\Delta }}_{1}^{\dagger }\rangle +2i\gamma \langle {\hat{\Delta }}_{1}^{\dagger }\rangle \langle {\hat{n}}_{1\downarrow }\rangle \\ & & -\,2i\gamma \langle {\hat{\Delta }}_{1}^{\dagger }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(52)

$$\begin{array}{rcl}i\frac{d\langle {\hat{\Delta }}_{1}\rangle }{dT} & = & {\varepsilon }_{1\uparrow }\langle {\hat{\Delta }}_{1}\rangle +{\varepsilon }_{1\downarrow }\langle {\hat{\Delta }}_{1}\rangle -{P}_{1}\langle {\hat{n}}_{1\downarrow }\rangle -{P}_{1}\langle {\hat{n}}_{1\uparrow }\rangle +{P}_{1}\\ & & -\,i\gamma \langle {\hat{\Delta }}_{1}\rangle +2i\gamma \langle {\hat{\Delta }}_{1}\rangle \langle {\hat{n}}_{1\downarrow }\rangle \\ & & -\,2i\gamma \langle {\hat{\Delta }}_{1}\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(53)

$$\begin{array}{rcl}i\frac{d\langle {\hat{\Delta }}_{2}^{\dagger }\rangle }{dT} & = & -{\varepsilon }_{2\uparrow }\langle {\hat{\Delta }}_{2}^{\dagger }\rangle -{\varepsilon }_{2\downarrow }\langle {\hat{\Delta }}_{2}^{\dagger }\rangle +{P}_{2}^{\ast }\langle {\hat{n}}_{2\uparrow }\rangle +{P}_{2}^{\ast }\langle {\hat{n}}_{2\downarrow }\rangle -{P}_{2}^{\ast }\\ & & -\,2i\gamma \langle {\hat{\Delta }}_{2}^{\dagger }\rangle -2i\gamma \langle {\hat{\Delta }}_{2}^{\dagger }\rangle \langle {\hat{n}}_{2\downarrow }\rangle \\ & & -\,2i\gamma \langle {\hat{\Delta }}_{2}^{\dagger }\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle ),\end{array}$$

(54)

$$\begin{array}{rcl}i\frac{d\langle {\hat{\Delta }}_{2}\rangle }{dT} & = & {\varepsilon }_{2\uparrow }\langle {\hat{\Delta }}_{2}\rangle +{\varepsilon }_{2\downarrow }\langle {\hat{\Delta }}_{2}\rangle -{P}_{2}\langle {\hat{n}}_{2\downarrow }\rangle -{P}_{2}\langle {\hat{n}}_{2\uparrow }\rangle +{P}_{2}\\ & & +\,i\gamma \langle {\hat{\Delta }}_{2}\rangle -2i\gamma \langle {\hat{\Delta }}_{2}\rangle \langle {\hat{n}}_{2\downarrow }\rangle \\ & & -\,2i\gamma \langle {\hat{\Delta }}_{2}\rangle \sum _{\sigma }\,(\langle {\hat{n}}_{1\sigma }\rangle -\langle {\hat{n}}_{2\sigma }\rangle )\mathrm{}.\end{array}$$

(55)