format: text

text:
# Josephson current through a correlated quantum level: Andreev states and 𝜋 junction behavior

## I Introduction

The observation of the Kondo effect in semiconducting quantum dots  has opened a new area of research in which electronic transport through a nanoscale strongly correlated system can be studied under controlled conditions. More recently, this effect has also been observed in carbon nanotubes coupled to metallic electrodes, a system that can behave in many respects as a quantum dot . These type of devices provide an almost ideal system to test different theoretical predictions. In this direction, a great theoretical interest has arisen in connection with the possibility of replacing the normal electrodes by superconducting ones. In fact, carbon nanotubes connected to superconducting leads, that are already being produced, can provide a physical realization of such a system .

In these type of systems an issue of fundamental relevance is the interplay between electron correlation effects and Andreev reflection processes, which provide the basic mechanism for transport between normal and superconducting regions. This interplay has been analyzed by several authors for the case of quantum dot coupled to a normal and a superconducting electrode . There have also been some works addressing the problem of the electron transport mediated by multiple Andreev reflections through a resonant level between two superconducting electrodes .

When both electrodes are superconducting a more basic question is how the electron-electron interactions in the dot would affect the Josephson current. This issue has received considerable attention in recent years . In particular, the discussion has been centered to a large extent around the appearance of a $`\pi `$-junction behavior induced by electron-electron interactions. The $`\pi `$-junction behavior, which consists in a reversal of the sign of the supercurrent under certain conditions, was first pointed out by Kulik  when analyzing the Josephson tunneling in the presence of magnetic impurities, and discussed afterwards by several authors . This issue is also related to the interplay between superconductivity and magnetism. $`\pi `$-junction behavior has in fact been recently observed in S-F-S Josephson junctions .

The existing theoretical works analyzing the Josephson effect in a single correlated level coupled to superconducting leads have adopted either a mean field description  or are restricted to some limiting situations . In Ref.  the current was obtained to the lowest order in the tunneling coupling, thereby neglecting the important physics associated with the Andreev bound states. On the other hand, in Refs.  the limit of infinite Coulomb repulsion was considered.

At this point we believe that further work is needed to understand the physical behavior of this system for a broad parameter range. This behavior is actually rather complex due to the several energy scales involved. In particular, as the Andreev bound states play a crucial role in the transport properties, a detailed analysis of electron correlation effects on these states seems desirable. In this paper we present a theoretical study of the Josephson effect in a single level quantum dot coupled to superconducting electrodes. Different approaches are considered. We first study the zero band-width limit in which the problem can be exactly diagonalized. Then, as an approximation to the full model we use a diagrammatic expansion of the self-energy associated with the Coulomb interaction in which the coupling to the superconducting leads is taken into account up to infinite order. This last ingredient is important for a proper description of the Andreev bound states.

In section II, we introduce the theoretical model and discuss the diagrammatic approximations for the electron self-energy. In section III we present a simplified version of this model in which the Coulomb interaction is replaced by an effective exchange field. This simple model already describes the $`\pi `$-junction transition and allows to understand the behavior of the Andreev states under a finite magnetization in the dot. Section IV is devoted to the analysis of the zero band-width limit in which the model becomes equivalent to a finite system which can be solved exactly. The exact solution is then used as a test of the approximations used for the electron self-energy. It turns out that the zero band-width limit can describe rather accurately most of the properties of the full model. The results for the full model are presented in section V. We first discuss two opposite limiting situations in which the coupling of the dot to the leads is either small or large compared to the superconducting gap. As in the zero band-width limit, in the case of weak coupling we show that the problem can be solved exactly. We also present some numerical results for the intermediate regime. The paper is closed by a brief summary of the main conclusions.

## II Model and theoretical approach

We describe a small quantum dot connected to superconducting electrodes by means of a modified single-level Anderson model, given by

$$\widehat{H}=\widehat{H}_L+\widehat{H}_R+\underset{\sigma }{}ϵ_0\widehat{n}_\sigma +U\widehat{n}_{}\widehat{n}_{}+\widehat{H}_T,$$
(1)
where $`\widehat{n}_\sigma =\widehat{c}_{0\sigma }^{}\widehat{c}_{0\sigma }`$, $`\widehat{H}_L`$ and $`\widehat{H}_R`$ represent the uncoupled superconducting leads; $`\widehat{H}_T=_{kL,R;\sigma }t_{0,k}\widehat{c}_{0\sigma }^{}\widehat{c}_{k,\sigma }+h.c.`$ describing the coupling between the dot level and the leads. Within this model the dot is represented by a single spin degenerate level with a repulsive Coulomb interaction described by the $`U`$-term in Eq. (1). We shall assume that the superconducting leads are well described by the BCS theory with a superconducting gap $`\mathrm{\Delta }`$ and that there is a fixed superconducting phase difference $`\varphi =\varphi _L\varphi _R`$ between both electrodes.

The relevant quantities like the current and the spectral densities can be expressed in terms of non-equilibrium Green functions . For the description of the superconducting state it is useful to introduce spinor field operators , which in a site representation are defined as:

$$\widehat{\psi _i}=\left(\begin{array}{c}c_i\\ c_i^{}\end{array}\right)\text{ , }\widehat{\psi }_i^{}=\left(\begin{array}{cc}c_i^{}& c_i\end{array}\right)$$
(2)
Then, the different correlation functions appearing in the Keldysh formalism adopt the standard causal form:

$$\widehat{G}_{ij}^{\alpha ,\beta }(t_\alpha ,t_\beta ^{})=i<\widehat{T}[\widehat{\psi }_i(t_\alpha )\widehat{\psi }_i^{}(t_\beta ^{})]>$$
(3)
where $`\widehat{T}`$ is the chronological ordering operator along the Keldysh time contour. The labels $`\alpha `$ and $`\beta `$ refer to the upper ($`\alpha +`$) and lower ($`\alpha `$) branches on this contour. The functions $`\widehat{G}_{ij}^+`$, which can be associated within this formalism with the electronic non-equilibrium distribution functions, are given by the (2x2) matrix:

$$\widehat{G}_{i,j}^+(t,t^{})=i\left(\begin{array}{cc}<c_j^{}(t^{})c_i(t)>& <c_j(t^{})c_i(t)>\\ <c_j^{}(t^{})c_i^{}(t)>& <c_j(t^{})c_i^{}(t)>\end{array}\right).$$
(4)
In terms of the Fourier transform matrix elements of $`\widehat{G}_{ij}^+(t,t^{})`$ one can write the charge and the induced order parameter on the dot, as well as the Josephson current:

$$<\widehat{n}_{}>=\frac{1}{2\pi i}_{\mathrm{}}^{\mathrm{}}𝑑\omega [\widehat{G}_{00}^+(\omega )]_{11},$$
(5)
$$<c_0^{}c_0^{}>=\frac{1}{2\pi i}_{\mathrm{}}^{\mathrm{}}𝑑\omega [\widehat{G}_{00}^+(\omega )]_{21},$$
(6)
$$I_{L,R}=\frac{e}{h}_{\mathrm{}}^{\mathrm{}}\underset{kL,R}{}d\omega \text{Tr}\left[\sigma _z\widehat{t}_{0,k}\widehat{G}_{k,0}^+(\omega )\sigma _z\widehat{t}_{k,0}\widehat{G}_{0,k}^+(\omega )\right],$$
(7)
where $`I_{L,(R)}`$ denotes the current between the left (right) lead and the dot, $`\sigma _z`$ is the usual Pauli matrix, and $`\widehat{t}_{0,k}`$ is the hopping matrix in the Nambu representation

$$\widehat{t}_{0,k}=\left(\widehat{t}_{k,0}\right)^{}=\left(\begin{array}{cc}t_{0,k}& 0\\ 0& t_{0,k}^{}\end{array}\right).$$
(8)
For the zero voltage case the calculation of the different $`\widehat{G}^+(\omega )`$ is particularly simple because the following relation holds:

$$\widehat{G}_{ij}^+(\omega )=f(\omega )[\widehat{G}_{ij}^a(\omega )\widehat{G}_{ij}^r(\omega )],$$
(9)
where $`f(\omega )`$ is the Fermi distribution function, and $`\widehat{G}^{a,r}`$ are the advanced and retarded Green functions. Therefore, the relevant quantity to be determined is the dot retarded Green function, which in a Nambu $`2\times 2`$ representation adopts the form

$$\widehat{G}_{00}^r(\omega )=\left[\omega \widehat{I}ϵ_0\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}^r(\omega )\widehat{\mathrm{\Gamma }}_L(\omega )\widehat{\mathrm{\Gamma }}_R(\omega )\right]^1,$$
(10)
where $`\widehat{\mathrm{\Gamma }}_{L,R}`$ are the tunneling rates given $`\widehat{\mathrm{\Gamma }}_k=\pi \rho _F\widehat{t}_{0,k}\widehat{g}^k\widehat{t}_{k,0}`$, with $`g_{11}^k=g_{22}^k=\omega /\sqrt{\mathrm{\Delta }^2\omega ^2}`$ and $`g_{12}^k=(g_{21}^k)^{}=\mathrm{\Delta }/\sqrt{\mathrm{\Delta }^2\omega ^2}e^{i\varphi _k}`$ with $`k=L,R`$ and $`\rho _F`$ is the normal density of states at the Fermi level (the chemical potential of the superconducting leads is taken as zero). The self-energy $`\widehat{\mathrm{\Sigma }}^r(\omega )`$ takes into account the effect of Coulomb interactions. To the lowest order in $`U`$ this is given by the Hartree-Fock Bogoliubov approximation: $`(\widehat{\mathrm{\Sigma }}^r)_{11}=U<\widehat{n}_{}>`$, $`(\widehat{\mathrm{\Sigma }}^r)_{22}=U<\widehat{n}_{}>`$; $`(\widehat{\mathrm{\Sigma }}^r)_{21}=(\widehat{\mathrm{\Sigma }}^r)_{12}^{}=U\widehat{c}_0^{}\widehat{c}_0^{}`$. We shall discuss in II A how correlation effects beyond this mean field approximation can be included.

This model has been analyzed within a mean field approximation in Ref. . In that work the mean field solution was simplified by neglecting the induced order parameter in the dot and imposing self-consistency only in the dot magnetization. However, the complete Hartree-Fock Bogoliubov solution requires the self-consistent determination of both the diagonal and non-diagonal charges in the dot. It should be noticed that the self-consistent determination of the induced complex order parameter is in principle necessary in order to ensure current conservation .

### A Inclusion of correlation effects

In order to go beyond the mean field approximation, we consider the diagrammatic expansion of the self-energy in terms of one-electron propagators in Nambu space. In Fig. 1 we show the corresponding second order diagrams in the electron-electron interaction. Due to the appearance of the anomalous propagators $`G_{12}`$ and $`G_{21}`$ in the superconducting state, there are additional diagrams to the one contributing in the normal state, labeled by 11(a) in Fig. 1, corresponding to the interaction of a quasi-particle with an electron hole pair with opposite spin.

The proper choice of the unperturbed one-electron Hamitonian over which the diagrammatic expansion is performed is an important issue. For the normal Anderson model in the symmetric case $`(ϵ_0=U/2)`$ the Hartree approximation which renormalizes the dot level as $`ϵ_0+U/2`$ is the adequate starting point for the perturbation theory as it automatically warrants charge consistency between the perturbed and the unperturbed situation. However, in a non-symmetric case perturbation over the Hartree field yields pathological results close to half-filling (see ). As discussed in Ref.  a better choice is to define an effective dot level $`ϵ_{eff}`$ in such a way that charge consistency between the effective and the interacting problems is achieved. The natural extension of this scheme to the superconducting case is to impose also consistency in the non-diagonal charge $`n_{12}`$ by introducing an effective local pairing potential $`\mathrm{\Delta }_{eff}`$ in the unperturbed Hamiltonian .

With these definitions, the dot Green functions in the unperturbed effective problem are given by (hereafter we omit the subscript $`00`$ in the dot Green functions)

$$\widehat{G}^{r(0)}(\omega )=\frac{1}{D(\omega )}\left(\begin{array}{cc}\omega ϵ_{eff}t_L^2g_{11}^Lt_R^2g_{11}^R& \mathrm{\Delta }_{eff}+t_L^2g_{12}^L+t_R^2g_{12}^R\\ \mathrm{\Delta }_{eff}^{}+t_L^2g_{21}^L+t_R^2g_{21}^R& \omega +ϵ_{eff}t_L^2g_{22}^Lt_R^2g_{22}^R\end{array}\right),$$
(11)
where $`D(\omega )`$ is the corresponding determinant. Notice that the Andreev states in the unperturbed problem are determined by the condition

$$D(\omega )=0,$$
(12)
which, in the general case, can have up to four solutions.

The different contributions to the self-energy represented in Fig. 1 can then be written in terms of the one-electron propagators in the following way

$`\mathrm{\Sigma }^{r(2)}(\omega )_{11,a}`$ $`=`$ $`{\displaystyle \frac{U^2}{(2\pi i)^3}}{\displaystyle 𝑑ϵ_1𝑑ϵ_2𝑑ϵ_3}`$ (14)
$`{\displaystyle \frac{G_{11}^{(0)+}(ϵ_1)G_{22}^{(0)+}(ϵ_2)G_{22}^{(0)+}(ϵ_3)+G_{11}^{(0)+}(ϵ_1)G_{22}^{(0)+}(ϵ_2)G_{22}^{(0)+}(ϵ_3)}{\omega ϵ_1ϵ_2+ϵ_3+i0^+}},`$
$`\mathrm{\Sigma }^{r(2)}(\omega )_{11,b}`$ $`=`$ $`{\displaystyle \frac{U^2}{(2\pi i)^3}}{\displaystyle 𝑑ϵ_1𝑑ϵ_2𝑑ϵ_3}`$ (16)
$`{\displaystyle \frac{G_{12}^{(0)+}(ϵ_1)G_{21}^{(0)+}(ϵ_2)G_{22}^{(0)+}(ϵ_3)+G_{12}^{(0)+}(ϵ_1)G_{21}^{(0)+}(ϵ_2)G_{22}^{(0)+}(ϵ_3)}{\omega ϵ_1ϵ_2+ϵ_3+i0^+}},`$
$`\mathrm{\Sigma }^{r(2)}(\omega )_{21,a}`$ $`=`$ $`{\displaystyle \frac{U^2}{(2\pi i)^3}}{\displaystyle 𝑑ϵ_1𝑑ϵ_2𝑑ϵ_3}`$ (18)
$`{\displaystyle \frac{G_{21}^{(0)+}(ϵ_1)G_{12}^{(0)+}(ϵ_2)G_{21}^{(0)+}(ϵ_3)+G_{21}^{(0)+}(ϵ_1)G_{12}^{(0)+}(ϵ_2)G_{21}^{(0)+}(ϵ_3)}{\omega ϵ_1ϵ_2+ϵ_3+i0^+}},`$
$`\mathrm{\Sigma }^{r(2)}(\omega )_{21,b}`$ $`=`$ $`{\displaystyle \frac{U^2}{(2\pi i)^3}}{\displaystyle 𝑑ϵ_1𝑑ϵ_2𝑑ϵ_3}`$ (20)
$`{\displaystyle \frac{G_{22}^{(0)+}(ϵ_1)G_{11}^{(0)+}(ϵ_2)G_{21}^{(0)+}(ϵ_3)+G_{22}^{(0)+}(ϵ_1)G_{11}^{(0)+}(ϵ_2)G_{21}^{(0)+}(ϵ_3)}{\omega ϵ_1ϵ_2+ϵ_3+i0^+}}.`$
Notice that $`G^{(0)+}/2\pi i`$ and $`G^{(0)+}/2\pi i`$ in the above equations correspond respectively to the occupied and unoccupied states in the dot spectral density.

The total self-energy is obtained by adding the contributions labeled by (a) and (b) in the above equations. The remaining self-energy components ($`\mathrm{\Sigma }_{22}^{r(2)}`$ and $`\mathrm{\Sigma }_{12}^{r(2)}`$) are obtained by similar expressions. In the non-magnetic case the diagonal self-energy components are related by $`\mathrm{\Sigma }_{22}^{r(2)}(\omega )=\mathrm{\Sigma }_{11}^{a(2)}(\omega )`$. This relation does not hold in the general magnetic case, except when $`ϵ_0=U/2`$ for which $`\mathrm{\Sigma }_{11}^{r(2)}(\omega )=\mathrm{\Sigma }_{22}^{r(2)}(\omega )`$.

In order to extend the range of validity of the second order self-energy an ansatz which interpolates correctly between the limits $`U/\mathrm{\Gamma }0`$ and $`U/\mathrm{\Gamma }\mathrm{}`$ can be used . Within this interpolative approach the matrix self-energy is given by the following expression

$$\widehat{\mathrm{\Sigma }}(\omega )=U\widehat{n}\widehat{\sigma }_z+\mathrm{\Delta }_d\widehat{\sigma }_x+\left[\widehat{I}\alpha \widehat{\mathrm{\Sigma }}^{(2)}\widehat{\sigma }_z\right]^1\widehat{\mathrm{\Sigma }}^{(2)}(\omega ),$$
(21)
where

$`\alpha ={\displaystyle \frac{ϵ_0+(1<\widehat{n}>)Uϵ_{eff}}{U^2<\widehat{n}>(1<\widehat{n}>)}},`$

and $`\widehat{\mathrm{\Sigma }^{(2)}}`$ is the second order self-energy.

This interpolative scheme has been used in many different contexts involving strongly correlated electrons . Notice, however, that for the electron-hole symmetric case the interpolative self-energy reduces to the second order one.

## III Transition from normal to $`\pi `$ junction behavior: a simple model

An exactly solvable model which can describe the transition from normal to $`\pi `$ junction behavior is obtained by substituting the interacting region by a single site with a local exchange field $`E_{ex}`$, in such a way that $`ϵ_{}=ϵ_0+E_{ex}`$ and $`ϵ_{}=ϵ_0E_{ex}`$. It should be noticed that this model is formally equivalent to a mean field solution of Hamiltonian (1) with the prescription $`E_{ex}=U(n_{}n_{})/2`$.

This model exhibits in general four bound states inside the superconducting gap (Andreev states) which give the dominant contribution to the current. This is determined by the derivative of the states below the Fermi energy with respect to the phase. The position and spectral weight of these states can be obtained from the retarded Green function (Eq. 10). In the limit $`\mathrm{\Delta }<<\mathrm{\Gamma }`$ and when $`\mathrm{\Gamma }_L=\mathrm{\Gamma }_R=\mathrm{\Gamma }`$, the Andreev states can be determined analytically by the expression

$`\left({\displaystyle \frac{\omega _\pm }{\mathrm{\Delta }}}\right)^2={\displaystyle \frac{\mathrm{cos}^2\varphi /2+2E^2+Z^2(Z^2+\mathrm{sin}^2\varphi /2)\pm 2XS(\varphi )}{Z^4+2(X^2+E^2)+1}},`$ (22)
where $`E=ϵ_0/2\mathrm{\Gamma }`$, $`X=E_{ex}/2\mathrm{\Gamma }`$ and $`Z^2=X^2E^2`$ and $`S(\varphi )`$ is given by

$`S(\varphi )=\sqrt{Z^2\mathrm{cos}^2\varphi /2+E^2+\mathrm{sin}^2\varphi /4}.`$

It can be verified that in the limit $`X0`$ Eq. (22) recovers the well known expression for the Andreev states in a single channel contact of transmission $`\tau =1/(1+E^2)`$, i.e. $`\omega /\mathrm{\Delta }=\pm \sqrt{1\tau \mathrm{sin}^2\varphi /2}`$ .

The qualitative behavior of the Andreev states can be easily analyzed in the particular case $`E=0`$, in which the position of the four states is given by

$$\omega _\pm ^\pm =\pm \mathrm{\Delta }\frac{\mathrm{cos}\varphi /2\pm X\sqrt{\mathrm{sin}^2\varphi /2+X^2}}{(X^2+1)}.$$
(23)
This expression clearly shows that the effect of a finite exchange field is to split the non-magnetic Andreev states at $`\pm \mathrm{\Delta }\mathrm{cos}\varphi /2`$ into four states. The transition to the $`\pi `$ junction behavior is associated with the progressive inter-crossing of the two ”inner” bound states. This crossing is illustrated in Fig. 2 for a generic situation. For small exchange field (Fig. 2 a), the inner states do not cross and the system is in the $`0`$ state. As $`X`$ increases, the inner states cross at $`\pi `$ and thus the current-phase relation exhibits a sign change around $`\varphi =\pi `$. Eventually, when the crossing between the inner states is complete the whole current-phase relation changes sign and the system is in the $`\pi `$ state. The state in the intermediate region is conventionally designed as $`0^{}`$ or $`\pi ^{}`$ depending on the relative stability of the minima of the inter-crossing states.

The boundaries between the different regions are straightforwardly obtained from Eq. (22): the curves $`X=E`$, $`X=(E+\sqrt{3+4E^2})/3`$ and $`X=\sqrt{1+E^2}`$ correspond to the $`00^{}`$, $`0^{}\pi ^{}`$ and $`\pi ^{}\pi `$ boundaries respectively. The full phase diagram of this simple model is illustrated in Fig. 3.

## IV Zero band-width limit

A special limit in which the model given by Eq. (1) becomes exactly solvable is the case where the band-width of the electrodes tends to zero. This is equivalent to neglect the high energy excitations in the superconducting electrodes, substituting the frequency dependent electrodes self-energy by an effective diagonal and non-diagonal level in Nambu space . In this limit the semi-infinite leads connected to the dot are replaced by an effective single site and the terms $`H_{L,R}`$ describing the leads in Eq. (1) are given by

$$H_{L,R}=\underset{\sigma }{}ϵ_{L,R}\widehat{c}_{L,R\sigma }^{}\widehat{c}_{L,R\sigma }+\mathrm{\Delta }_{L,R}\widehat{c}_{L,R}\widehat{c}_{L,R}+\mathrm{\Delta }_{L,R}^{}\widehat{c}_{L,R}^{}\widehat{c}_{L,R}^{}.$$
(24)
The system thus becomes a sort of “superconducting molecule” with a finite number of electron configurations which can be diagonalized exactly. This solution can be useful as a test of the accuracy of the proposed approximations for the infinite system. In spite of its simplicity this limit also gives a qualitative description of the behavior of the full model, as will be shown below.

The exact solution in this limit is obtained by considering electronic configurations with all possible total number of electrons, the pairing interaction connecting configurations that differ in two electrons. The eigenstates can be classified into those arising from even or odd number configurations. In the even case the ground state has zero total spin, while in the odd case the ground state is two-fold degenerate corresponding to $`S_z=\pm 1/2`$.

In Fig. 4 we show the evolution of the ground state energy for the even and odd cases for fixed $`ϵ_0`$ and increasing $`U`$. As can be observed, in the even case the ground state energy exhibits always a minimum for $`\varphi =0`$, while in the odd case the minimum appears for $`\varphi =\pi `$. This is to be expected from the fact that the system magnetization is finite in the odd case. For small values of $`U`$ the even case is more stable for all values of the phase (Fig. 4 a), while the opposite situation is found for large enough $`U`$ (Fig. 4 d). In the intermediate region (Fig. 4 b and c) the even and odd energy levels cross. This level crossing corresponds to the transition between 0 and $`\pi `$ junction behavior in the full model. The $`\pi `$ state thus corresponds to a situation in which the dot adquires a finite magnetization.

The exact ground state energy of this simple model can now be used to check the accuracy of the different approximations discussed in Sect. II. For the sake of simplicity we shall restrict this comparison to the $`ϵ_0=U/2`$ case. In Fig. 5 we show the comparison between the exact, the mean field, and the second-order self-energy results for the ground state energy for different values of $`U/\mathrm{\Delta }`$. At small values of $`U/\mathrm{\Delta }`$ (Fig. 5a), the energy of the mean field approximation lies slightly above the exact result both in the non-magnetic and in the magnetic case. When increasing $`U/\mathrm{\Delta }`$ (Fig. 5b), the non-magnetic mean field solution increasingly deviates from the exact result. In contrast, in the magnetic case, the deviation between the exact and the mean field solution first increases in the small $`U`$ range while it progressively decreases for large $`U`$. This is due to the fact that the exact ground state corresponds to a localized spin at the central site for $`U\mathrm{}`$, which is correctly described by the magnetic mean field solution.

The inclusion of the second-order self-energy substantially improves the results of the mean field approximation both in the magnetic and in the non-magnetic cases for small and moderate values of $`U`$. In the non-magnetic case the improvement is considerable in the whole range of values of $`U`$ (Fig. 5a and 5b). On the other hand, in the magnetic case the improvement due to the inclusion of correlation effects progressively becomes less important as the mean field tends to the exact solution in the large $`U`$ limit.

Finally, it is interesting to analyze the phase-diagram of the model in the zero band-width limit, which turns out to contain the essential features of the full model with the four types of solutions: 0, $`0^{}`$, $`\pi ^{}`$ and $`\pi `$ appearing in the previous section. In this discrete model the different phases can be identified by the phase-dependence of the ground state energy.

The exact phase diagram for the zero band-width limit is shown in Fig. 6a. Only the region $`U>0`$ and $`ϵ<0`$ is shown as for all other regions the only possible phase is $`0`$ type. Roughly, the $`\pi `$ state is found for $`U/\mathrm{\Gamma }`$ and $`ϵ_0/\mathrm{\Gamma }`$ sufficiently large, while the $`0`$ state appears either for $`U/\mathrm{\Gamma }`$ or $`ϵ_0/\mathrm{\Gamma }`$ sufficiently small. The $`0^{}`$ and the $`\pi ^{}`$ behavior is found in the intermediate regime. It is interesting to notice that for $`ϵ_0/\mathrm{\Gamma }<1`$ and $`U/\mathrm{\Gamma }\mathrm{}`$ the $`0`$ state is always the more stable. An extrapolation of this result to the full model would imply the absence of $`\pi `$ junction behavior in the so called mixed valence regime, as has been predicted in Ref. .

For comparison, the mean field diagram is also shown in Fig. 5b. Notice that although the overall behavior is captured in the mean field solution, correlation effects displace the $`\pi `$ region to much larger values of $`U/\mathrm{\Gamma }`$ and $`ϵ_0/\mathrm{\Gamma }`$. Neglecting the induced order parameter in the dot as done in Ref.  does not qualitatively change the phase diagram although the full self-consistent solution is somewhat closer to the exact result.

## V Results for the full model

The mean field approximation for the full model has been analyzed in Ref. . Although in this work self-consistency in the induced order parameter was neglected, its effects on the total energy is small as already discussed for the zero band-width limit. In the present section we will concentrate in discussing the effects of electronic correlations beyond the mean field solution. We will consider the electron-hole symmetric situation ($`ϵ_0=U/2`$) in which electron correlation effects are expected to be more important.

For a fixed value of $`U`$, the physical behavior of the model is controlled by the dimensionless parameter $`\mathrm{\Delta }/\mathrm{\Gamma }`$. Is is interesting to analyze in detail the two opposite limits $`\mathrm{\Delta }/\mathrm{\Gamma }1`$ and $`\mathrm{\Delta }/\mathrm{\Gamma }1`$

$`\mathrm{\Delta }/\mathrm{\Gamma }1`$ limit: This situation corresponds to a dot very weakly coupled to the leads. In this limit the exact solution can be obtained due to the fact that the problem can be mapped into a two-level Hamiltonian describing the dynamics of the Andreev states. This can be shown by first considering the non-interacting ($`U=0`$) situation. In this case, the spectral density of the dot exhibits bound states at energies $`\omega _s=\pm 2\mathrm{\Gamma }\mathrm{cos}\varphi /2`$ and the continuous part becomes negligible. Consequently, the retarded Green function at the dot can be simply written as

$$\widehat{G}^r(\omega )\left(\begin{array}{cc}\omega & \omega _s\\ \omega _s& \omega \end{array}\right)^1,$$
(25)
indicating that the corresponding effective Hamiltonian is given by $`\widehat{H}^{eff}=\omega _s[\widehat{c}_{}\widehat{c}_{}+\widehat{c}_{}^{}\widehat{c}_{}^{}]`$. When introducing a finite $`U`$, the diagonal level at $`U/2`$ is canceled by the Hartree term and the remaining part of the Hamiltonian associated with the Coulomb interaction becomes $`U(n_{}1/2)(n_{}1/2)U/4`$. This term vanishes for the case of an empty and a doubly occupied dot, while it yields an energy $`U/2`$ for the single electron case. As a consequence, the ground state of the system is either the symmetric combination of empty and doubly occupied configurations (with energy $`\omega _s`$) or the doubly degenerate single occupied state (with energy $`U/2`$). Thus, the transition to the magnetic $`\pi `$ state takes place when $`U/2>2\mathrm{\Gamma }`$.

$`\mathrm{\Delta }/\mathrm{\Gamma }1`$ limit: In this opposite limit one would expect to recover gradually the properties of a normal system. In particular, for $`U>\mathrm{\Gamma }`$ the features associated with the Kondo effect should emerge. Fig. 7 shows the dot spectral density in this regime for increasing values of $`U`$ obtained in the second-order self-energy approximation. For comparison the spectral density in the normal state is also shown. As can be observed, the spectral density is similar to the one found in the normal state except for the superimposed features in the gap region. The overall shape evolves as in the normal case from a single Lorentzian broad resonance for $`U<\mathrm{\Gamma }`$ to the three peaked structure characteristic of the Kondo regime when $`U>\mathrm{\Gamma }`$. In this regime the relevant energy scale is set by the Kondo temperature $`T_K`$ which essentially measures the width of the Kondo resonance in the normal state. Within the second-order self-energy approximation $`T_K\mathrm{\Gamma }/(1\alpha _0)`$, where $`\alpha _0=\frac{\mathrm{\Sigma }}{\omega }(0)`$. In the symmetric case

$$\alpha _0=\left(\frac{U}{2\pi \mathrm{\Gamma }}\right)^2\left(3\frac{\pi ^2}{4}\right),$$
(26)
which coincides with the perturbative result of Ref.  for the Anderson model. Although this approach fails to give the exponential decay of the Kondo temperature with $`U`$ it gives a rather good description of the spectral density in the moderate $`U/\mathrm{\Gamma }`$ range .

The coexistence of the Kondo and Josephson effects is to be expected as far as $`T_K>\mathrm{\Delta }`$. When $`U`$ is further increased the system should evolved into the magnetic $`\pi `$ state with the suppression of the Kondo effect.

Let us analyze the self-energy in this limit. The effective one-electron problem in this case is characterized by the presence of Andreev states which, from Eq. (11), are located at $`\omega _s=\mathrm{\Delta }\sqrt{1\tau \mathrm{sin}^2\varphi /2}`$, where $`\tau =4\mathrm{\Gamma }^2/(ϵ_{eff}^2+4\mathrm{\Gamma }^2)`$ is the normal transmission in the effective problem (in the electron-hole symmetric case considered here $`ϵ_{eff}=0`$ and $`\tau =1`$), like in the case of a single channel point-contact. The weight of the Andreev states at the dot site decreases as $`\mathrm{\Delta }/\mathrm{\Gamma }`$ according to the expression $`\mathrm{\Delta }|\mathrm{sin}\varphi /2|/4\mathrm{\Gamma }`$. Also the induced order parameter tends to zero as $`\mathrm{\Delta }/\mathrm{\Gamma }`$ in this limit.

The electron self-energy can then be evaluated retaining only terms of order $`\mathrm{\Delta }/\mathrm{\Gamma }`$. The contributions labeled as 11(b) and 21(a) can be altogether neglected as they involve more than one anomalous propagator and are thus of order $`(\mathrm{\Delta }/\mathrm{\Gamma })^2`$. On the other hand, from the general expressions for the second-order self-energy, Eqs. (14,20), it should be noticed that there are three types of contributions to $`\mathrm{\Sigma }^{(2)}`$: one involving only the discrete part of the one-electron spectral densities, another one involving only the continuous part; and terms in which both the localized states and the continuous spectrum are mixed. The first of these contributions is of order $`(\mathrm{\Delta }/\mathrm{\Gamma })^3`$ and can be neglected. The resulting expression for $`\mathrm{\Sigma }^{2)}`$ up to first order in $`\mathrm{\Delta }/\mathrm{\Gamma }`$ for $`|\omega |<\mathrm{\Delta }`$ is

$`\mathrm{\Sigma }_{11}^{(2)}(\omega )`$ $``$ $`\left({\displaystyle \frac{U}{2\pi \mathrm{\Gamma }}}\right)^2\left(3{\displaystyle \frac{\pi ^2}{4}}\right)\omega `$ (27)
$`\mathrm{\Sigma }_{12}^{(2)}(\omega )`$ $``$ $`\left({\displaystyle \frac{U}{2\pi \mathrm{\Gamma }}}\right)^2\left[|\mathrm{sin}(\varphi /2)|+2(1+{\displaystyle \frac{1}{\pi }})\mathrm{cos}(\varphi /2)\right]\mathrm{\Delta }.`$ (28)
These expressions for the self-energies allow to determine the renormalization of the states inside the gap due to the interactions. For moderate values of $`U`$ ($`U<10\mathrm{\Gamma }`$) the renormalized states have approximately the same phase-dependence as in the non-interacting case (i.e. $`\mathrm{cos}(\varphi /2)`$) but with a narrower dispersion given by $`\stackrel{~}{\omega }_s(0)\mathrm{\Delta }(1(U/U_0)^2)`$, where $`(U_0/\mathrm{\Gamma })^2=(\mathrm{\Gamma }/\mathrm{\Delta })\pi ^2/(2\pi +2)`$. For larger values of $`U`$ the phase-dependence of the states starts to deviate from this simple law. This is illustrated in Fig. 8.

The evolution of the renormalized Andreev states indicate that the critical currents are suppressed as $`(1(U/U_0)^2)`$ in this limit (for moderate values of $`U`$).

One can summarize the results for the full model in the symmetric case by discussing the phase diagram shown in Fig. 9. In this figure we compare the results of the mean field and the second-order self-energy approximations for the critical $`U`$, $`U_c`$, defining the transition to the $`\pi `$ state. As can be observed both approximations yield the same result in the $`\mathrm{\Gamma }/\mathrm{\Delta }1`$ limit. A close inspection of this regime, illustrated by the inset in Fig. 9, shows that the transition takes place for $`U_c4\mathrm{\Gamma }`$ when $`\mathrm{\Gamma }0`$ in agreement with the exact result. As $`\mathrm{\Gamma }/\mathrm{\Delta }`$ increases the predictions of the two approaches start to deviate. The mean field predicts an almost linear relation between $`U_c`$ and $`\mathrm{\Gamma }`$ although with a larger slope than in the $`\mathrm{\Gamma }\mathrm{\Delta }`$ limit. On the other hand, the second order self-energy predicts a faster increase of $`U_c`$ with $`\mathrm{\Gamma }`$. Notice that in the normal state ($`\mathrm{\Delta }0`$) there should be no transition into a magnetic state for finite $`U`$. We also show as a reference the line where $`T_K`$ in the normal state, given by the second-order self-energy approximation, matches $`\mathrm{\Delta }`$. It can be observed that the criterium $`T_K\mathrm{\Delta }`$ for the transition to the $`\pi `$ state is a rather good for $`\mathrm{\Gamma }`$ sufficiently large.

## VI Conclusions and final remarks

We have presented a theoretical analysis of the Josephson transport through a strongly correlated level coupled to superconducting electrodes. The analysis has been specially focused in the effect of electron correlations on the subgap Andreev states. These states determine to a large extent the current-phase relation in this system. The transition to a $`\pi `$ state can be understood as a result of the intercrossing of the subgap states induced by an increasing Coulomb interaction. Within this model this transition corresponds to a truly quantum phase transition in which the ground state becomes degenerate having a localized magnetic moment at the dot level as already noticed in Ref. . It is worth noticing that this situation cannot exist in the absence of superconductivity. In fact, this behavior can be traced to the suppression of the Kondo effect (in which the electrons at the dot level couple to a singlet) due to the absence of low energy excitations in the superconducting leads.

In the present analysis we have used different approximation methods. In order to get insight on the behavior of the Andreev states in the transition to the $`\pi `$ state we have first analyzed a simple mean field model in which the electron interactions are represented by a local exchange field. We have also studied the zero band-width limit which allows for an exact diagonalization. It has been shown that the study of this limit already illustrates the different types of behaviors that can be found in the full model. These results have also been used to determine the accuracy of the self-energy approach, showing that it improves considerably the results of the mean field approximation for moderate values of the Coulomb interaction. Finally, we have presented results for the full model in different regimes. In the limit $`\mathrm{\Delta }\mathrm{\Gamma }`$ we have shown that the problem can be solved exactly, its dynamics being described by a two level Hamiltonian corresponding to the Andreev states. On the other hand, for $`\mathrm{\Gamma }\mathrm{\Delta }`$ there is a coexistence of Kondo and Josephson effects for $`T_K>\mathrm{\Delta }`$, the main effect of electron correlations being included in a renormalization of the critical current.

The present work constitutes a first step in the study of the transport properties of a quantum dot coupled to superconducting leads in a non-equilibrium situation, i.e. with an applied bias voltage. This would allow to analyze an experimental situation like the one of Ref. . Work along this line is under progress.


text:
# 𝐻₂⁺ molecular ion in a strong magnetic field: ground state

## I Introduction

Many years have passed since the moment when theoretical qualitative arguments were given that show that in the presence of a strong magnetic field the physics of atoms and molecules exhibits a wealth of new, unexpected phenomena even for the simplest systems Kadomtsev:1971 ; Ruderman:1971 . In particular, a chance that unusual chemical compounds may be formed which do not exist without magnetic field was mentioned. In practice, the atmosphere of neutron stars, which is characterized by the presence of enormous magnetic fields $`10^{12}10^{13}G`$, as well as other astronomical objects carrying large magnetic fields ($`>10^8G`$) provide a valuable paradigm where this physics could be realized. Recently, the experimental data collected by the Chandra X-ray observatory revealed certain irregularities in the spectrum of an isolated neutron star 1E1207.4-5209. These irregularities can be interpreted as absorption features at $`0.7`$ KeV and $`1.4`$ KeV of possible atomic or molecular nature Sandal:2002 .

One of the first general features observed in standard atomic and molecular systems placed in a strong magnetic field is an increase of both total and binding energies, accompanied by a drastic shrinking of the electron localization length in both the longitudinal and transverse directions. Naturally, this leads to a decrease of the equilibrium distance with magnetic field growth. This behavior can be considered to be a consequence of the fact that for large magnetic fields the electron cloud takes a needle-like form extended along the magnetic field direction and the system becomes effectively quasi-one-dimensional Ruderman:1971 . It is obvious that the phenomenon of quasi-one-dimensionality enhances the stability of standard atomic and molecular systems from the electrostatic point of view. In particular, molecules become elongated along the magnetic line forming a type of linear molecular polymer (for details see the review articles Garstang:1977 ; Lai:2001 ). It also hints at the occurrence of exotic atomic and molecular systems which do not exist in the absence of a magnetic field. Motivated by these simple observations it was shown in Refs.Turbiner:1999 ; Lopez-Tur:2000  that three and even four protons can be bound by one electron. This shows that exotic one-electron molecular systems $`H_3^{2+}`$ and $`H_4^{3+}`$ can exist in sufficiently strong magnetic fields in the form of linear polymers. However, the situation becomes much less clear (and also much less studied) when the nuclei are not aligned with the magnetic field direction, and thus in general do not form a linear system. Obviously, such a study would be important for understanding the kinetics of a gas of molecules in the presence of a strong magnetic field. As a first step towards such a study, even the simplest molecules in different spatial configurations deserve attention. Recently, a certain spatial configuration of $`H_3^{2+}`$ was studied in detail Lopez-Tur:2002 . It was shown that in the range of magnetic fields $`10^8<B<10^{11}G`$ the system $`(pppe)`$, with the protons forming an equilateral triangle perpendicular to the magnetic lines, has a well-pronounced minimum in the total energy for a certain size of triangle. The goal of the present work is to attempt for the first time to carry out an extensive quantitative investigation of the ground state of $`H_2^+`$ in the framework of a single approach in its entire complexity: a wide range of magnetic field strengths ($`04.414\times 10^{13}G`$), arbitrary (but fixed) orientation of the molecular axis with respect to the magnetic line and arbitrary internuclear distances. We are going to carry out this study in the Born-Oppenheimer approximation at zero order – assuming protons to be infinitely heavy charged centers. In principle, when the molecular axis is perpendicular to the magnetic line the system $`(ppe)`$ acquires extra stability from the electrostatic point of view. Electrostatic repulsion of the classical protons is compensated for by the Lorentz force acting on them.

It is well known that the molecular ion $`H_2^+`$ is the most stable one-electron molecular system in the absence of a magnetic field. It remains so in the presence of a constant magnetic field unless $`B10^{13}G`$, where the exotic ion $`H_3^{2+}`$ appears to be the most bound (see Lopez-Tur:2000 ). The ion $`H_2^+`$ has been widely studied, both with and without the presence of a magnetic field, due to its importance in astrophysics, atomic and molecular physics, solid state and plasma physics (see Garstang:1977 -Potekhin  and references therein). The majority of the previous studies were focused on the case of the parallel configuration, where the angle between the molecular axis and the magnetic field direction is zero, $`\theta =0^o`$. The only exception is Ref.Schmelcher , where a detailed quantitative analysis was performed for any $`\theta `$ but for a single magnetic field $`B=1a.u.`$. Previous studies were based on various numerical techniques, but the overwhelming majority used different versions of the variational method, including the Thomas-Fermi approach. As a rule, in these studies the nuclear motion was separated from the electronic motion using the Born-Oppenheimer approximation at zero order (see above). It was observed at the quantitative level that magnetic field growth is always accompanied by an increase in the total and binding energies, as well as a shrinking of the equilibrium distance. As a consequence it led to a striking conclusion about the drastic increase in the probability of nuclear fusion for $`H_2^+`$ in the presence of a strong magnetic field Kher .

In the present study we will also use a variational method. Our consideration will be limited by a study of the $`1_g`$-state, which realizes the ground state of the system if the bound state exists 111After a straightforward separation of the spin part of wavefunction, the original Schroedinger equation becomes a scalar Schroedinger equation. Then it can be stated that a nodeless eigenfunction corresponds to the ground state (Perron theorem). We will construct state-of-the-art, non-straightforward, ‘adequate’ trial functions consistent with a variationally optimized choice of vector potential. We should stress that a proper choice of the form of the vector potential is one of the crucial points which guarantee the adequacy and reliability of the consideration. In particular, a proper position of the gauge origin, where the vector potential vanishes, is drastically important, especially for large interproton distances. For the parallel configuration, $`\theta =0^o`$ the present work can be considered as an extension (and also an development) of our previous work Lopez:1997 . It is necessary to emphasize that we encounter several new physical phenomena which occur when the molecular axis deviates from the magnetic field direction. If the magnetic field is sufficiently strong, $`B10^{11}G`$ and inclination $`\theta `$ is larger than a certain critical angle the ion $`H_2^+`$ does not exist in the contrary to a prediction in Refs.Larsen , Kher , Turbiner:2001 . This prediction was based on an improper gauge dependence of the trial functions, which caused a significant loss of accuracy and finally led to a qualitatively incorrect result. We find that in the weak field regime the $`(ppe)`$ system in the equilibrium position at any inclination, the electronic distribution peaks at the positions of the protons. While at large magnetic fields the electronic distribution is characterized by single peak at the midpoint between two protons. This change from a two-peak to a one-peak configuration appears around $`B10^{10}10^{11}G`$ with a slight dependence on the inclination angle $`\theta `$. From a physical point of view the former means that the electron prefers to stay in the vicinity of a proton. This can be interpreted as a dominance of the $`H`$-atom plus proton configuration. The latter situation implies that the electron is ‘shared’ by both protons and hence such a separation to $`H`$-atom plus proton is irrelevant. Therefore, we can call the two-peak situation ‘ionic’ coupling, while the one-peak case can be designated as ‘covalent’ coupling, although this definition differs from that widely accepted in textbooks (see, for example LL ). Thus, we can conclude that a new phenomenon appears - as the magnetic field grows the type of coupling changes from ‘ionic’ to ‘covalent’. At large internuclear distances the electron is always attached to one of the charged centers, so the coupling is ‘ionic’.

One particular goal of our study is to investigate a process of dissociation of the $`(ppe)`$ system: $`H_2^+H+p`$ which appears with increase of interproton distance. It is clear from a physical point of view that at large distances the electronic distribution should be first of the two-peak type and then should change at asymptotically large distances, to a single-peak one, but with a peak at the position of one of the protons. Somehow this process breaks permutation symmetry and we are not aware of any attempt to describe it. In our analysis this phenomenon appears as a consequence of a change of position of the gauge origin with increase of interproton distance.

From the physical point of view it is quite interesting to note how the $`(ppe)`$ system behaves at very large interproton distances. This domain is modelled by an $`H`$-atom plus proton interaction. The interaction corresponds to (magnetic-field-inspired-quadrupole) + charge interaction and is dominant comparing to the standard Van der Waals force. For small inclinations the above interaction is attractive as in the Van der Waals case, but becomes repulsive for large inclinations. This implies that the potential curves approach to asymptotic value of total energy at the large interproton distances from above contrary to the Van der Waals case.

The Hamiltonian which describes two infinitely heavy protons and one electron placed in a uniform constant magnetic field directed along the $`z`$axis, $`𝐁=(0,0,B)`$ is given by (see e.g. LL )

$$=\widehat{p}^2+\frac{2}{R}\frac{2}{r_1}\frac{2}{r_2}+(\widehat{p}𝒜+𝒜\widehat{p})+𝒜^2,$$
(1)
(see Fig.1 for notations), where $`\widehat{p}=i`$ is the momentum, $`𝒜`$ is a vector potential, which corresponds to the magnetic field $`𝐁`$. Hence, the total energy $`E_T`$ of $`H_2^+`$ is defined as the total electronic energy plus the Coulomb energy of proton repulsion. The binding energy is defined as an affinity to have the electron at infinity, $`E_b=BE_T`$. The dissociation energy is defined as affinity to have a proton at infinity, $`E_d=E_HE_T`$, where $`E_H`$ is the total energy of the hydrogen atom in magnetic field $`B`$.

Atomic units are used throughout ($`\mathrm{}`$=$`m_e`$=$`e`$=1) albeit energies are expressed in Rydbergs (Ry). Sometimes, the magnetic field $`B`$ is given in a.u. with $`B_0=2.35\times 10^9G`$ 222In absence of convention, some results presented in literature are obtained for $`B_0=2.3505\times 10^9G`$. Thus, making a comparison of the results obtained by different authors, this fact should be taken into account.

## II Optimization of vector potential

It is well known that the vector potential for a given magnetic field, even in the Coulomb gauge $`(𝒜)=0`$, is defined ambiguously, up to a gradient of an arbitrary function. This gives rise a feature of gauge invariance: the Hermitian Hamiltonian is gauge-covariant, while the eigenenergies and other observables are gauge-independent. However, since we are going to use an approximate method for solving the Schroedinger equation with the Hamiltonian (1), our approximation of eigenenergies can well be gauge-dependent (only the exact ones are gauge-independent). Hence one can choose the form of the vector potential in a certain optimal way. In particular, if the variational method is used, the vector potential can be considered as a variational function and be chosen by a procedure of minimization.

Let us consider a certain one-parameter family of vector potentials corresponding to a constant magnetic field $`𝐁=(0,0,B)`$

$$𝒜=B((\xi 1)y,\xi x,0),$$
(2)
where $`\xi `$ is parameter, in the Coulomb gauge. The position of the gauge center or gauge origin, where $`𝒜(x,y,z)=0`$, is defined by $`x=y=0`$, with $`z`$ arbitrary. For simplicity we fix $`z=0`$. If $`\xi =1/2`$ we get the well-known and widely used gauge which is called symmetric or circular. If $`\xi =0`$ or 1, we get the asymmetric or Landau gauge (see LL ). By substituting (2) into (1) we arrive at a Hamiltonian of the form

$`=^2+{\displaystyle \frac{2}{R}}{\displaystyle \frac{2}{r_1}}{\displaystyle \frac{2}{r_2}}2iB[(\xi 1)y_x+\xi x_y]+B^2[\xi ^2x^2+(1\xi )^2y^2],`$ (3)
where $`R`$ is the interproton distance (see Fig. 1).

It is evident that for small interproton distances, $`R`$, the electron prefers to be near the mid-point between the two protons (coherent interaction with the protons). In the opposite limit, $`R`$ large, the electron is situated near one of the protons (this is an incoherent situation - the electron selects and then interacts essentially with one proton). This fact, together with naive symmetry arguments, leads us to a natural assumption that the gauge center is situated on a line connecting the protons. Therefore, the coordinates of mid-point between protons are

$$Y=\frac{Rd}{2}\mathrm{sin}\theta ,Z=\frac{Rd}{2}\mathrm{cos}\theta ,$$
(4)
(see Fig.1), where $`d`$ is a parameter. Thus, the position of the gauge center is effectively measured by the parameter $`d`$ – a relative distance between the middle of the line connecting the protons and the gauge center. If the mid-point coincides with the gauge center, then $`d=0`$. On other hand, if the position of a proton coincides with the gauge center, then $`d=1`$ or $`d=1`$. Hence the parameter $`d`$ makes sense as a parameter characterizing a gauge.

The idea of choosing an optimal (convenient) gauge has been widely exploited in quantum field theory calculations. It has also been discussed in quantum mechanics and, in particular, in connection to the present problem. Perhaps, the first constructive (and remarkable) attempt to realize the idea of an optimal gauge was made in the eighties by Larsen Larsen . In his variational study of the ground state of the $`H_2^+`$ molecular ion it was explicitly shown that for a given fixed trial function the gauge dependence of the energy can be quite significant. Furthermore, even an oversimplified optimization procedure improves the accuracy of the numerical results 333For review of different optimization procedures of vector potential see, for instance, Kappes:1994  and references therein.

Our present aim is to study the ground state of (1) or, more concretely, (3). We propose a different way of optimization of vector potential than those discussed by previous authors. It can be easily demonstrated that for a one-electron system there always exists a certain gauge for which the ground state eigenfunction is a real function. Let us fix a vector potential in (1). Assume that we have solved the spectral problem exactly and have found the ground state eigenfunction. In general, it is a certain complex function with a non-trivial, coordinate-dependent phase. Treating this phase as a gauge phase and then gauging it away, finally results in a new vector potential. This vector potential has the property we want – the ground state eigenfunction of the Hamiltonian (1) is real. It is obvious that similar considerations are valid for any excited state. In general, for a given eigenstate there exists a certain gauge in which the eigenfunction is real. For different eigenstates these gauges can be different. It is obvious that similar situation occurs for any one-electron system in a magnetic field.

Dealing with real trial functions has an immediate advantage: the expectation value of the terms proportional to $`𝒜`$ in (1) (or $`B`$ in (3)) vanishes when it is taken over any real, normalizable function. Thus, without loss of generality, the term $`B`$ in (3) can be omitted. Thus, we can use real trial functions with explicit dependence on the gauge parameters $`\xi `$ and $`d`$. These parameters are fixed by performing a variational optimization of the energy. Therefore, as a result of the minimization we find both a variational energy and a gauge for which the ground state eigenfunction is real, as well as the corresponding Hamiltonian. One can easily show that for a system possessing axial (rotational) symmetry 444 This is the case whenever the magnetic field is directed along the molecular axis (parallel configuration) the optimal gauge is the symmetric gauge $`\xi =1/2`$ with arbitrary $`d`$. This is precisely the gauge which has been overwhelmingly used (without any explanations) in the majority of the previous research on $`H_2^+`$ in the parallel configuration Kadomtsev:1971  \- Potekhin . However, this is not the case if $`\theta 0^{}`$. For the symmetric gauge the exact eigenfunction now becomes complex, therefore complex trial functions must be used. But following the recipe proposed above we can avoid complex trial functions by adjusting the gauge in such a way the eigenfunction remains real. This justifies the use of real trial functions. Our results (see Section IV) lead to the conclusion that for the ground state the optimal gauge parameter varies in the interval $`\xi [0.5,\mathrm{\hspace{0.17em}1}]`$.

## III Choosing trial functions

The choice of trial functions contains two important ingredients: (i) a search for the gauge leading to the real, exact ground state eigenfunction and (ii) performance of a variational calculation based on real trial functions. The main assumption is that a gauge corresponding to a real, exact ground state eigenfunction is of the type (2) (or somehow is close to it) 555It can be formulated as a problem – for a fixed value of $`B`$ and a given inclination to find a gauge for which the ground state eigenfunction is real.. In other words, one can say that we are looking for a gauge of type (2) which admits the best possible approximation of the ground state eigenfunction by real functions. Finally, in regard to our problem, the following recipe of variational study is used: *As the first step, we construct an adequate variational real trial function $`\mathrm{\Psi }_0`$ Tur , for which the potential $`V_0=\frac{\mathrm{\Delta }\mathrm{\Psi }_0}{\mathrm{\Psi }_0}`$ reproduces the original potential near Coulomb singularities and at large distances, where $`\xi `$ and $`d`$ would appear as parameters. The trial function should support symmetries of an original problem. We then perform a minimization of the energy functional by treating the free parameters of the trial function and $`\xi ,d`$ on the same footing.* In particular, such an approach enables us to find the *optimal* form of the Hamiltonian as a function of $`\xi ,d`$.

The Hamiltonian (1) gives rise to different symmetry properties depending on the orientation of the magnetic field with respect to the internuclear axis. The most symmetric situation corresponds to $`\theta =0^{}`$, where invariance under permutation of the (identical) charged centers $`P:(12)`$ together with $`P_z:(zz)`$ holds. Since the angular momentum projection $`\mathrm{}_z=m`$ is conserved, $`P_z`$ accounts also for the degeneracy $`mm`$. Thus, we classify the states as $`1\sigma _{g,u},2\sigma _{g,u},\mathrm{}1\pi _{g,u},2\pi _{g,u}\mathrm{}1\delta _{g,u},2\delta _{g,u}\mathrm{}`$, where the numbers $`1,2,\mathrm{}`$ refer to the electronic states in increasing order of energy. The labels $`\sigma ,\pi ,\delta \mathrm{}`$ are used to denote $`|m|=0,1,2\mathrm{}`$, respectively, the label $`g`$ ($`u`$) gerade (ungerade) is assigned to the states of even (odd) parity $`P`$ of the system. At $`\theta =90^o`$ the Hamiltonian still remains invariant under the parity operations $`P`$ and $`P_z`$, while the angular momentum projection is no longer conserved and $`m`$ is no longer a quantum number. The classification in this case is $`1_{g,u}^\pm ,2_{g,u}^\pm ,\mathrm{}`$, where the sign $`+()`$ is used to denote even (odd) $`z`$-parity. Eventually, for arbitrary orientation, only parity under permutations $`P`$ is conserved. In general, we refer to the lowest gerade and ungerade states in our study as $`1_g`$ and $`1_u`$. This is the only unified notation which make sense for all orientations $`0^{}\theta 90^{}`$.

The above recipe (for the symmetric gauge where $`\xi =1/2,d=0`$) was successfully applied in a study of the $`H_2^+`$-ion in a magnetic field for the parallel configuration $`\theta =0^o`$ Lopez:1997  and also for general one-electron linear systems aligned along the magnetic field Lopez-Tur:2000 . In particular, this led to the prediction of the existence of the exotic ions $`H_3^{2+}`$ at $`B10^{10}G`$ and in a linear configuration $`H_4^{3+}`$ at $`B10^{13}G`$ Turbiner:1999 ; Lopez-Tur:2000 . Recently, this recipe was used for the first time to make a detailed study of the spatial configuration $`H_3^{2+}`$ Lopez-Tur:2002 . It was demonstrated that inconsistency between the form of vector potential and a choice of trial functions can lead to non-trivial artifacts like existence of spurious bound states (see Lopez:2000 ).

One of the simplest trial functions for $`1_g`$ state which meets the requirements of our criterion of adequacy is

$`\mathrm{\Psi }_1=e^{\alpha _1(r_1+r_2)}e^{B[\beta _{1x}\xi x^2+\beta _{1y}(1\xi )y^2]}`$ (5)
(cf. Lopez:1997 ; Turbiner:2001 ), where $`\alpha _1`$, $`\beta _{1x}`$ and $`\beta _{1y}`$ are variational parameters and $`\xi `$ is the parameter of the gauge (2). The first factor in the function (5), being symmetric under interchange of the charge centers $`r_1r_2`$, corresponds to the product of two $`1s`$-Coulomb orbitals centered on each proton. It is nothing but the celebrated Heitler-London approximation for the ground state $`1\sigma _g`$. The second factor is the lowest Landau orbital corresponding to the vector potential of the form Eq. (2). So, the function (5) can be considered as a modification of the free field Heitler-London function. Following the experience gained in studies of $`H_2^+`$ without a magnetic field it is natural to assume that Eq. (5) is adequate to describe interproton distances near equilibrium. This assumption will be checked (and eventually confirmed) a posteriori, after making concrete calculations (see Section IV).

The function (5) is an exact eigenfunction in the potential

$$V_1=\frac{^2\mathrm{\Psi }_1}{\mathrm{\Psi }_1}=2\alpha _1^22B[\beta _{1x}\xi +\beta _{1y}(1\xi )]+4B^2[\beta _{1x}^2\xi ^2x^2+\beta _{1y}^2(1\xi )^2y^2]+2\alpha _1^2(\widehat{n}_1\widehat{n}_2)$$

$$+4\alpha _1B\left[\frac{\beta _{1x}\xi x^2+\beta _{1y}(1\xi )y(yy_1)}{r_1}+\frac{\beta _{1x}\xi x^2+\beta _{1y}(1\xi )y(yy_2)}{r_2}\right]2\alpha _1\left[\frac{1}{r_1}+\frac{1}{r_2}\right],$$

where $`y_{1,2}`$ are the $`y`$-coordinates of protons (see Fig.1). The potential $`V_1`$ reproduces the functional behavior of the original potential (3) near Coulombic singularities and at large distances. These singularities are reproduced exactly when $`\beta _{1x}=\beta _{1y}=1/2`$ and $`\alpha _1=1`$.

One can construct another trial function which meets the requirements of our criterion of adequacy as well,

$$\mathrm{\Psi }_2=\left(e^{\alpha _2r_1}+e^{\alpha _2r_2}\right)e^{B[\beta _{2x}\xi x^2+\beta _{2y}(1\xi )y^2]},$$
(6)
(cf. Lopez:1997 ; Turbiner:2001 ). This is the celebrated Hund-Mulliken function of the free field case multiplied by the lowest Landau orbital, where $`\alpha _2`$, $`\beta _{2x}`$ and $`\beta _{2y}`$ are variational parameters. From a physical point of view this function has to describe the interaction between a hydrogen atom and a proton (charge center), and, in particular, models the possible dissociation mode of $`H_2^+`$ into a hydrogen atom plus proton. Thus, one can naturally expect that for sufficiently large internuclear distances $`R`$ this function prevails, giving a dominant contribution. Again this assumption will be checked a posteriori, by concrete calculations (see Section IV).

There are two natural ways to incorporate the behavior of the system in both regimes – near equilibrium and at large distances – into a single trial function. It is to make a linear or a nonlinear interpolation. The linear interpolation is given by a linear superposition

$$\mathrm{\Psi }_{3a}=A_1\mathrm{\Psi }_1+A_2\mathrm{\Psi }_2,$$
(7)
where $`A_1`$ or $`A_2`$ are parameters and one of them is kept fixed by a normalization condition. In turn, the simplest nonlinear interpolation is of the form

$$\mathrm{\Psi }_{3b}=\left(e^{\alpha _3r_1\alpha _4r_2}+e^{\alpha _3r_2\alpha _4r_1}\right)e^{B[\beta _{3x}\xi x^2+\beta _{3y}(1\xi )y^2]},$$
(8)
(cf. Lopez:1997 ; Turbiner:2001 ), where $`\alpha _3`$, $`\alpha _4`$, $`\beta _{3x}`$ and $`\beta _{3y}`$ are variational parameters. This is a Guillemin-Zener function for the free field case multiplied by the lowest Landau orbital. If $`\alpha _3=\alpha _4`$, the function (8) coincides with (5). If $`\alpha _4=0`$, the function (8) coincides with (6).

The most general Ansatz is a linear superposition of the trial functions (7) and (8),

$$\mathrm{\Psi }=A_1\mathrm{\Psi }_1+A_2\mathrm{\Psi }_2+A_3\mathrm{\Psi }_{3b},$$
(9)
where we fix one of the $`A`$’s and let all the other parameters vary. Finally, the total number of variational parameters in (9), including $`R,\xi ,d`$, is fifteen for the ground state. For the parallel configuration, $`\theta =0^{}`$, the parameters $`\xi =0.5,d=0`$ are fixed in advance and also $`\beta _{1x}=\beta _{1y},\beta _{2x}=\beta _{2y},\beta _{3x}=\beta _{3y}`$. Hence the number of free parameters is reduced to ten for the ground state. Finally, with the function (9) we intend to describe the ground state for all magnetic fields where non-relativistic consideration is valid, $`B4.414\times 10^{13}G`$, and for all orientations of the molecular axis.

Calculations were performed using the minimization package MINUIT from CERN-LIB. Numerical integrations were carried out with a relative accuracy of $`10^7`$ by use of the adaptive NAG-LIB (D01FCF) routine. All calculations were performed on a PC Pentium-III $`800`$ MHz.

## IV Results

We carry out a variational study of the system $`(ppe)`$ with infinitely heavy protons in the range of magnetic fields $`0<B<4.414\times 10^{13}G`$, inclinations $`0^o90^o`$, for a wide range of interproton distances. For magnetic fields $`B<10^{11}G`$ the system displays a well-pronounced minimum in the total energy at all inclinations. However, for $`B>10^{11}G`$ at large inclinations the minimum in the total energy disappears, while for small inclinations a minimum continues to exist. This picture describes the domain of existence of the molecular ion $`H_2^+`$. In general, we confirm a qualitative results by Khersonskij Kher  about the non-existence of a minimum at finite distances on the total energy surfaces of the system $`(ppe)`$ at sufficiently strong magnetic fields for some far from parallel orientations. It is worth mentioning that the variational study in Kher  was carried out with a trial function somewhat similar to that of Eq.(6), which, however, does not fully fulfill our criterion of adequacy. The potential corresponding to this function reproduces correctly the original potential near Coulomb singularities and $`\rho ^2`$ growth at large distances. However, it generates growing terms $`\rho `$ which implies a reduction of the rate of convergence of a perturbation theory for which the variational energy represents the first two terms (see the discussion in Tur ). Also, this trial function is not satisfactory from the point of view of gauge invariance. However, in spite of all the above-mentioned deficiencies it led to qualitatively correct picture.

In Figs. 2-5 the total energy $`E_T`$ of the $`(ppe)`$ system as a function of interproton distance $`R`$ for several values of the magnetic field strength and different values of inclination $`\theta `$ is shown. For magnetic fields $`B10^{11}G`$ and for all inclinations $`0^o90^o`$, each plot displays a well-pronounced minimum at $`R=R_{eq}`$, manifesting the existence of the molecular system $`H_2^+`$. For $`B=1a.u.`$ and $`R3.5a.u.`$ (see Fig. 2) our results are similar to the results of Wille:1988 ; Schmelcher  – for fixed $`R`$ the potential energy $`E_T`$ grows with inclination. In general, at large $`R>R_{eq}`$ and for $`\theta >0^o`$ all the curves behave alike: they reveal a maximum and then tend (from above) to the total energy of the hydrogen atom. For $`\theta =0^o`$ the potential curves approach to the asymptotics from below, displaying in general a behavior similar to the field-free case, to a Van der Waals-force-inspired behavior. This behavior is related to the fact that at large $`R`$ the configuration $`H`$-atom + proton appears. The $`H`$-atom has quadrupole moment, $`QB^2`$ (see Tur ; Lozovik ; Turbiner:1983 ; Potekhin ). Hence at large distances the total energy is defined by a quadrupole moment - charge interaction

$$E_T=\frac{eQ(B)P_2(\mathrm{cos}\theta )}{R^3},$$
(10)
where $`P_2`$ is the second Legendre polynomial. At small inclinations $`P_2(\mathrm{cos}\theta )`$ is positive, the total energy is negative, thus corresponding to attraction between the quadrupole and the charge. Therefore, the total energy curve approaches to the asymptotics from below. For large inclinations $`P_2(\mathrm{cos}\theta )`$ is negative and the total energy is positive. Thus, this corresponds to repulsion between quadrupole and charge, and implies an existence of maximum of the total energy for large interproton distances $`R>R_{eq}`$. We observe the maximum in all Figs.2-5. It is worth mentioning that in the calculations Wille:1988 ; Schmelcher  for $`B=1a.u.`$ and $`\theta =90^o`$ (and other inclinations) the maximum was not observed in contradiction to our predictions (see Fig. 2 and also below Fig. 9). Looking at Fig. 2 it is interesting to compare a rate with which potential curves are approaching to the asymptotic total energy at large $`R`$. This asymptotic energy is equal to the total energy of the hydrogen atom, $`E_H=0.6623Ry`$, while $`E_T^{\theta =0^o}(R=8a.u.)=0.6647`$ (from below), $`E_T^{\theta =45^o}(R=8a.u.)=0.6576`$ (from above), $`E_T^{\theta =90^o}(R=8a.u.)=0.6620`$ (from above). Thus, any deviation does not exceed $`1\%`$. There exists a different manner of viewing these results. It can be treated as a demonstration of the quality of our trial function (9) but for the calculation of the total energy of the atom (!).

However, the situation is drastically different for $`B>10^{11}G`$, see Figs. 4-5. There exists a certain critical angle $`\theta _{cr}`$, such that for $`\theta <\theta _{cr}`$ the situation remains similar to that given above – each potential curve is characterized by a well-pronounced minimum at finite $`R`$. With increase of the inclination, at $`\theta \theta _{cr}`$ the minimum in the total energy first becomes very shallow with $`E_T>E_H`$ and ceases to exist at all. We were unable to localize with a confidence the domain in $`R`$ which correspond to a shallow minimum which leads to the possible dissociation $`H_2^+H+p`$ that was predicted in Larsen  as well as being discussed in our previous work Turbiner:2001 . We consider that the prediction of dissociation for large inclinations emerged as an artifact of an improper choice of the gauge fixing (see the discussion above). A detailed study of the transition domain (existence $``$ non-existence) of $`H_2^+`$ is beyond of scope of the present article. In any case, such a study requires much more accurate quantitative techniques as well as a sophisticated qualitative analysis. Schematically, the situation is illustrated in Fig. 6.

It is quite interesting to explore the variation of the vector potential (2) for $`\theta 0^o`$, in particular, the position of the gauge center as a function of interproton distance $`R`$ and magnetic field 666At $`\theta =0^o`$ (parallel configuration) the vector potential (2) remains unchanged, since $`\xi =1/2`$.. In Figs. 7 a,b for $`B=1a.u.`$ and Figs. 8 a,b for $`10^{12}G`$, correspondingly, both the $`\xi `$ and $`d`$ dependence are presented (see (2) and discussion in Sec.III). This dependence is very similar for all magnetic fields studied. It is worthwhile to emphasize that for all the potential curves given the minimum (in other words, the equilibrium position) at $`R=R_{eq}`$ somehow corresponds to a gauge close to the symmetric gauge: $`\xi 1/2`$ 777The value of $`\xi `$ grows with $`B`$ (see Figs.7a,b and below Tables II-III) and $`Y=Z=0(d=0)`$. A similar situation holds for small interproton distances, $`R<R_{eq}`$. However, for large $`R,R>R_{eq}`$ the parameter $`\xi `$ grows smoothly, reaching a maximum near the maximum of the potential curve which we denote by $`R=R_{cr}`$. It then falls sharply to the value $`\xi 1/2`$. In turn, the parameter $`d`$ remains equal 0 up to $`R=R_{cr}`$ (which means the gauge center coincides to the mid-point between protons), then sharply jumps to 1 (gauge center coincides with the position of a proton), displaying a behavior similar to a phase transition. It is indeed a type of phase transition behavior stemming from symmetry breaking: from the domain $`R<R_{cr}`$, where the permutation symmetry of the protons holds and where the protons are indistinguishable, to the domain $`R>R_{cr}`$, where this symmetry does not exist and electron is attached to one particular proton. Such a type of ‘phase transitions’ is typical in chemistry and is called a ‘chemical reaction’. Hence the parameter $`R_{cr}`$ characterizes a distance at which the chemical reaction $`H_2^+H+p`$ starts. Somewhat similar behavior of the gauge parameters has appeared in the study of the exotic $`H_3^{2+}`$-ion Lopez-Tur:2002 .

In Figs. 9-10 the behaviors of $`R_{eq},R_{max},R_{cr}`$ vs inclination at $`B=1a.u.`$ and $`10^{12}G`$ are displayed. The behavior of $`R_{eq}`$ vs $`\theta `$ demonstrates almost no dependence on $`\theta `$ in contrast to both $`R_{max}`$ and $`R_{cr}`$ which drastically decrease with the growth of $`\theta `$. We do not have a reliable physical explanation of this behavior.

The total energy dependence of $`H_2^+`$ (at $`R=R_{eq}`$) as a function of the inclination angle $`\theta `$ for different magnetic fields is shown in Fig. 11. The dotted line corresponds to the $`H`$-atom total energy in the corresponding magnetic field. For weak magnetic fields the hydrogen atom total energy is always higher than that of the $`H_2^+`$-ion. However, for $`B2\times 10^{11}G`$ the situation changes – a minimum of the $`H_2^+`$ total energy for angles $`\theta >\theta _{cr}`$ does not exist any more. Surprisingly, $`\theta _{cr}`$ corresponds approximately to the moment when the total energy of the $`H`$-atom becomes equal to the total energy of the $`H_2^+`$-ion. If the form of vector potential (2) is kept fixed with $`\xi =1/2`$ and $`Y=Z=0(d=0)`$, then a spurious minimum appears; its position is displayed by the dotted curve. However, if the gauge center parameters are released this minimum disappears (see the discussion above). It was an underlying reason for the erroneous statement about the existence of the unstable $`H_2^+`$ ion in this domain with a possibility to dissociate $`H_2^+H+p`$ (see Turbiner:2001 ). For all magnetic fields studied the total energy is minimal at $`\theta =0^o`$ (parallel configuration) and then increases monotonically with inclination in complete agreement with statements of other authors Wille:1988 ; Kher ; Larsen ; Schmelcher .

In a similar way the binding energy $`E_b=BE_T`$, as well as the dissociation energy (affinity to a hydrogen atom) $`E_d=E_HE_T`$ as a function of $`\theta `$ always decreases when changing from the parallel to the perpendicular configuration (see Fig. 11). Such behavior holds for all values of the magnetic field strength studied. Thus, we can draw the conclusion that the molecular ion becomes less and less stable monotonically as a function of inclination angle. This confirms the statement made in Kher ; Wille:1988 ; Larsen ; Schmelcher , that the highest molecular stability of the $`1_g`$ state of $`H_2^+`$ occurs for the parallel configuration. Thus, the $`H_2^+`$ molecular ion is the most stable in parallel configuration.

We extend the validity of this statement to magnetic field strengths $`10^{13}<B4.414\times 10^{13}G`$. It is worth emphasizing that the rate of increase of binding energy $`E_b`$ with magnetic field growth depends on the inclination – it slows down with increased inclination. This effect implies that the $`H_2^+`$-ion in the parallel configuration becomes more and more stable against rotations – the energy of the lowest rotational state increases rapidly with magnetic field (see Table V below and the discussion there).

Regarding the interproton equilibrium distance $`R_{eq}`$, one would naively expect that it would always decrease with inclination (see Fig. 12). Indeed, for all the magnetic fields studied we observe that $`R_{eq}`$ at $`\theta =0^{}`$ is larger than for any $`\theta 0^{}`$ (see below, Tabs. I, II, III). This can be explained as a natural consequence of the much more drastic shrinking of the electronic cloud in the direction transverse to the magnetic field than in the longitudinal one. Actually, for magnetic fields $`B10^{12}G`$ the equilibrium distance $`R_{eq}`$ decreases monotonically with inclination growth until it reaches $`\theta _{cr}`$, as seen in Fig. 12. As mentioned above, if the parameters of the vector potential (2) are kept fixed, $`\xi =1/2`$ and $`Y=Z=0(d=0)`$, a spurious minimum appears and generates anomalous (spurious) $`R_{eq}`$ behavior for $`\theta >\theta _{cr}`$(see Turbiner:2001 ).

In Tabs. I, II and III the numerical results for the total energy $`E_T`$, binding energy $`E_b`$ and equilibrium distance $`R_{eq}`$ are displayed for $`\theta `$=0, $`45^{}`$ and $`90^{}`$, respectively. As seen in Table I, our results for $`\theta =0^{}`$ lead to the largest binding energies for $`B>10^{11}G`$ in comparison with previous calculations. For $`B10^{11}G`$, our binding energies for the parallel configuration appear to be very close (of the order of $`10^410^5`$ in relative deviation) to the variational results by Wille Wille:1988 , which are the most accurate so far in this region of magnetic field strengths. Those results are based on the use of a trial function in the form of a linear superposition of $`500`$ Hylleraas type functions. It is quite striking that our simple trial function (8) with ten variational parameters gives comparable (for $`B10^{11}G`$) or even better (for $`B>10^{11}G`$) accuracy. It is important to reveal the reason why the trial function Wille:1988  fails to be increasingly inaccurate with magnetic field growth for $`B>10^{11}G`$. An explanation of this inaccuracy is related to the fact that in the $`(x,y)`$\- directions the exact wave function decays asymptotically as a Gaussian function, unlike the Hylleraas functions which decay as the exponential of a linear function. The potential corresponding to the function Wille:1988  reproduces correctly the original potential near Coulomb singularities but fails to reproduce $`\rho ^2`$-growth at large distances. This implies a zero radius of convergence of the perturbation theory for which the variational energy represents the first two terms (see discussion in Tur ).

The results for $`\theta =45^{}`$ are shown in Table II, where a gradual shortening of the equilibrium distance is accompanied by an increase of total and binding energies with magnetic field. It is worth noting that the parameter $`\xi `$ evolves from about $`0.5`$ to $`0.93`$ with magnetic field growth, thus changing from the symmetric gauge for weak fields to an almost asymmetric one for strong ones. This phenomenon takes place for all orientations $`0<\theta <\theta _{cr}`$, becoming more and more pronounced with increasing inclination angle (see below). We are unaware of any other calculations for $`\theta =45^{}`$ to compare ours with.

For the perpendicular configuration $`\theta =90^{}`$, the results are presented in Table III. Similar to what appeared for the parallel configuration (see above) our results are again slightly less accurate than those of Wille Wille:1988  for $`B10^{10}G`$, but becoming the most accurate results for stronger fields. In particular, it indicates that the domain of applicability of a trial function in the form of a superposition of Hylleraas type functions, becomes smaller as the inclination grows. The results reported by Larsen Larsen  and by Kappes-Schmelcher Schmelcher  are slightly worse than ours, although the difference is very small. The evolution of the gauge parameters follow a similar trend, as was observed at $`\theta =45^{}`$. In particular, $`\xi `$ varies from $`\xi =0.64`$ to $`\xi =0.98`$ with magnetic field growth from $`B=10^9G`$ to $`B2\times 10^{11}G`$ 888$`\xi =0.5`$ at $`B=0`$. We should emphasize that the results of Larsen Larsen  and Wille Wille:1988  for $`B>10^{11}G`$ do not seem relevant because of loss of accuracy, since the $`H_2^+`$ ion does not exist in this region.

In order to characterize the electronic distribution of $`H_2^+`$ for different orientations we have calculated the expectation values of the transverse $`<\rho >`$ and longitudinal $`<|z|>`$ sizes of the electronic cloud (see Table IV). Their ratio is always limited,

$$\frac{<\rho >}{<|z|>}<1,$$

and quickly decreases with magnetic field growth, especially for small inclination angles. This reflects the fact that the electronic cloud has a more and more pronounced needle-like form oriented along the magnetic line, as was predicted in the classical papers \[Kadomtsev:1971, -Ruderman:1971, \]. The behavior of $`<\rho >`$ itself does not display any unusual properties, smoothly decreasing with magnetic field, quickly approaching the cyclotron radius for small inclinations and large magnetic fields. In turn, the $`<|z|>`$ monotonically decreases with inclination growth.

As already mentioned, the results of our analysis of the parallel configuration of $`H_2^+`$ turned out to be optimal for all magnetic fields studied, being characterized by the smallest total energy. Therefore, it makes sense to study the lowest vibrational and also the lowest rotational state (see Table V). In order to do this we separate the nuclear motion along the molecular axis near equilibrium in the parallel configuration (vibrational motion) and deviation in $`\theta `$ of the molecular axis from $`\theta =0^{}`$ (rotational motion). The vicinity of the minimum of the potential surface $`E(\theta ,R)`$ at $`\theta =0^{},R=R_{eq}`$ is approximated by a quadratic potential, and hence we arrive at a two-dimensional harmonic oscillator problem in the $`(R,\theta )`$-plane. Corresponding curvatures near the minimum define the vibrational and rotational energies (for precise definitions and discussion see, for example, Larsen ). We did not carry out a detailed numerical analysis, making only rough estimates of the order of $`20\%`$. For example, at B$`=10^{12}G`$ we obtain $`E_{vib}=0.276Ry`$ in comparison with $`E_{vib}=0.259Ry`$ given in Lopez-Tur:2000 , where a detailed variational analysis of the potential electronic curves was performed. Our estimates for the energy, $`E_{vib}`$, of the lowest vibrational state are in reasonable agreement with previous studies. In particular, we confirm a general trend of the considerable increase of vibrational frequency with the growth of $`B`$ indicated for the first time by Larsen Larsen . The dependence of the energy on the magnetic field is much more pronounced for the lowest rotational state – it grows much faster than the vibrational one with magnetic field increase. This implies that the $`H_2^+`$-ion in the parallel configuration becomes more stable for larger magnetic fields (see the discussion above). From a quantitative point of view the results obtained by different authors are not in good agreement. It is worth mentioning that our results agree for large magnetic fields $`10a.u.`$ with results by Le Guillou et al. Legui , obtained in the framework of the so called ‘improved static approximation’, but deviate drastically at $`B=1a.u.`$, being quite close to the results of Larsen Larsen  and Wille Wille:1987 . As for the energy of the lowest rotational state, our results are in good agreement with those obtained by other authors (see Table V).

In Fig. 13 we show the electronic distributions $`𝑑y|\psi (x,y,z)|^2`$ for magnetic fields $`10^9,10^{10},10^{11},10^{12}G`$ and different orientations for $`H_2^+`$ in the equilibrium configuration, $`R=R_{eq}`$. It was already found explicitly Lopez:1997  that at $`\theta =0^o`$ with magnetic field increase there is a change from ‘ionic’ (two-peak electronic distribution) to ‘covalent’ coupling (one-peak distribution). We find that a similar phenomenon holds for all inclinations. If for $`B=10^9G`$, all electronic distributions are characterized by two peaks for all inclinations, then for $`B=10^{12}G`$ all distributions have a single sharp peak. The ‘sharpness’ of the peak grows with magnetic field. Fig. 13 also demonstrates how the change of the type of coupling appears for different inclinations – for larger inclinations a transition $`(twopeaks)(onepeak)`$ appears for smaller magnetic fields. It seems natural that for the perpendicular configuration $`\theta =90^{}`$, where the equilibrium distance is the smallest, this change appears for even smaller magnetic field.

In Figs. 14-18 we present the evolution of the electronic distributions as a function of interproton distance $`R`$, for inclinations $`0^o,45^o`$ at $`B=1a.u.`$ and $`10^{12}G`$ together with the $`R`$-dependence for the inclination $`90^o`$ at $`B=1a.u.`$. The values of magnetic fields are chosen to illustrate in the most explicit way the situation. In all figures a similar picture is seen. Namely, at not very large magnetic fields $`B10^{11}G`$ and for all inclinations the electronic distribution at small $`R<R_{cr}`$ is permutationally symmetric and evolves with increase of $`R`$ from a one-peak to a two-peak picture with more and more clearly pronounced, separated peaks. Then for $`R=R_{cr}`$ this symmetry becomes broken and the electron randomly chooses one of protons and prefers to stay in its vicinity. For $`RR_{cr}`$ the electronic distribution becomes totally asymmetric, the electron looses its memory of the second proton. This signals that the chemical reaction $`H_2^+H+p`$ has already happened. For larger magnetic fields $`B10^{11}G`$ for $`R<R_{cr}`$ the electronic distribution is always single-peaked, a transition from a one-peak to a two-peak picture occurs for $`R>R_{cr}`$, where the electronic distribution is already asymmetric.

To complete the study of the $`1_g`$ state we show in Fig. 19 the behavior of the variational parameters of the trial function (9) as a function of the magnetic field strength for the optimal (parallel) configuration, $`\theta =0^o`$. In general, the behavior of the parameters is rather smooth and very slowly-changing, even though the magnetic field changes by several orders of magnitude. This is in drastic contrast with the results of Kappes-Schmelcher Kappes:1994  (see Fig. 1 in this paper). In our opinion such behavior of the parameters of our trial function (9) reflects the level of adequacy (or, in other words, indicates the quality) of the trial function. In practice, the parameters can be approximated by the spline method and then can be used to study magnetic field strengths other than those presented here.

## V Conclusion

We have carried out an accurate, non-relativistic calculation in the Born-Oppenheimer approximation for the lowest state of the $`H_2^+`$ molecular ion for different orientations of the magnetic field with respect to the molecular axis. We studied constant uniform magnetic fields ranging from zero up to $`B=4.414\times 10^{13}G`$, where non-relativistic consideration holds.

For all magnetic fields studied there exist a region of inclinations for which a well-pronounced minimum in the total energy surface for the $`1_g`$ state of the system $`(ppe)`$ is found. This shows the existence of the $`H_2^+`$ molecular ion for magnetic fields $`B=04.414\times 10^{13}G`$. The smallest total energy is always found to correspond to the parallel configuration, $`\theta =0^{}`$, where protons are situated along the magnetic line. The total energy increases, while the binding energy decreases monotonically as the inclination angle grows. The rate of total energy increase as well as binding energy decrease is seen to be always maximal for the parallel configuration. The equilibrium distance exhibits quite natural behavior as a function of the orientation angle $`\theta `$ – for fixed magnetic field the shorter equilibrium distance always corresponds to the larger $`\theta `$.

Confirming the qualitative observations made by Khersonskij Kher  for the $`1_g`$ state in the contrast to statements in Larsen ; Wille:1988 , we demonstrate accurately that the $`H_2^+`$-ion does not exist at a certain range of orientations for magnetic fields $`B2\times 10^{11}G`$. As the magnetic field increases the region of inclinations where $`H_2^+`$ does not exist is seen to broaden, reaching rather large domain $`25^{}\theta 90^{}`$ for $`B=4.414\times 10^{13}G`$.

We find that the electronic distributions for $`H_2^+`$ in the equilibrium position are qualitatively different for weak and large magnetic fields. In the domain $`B<10^{10}G`$ the electronic distribution for any inclination has a two-peak form, peaking near the position of each proton. On the contrary for $`B>10^{11}G`$ the electronic distribution always has a single peak form with the peak near the midpoint between the protons for any inclination. This implies a physically different structure for the ground state - for weak fields the ground state can be modelled as a ‘superposition’ of hydrogen atom and proton, while for strong fields such modelling is not appropriate.

Unlike standard potential curves for molecular systems in the field-free case, we observe for $`\theta >0^o`$ that each curve has a maximum and approaches to the asymptotics at $`R\mathrm{}`$ from above. The electronic distribution evolves with $`R`$ from a one-peak form at small $`R`$ to a two-peak one at large $`R`$. There exists a certain critical $`R_{cr}`$ at which one of peaks starts to diminish, manifesting a breaking of permutation symmetry between the protons and simultaneously the beginning of the chemical reaction $`H_2^+H+p`$.

Combining all the above-mentioned observations we conclude that for magnetic fields of the order of magnitude $`B10^{11}G`$ some qualitative changes in the behavior of the $`H_2^+`$ ion take place. The behavior of the variational parameters also favors this conclusion. This hints at the appearance of a new scale in the problem. It might be interpreted as a signal of a transition to the domain of developed quantum chaos (see, for example, chaos ).

###### Acknowledgements.

The authors wish to thank B.I. Ivlev and M.I. Eides for useful conversations and interest in the subject. We thank M. Ryan for a careful reading of the manuscript. This work was supported in part by DGAPA Grant # IN120199 (México) and CONACyT Grant 36650-E (México). AT thanks the University Program FENOMEC for financial support.