domain: chim phys chim.theo sde math

text:
Alexander Alijah
email: alexander.alijah@univ.reims.fr

Juan Carlos

López Vieyra
email: vieyra@nucleares.unam.mx

Daniel J Nader
email: daniel.nader@correo.nucleares.unam.mx

Alexander V Turbiner
email: turbiner@nucleares.unam.mx

Héctor Medel Cobaxin

Carlos López Vieyra

The hydrogen molecule H2 in inclined configuration in a weak magnetic field

Keywords: variational method, weak magnetic field, critical magnetic field, magnetic susceptibility, ro-vibrational states

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

More than fifty years have passed since it was predicted that extremely strong magnetic fields up to B = 10 14 -10 16 G (B ∼ 4 × 10 4-6 a.u.), which are by far beyond those that can be reached in the laboratory, could exist in neutron stars remnant of a supernova explosion as effect of the magnetic field flux compression [START_REF] Woltjer | X-rays and type I supernova remnants[END_REF] (see also [START_REF] Pacini | Energy emission from a neutron star[END_REF][START_REF] Gold | Rotating neutron stars as the origin of the pulsating radio sources[END_REF][START_REF] Goldreich | Pulsar electrodynamics[END_REF]). As for magnetized white dwarfs, the surface magnetic field can reach B ∼ 10 9 G (see e.g. [START_REF] García-Berro | Magnetic white dwarfs: Observations, theory and future prospects[END_REF] and references therein). Soon afterwards it was recognized that the structure of atoms and molecules might be qualitatively different under strong magnetic fields B B 0 (B 0 = 1 a.u. ≡ 2.35×10 9 G = 2.35×10 5 T) [START_REF] Ruderman | Matter in superstrong magnetic fields: the surface of a neutron star[END_REF][START_REF] Kadomtsev | Molecules in an ultrastrong magnetic field[END_REF][START_REF] Kadomtsev | Matter in a superstrong magnetic field[END_REF] from the field-free case. The electronic clouds assume a well-pronounced cigar-like form, and molecules become oriented along the magnetic line. Eventually, the problem becomes quasi-one-dimensional, where longitudinal and transverse motions of the electrons are almost separated. This gives hope to develop an analytical theory in the domain of very strong magnetic fields. The situation gets much more complicated in the domain of intermediate magnetic fields, say, of order of B ∼ 10 -1 a.u., where quadratic corrections to the linear Zeeman effect become significant. This domain is 'slightly' above the magnetic fields reachable in the laboratory. In this case we do not see hope to develop analytic approaches. We will call the fields 0.01 B 1 a.u. the intermediate magnetic fields.

Due to mainly technical difficulties in solving the Schrödinger equation in the presence of intermediate and strong magnetic fields, only a relatively small number of simple atomic and molecular systems has been studied. Naturally, the hydrogen atom H and the hydrogen molecular ion H + 2 are the most studied systems, see e.g. [START_REF] Kravchenko | Exact solution for a hydrogen atom in a magnetic field of arbitrary strength[END_REF] and [START_REF] Turbiner | One-electron molecular systems in a strong magnetic field[END_REF], respectively, and references therein. The first quantitative study of the H 2 molecule was carried out by one of the present authors in 1983 [START_REF] Turbiner | Hydrogen molecule in a strong magnetic field[END_REF]. In the majority of studies of molecules and molecular ions all non-adiabatic terms in the Hamiltonian are neglected by assuming an infinite nuclear mass (what is usually called Born-Oppenheimer (BO) approximation of zeroth order). The fact is that in both the H and H + 2 systems the binding energy grows dramatically with an increase of the magnetic field strength, it hints at the possible existence of unusual chemical species in strong magnetic fields. Other simple, traditional such as H + 3 [START_REF] Turbiner | H + 3 molecular ion in a magnetic field: Linear parallel configuration[END_REF][START_REF] Medel Cobaxin | Vibrating H + 3 in a uniform magnetic field[END_REF], and exotic compounds mainly formed by protons and/or α-particles (helium nuclei) and one or two electrons have been studied to a certain degree. For a discussion, see [START_REF] Turbiner | One-electron molecular systems in a strong magnetic field[END_REF] for one-electron systems and [START_REF] Turbiner | Charged hydrogenic, helium, and helium-hydrogenic molecular chains in a strong magnetic field[END_REF] for two-electron systems.

Recently, a detailed study of the H + 2 molecular ion in inclined configuration (when the molecular axis and the magnetic line form some non-zero angle) was carried out for intermediate and strong magnetic fields [START_REF] Turbiner | One-electron molecular systems in a strong magnetic field[END_REF][START_REF] Medel Cobaxin | H 2 + in a weak magnetic field[END_REF]. It was shown that for the ground state the optimal configuration is always parallel, where the molecular axis and magnetic field direction coincide. The spectra of rovibrational states was exhaustively studied.

As for the H 2 molecule, it was found long ago that the minimal energy (ground) state evolves with magnetic field strength being realized by different states depending on the strength of the magnetic field, see [START_REF] Detmer | Hydrogen molecule in magnetic fields: The ground states of the Σmanifold of the parallel configuration[END_REF][START_REF] Detmer | Hydrogen molecule in a magnetic field: The lowest states of the Π-manifold and the global ground state of the parallel configuration[END_REF] and references therein. At zero and weak magnetic fields, the H 2 ground state is realized by the spin-singlet S = 0, 1 Σ g state in parallel configuration, but with the magnetic field strength increasing to above the critical field strength of B cr = 0.178 a.u., see below, the ground state changes to a spin-triplet S = 1, 3 Σ u , state which is a repulsive state (!). It corresponds to two hydrogen atoms at large distances with electron spins antiparallel to the magnetic field, hence, the hydrogen molecule does not exist as a compact system. It is worth noting that this value of the critical magnetic field was calculated accurately in present paper and confirms the rough estimate B cr 0.2 a.u. from [START_REF] Detmer | Hydrogen molecule in magnetic fields: The ground states of the Σmanifold of the parallel configuration[END_REF]. Nonetheless, for even stronger magnetic fields, B 12 a.u., the ground state is realized by a spin triplet S = 1, 3 Π u state, see [START_REF] Detmer | Hydrogen molecule in a magnetic field: The lowest states of the Π-manifold and the global ground state of the parallel configuration[END_REF] and references therein. A similar behavior is observed in the case of the linear H + 3 molecular ion in strong magnetic fields: the ground state evolves from the spin-singlet 1 Σ g state for weak magnetic fields B 5×10 8 G ( 0.2 a.u.) to a weakly bound spin-triplet 3 Σ u state for intermediate and strong fields and, eventually, to a spin-triplet 3 Π u state for magnetic fields B 5 × 10 10 G ( 21 a.u.) [START_REF] Turbiner | H + 3 molecular ion in a magnetic field: Linear parallel configuration[END_REF]. In such studies the parallel configuration of the molecular axis and the magnetic field direction is explicitly assumed.

Non-aligned configurations, where the molecular axis is not parallel to the magnetic field direction, have received much less attention. This is due to the fact that such configurations require a much larger computational effort to reach the accuracies obtained in the parallel case. The present authors are not aware of any studies of inclined configurations for the H 2 molecule for B 0.2 a.u.

The goal of this paper is to study the hydrogen molecule H 2 arbitrarily oriented i.e. with the molecular axis forming an angle θ with respect to the direction of a uniform magnetic field B in lowest spin-singlet state 1 g . The magnetic field strengths of interest in this work are chosen to be B = 0, 0.1, 0.175 and 0.2 a.u. (equivalently, 0, 2.35 × 10 8 , 4.1 × 10 8 and 4.7 × 10 8 G), where the ground state is realized by the spin-singlet state 1 Σ g at θ = 0 for B < B cr . We use the variational method with trial functions designed following a criterion of physical adequacy [START_REF] Turbiner | The eigenvalue spectrum in quantum mechanics and the nonlinearization procedure[END_REF][START_REF] Turbiner | One-electron molecular systems in a strong magnetic field[END_REF]. Three inclinations θ = 0 • , 45 • and 90 • will be studied in detail and the potential energy curves for each inclination and each magnetic field will be constructed. The two-dimensional potential energy surfaces are obtained by interpolation in the θ coordinate. This allows us to calculate for the first time the lowest rovibrational levels of the H 2 molecule in weak and intermediate magnetic fields, where this molecule exists as a compact object. A study of the magnetic susceptibility of the H 2 molecule is also performed. We will follow in presentation our previous work on H + 2 in weak and intermediate magnetic fields [START_REF] Medel Cobaxin | H 2 + in a weak magnetic field[END_REF]. Atomic units will be used through the text.

The Hamiltonian and generalities

We consider the hydrogen molecule H 2 interacting with an external magnetic field B. The origin of coordinates is chosen in the midpoint of the line connecting the nuclei (molecular axis). The molecular axis in turn forms an angle θ with respect to the magnetic field direction (chosen to coincide with the z-axis). A convenient gauge which describes a magnetic field oriented parallel to the z-axis, is the linear gauge

 = B[(ξ -1)y, ξx, 0], ( 1
)
where ξ is a parameter. If ξ = 0 the linear gauge is reduced to the Landau gauge, and if ξ = 1/2 then the symmetric gauge is obtained. In approximate variational calculations the parameter ξ is considered as an extra variational parameter.

Since the nucleus mass is by far larger than the electron mass, we can neglect all non-adiabatic coupling terms in the Hamiltonian to obtain the order zero BO approximation. Thus, the electronic Hamiltonian in atomic units ( = m e = c = 1) is given by Ĥe = -1 2

2 i=1 ∇ 2 i -iB 2 i=1 ((ξ -1)y i ∂ x i + ξx i ∂ y i ) + S • B + 1 2 B 2 2 i=1 ξ 2 x 2 i + (ξ -1) 2 y 2 i - 2 i=1 1 r ia + 1 r ib + 1 r 12 + 1 R , (2)
where ∇ i is the Laplacian operator with respect to the coordinates of the i-th electron r = (x i , y i , z i ), r ia,ib are the distances between the i-th electron and the nuclei a or b, respectively, r ij is the distance between the electrons and R is the distance between the nuclei. As usual, the contribution to the energy due to the Coulomb interaction between the nuclei (1/R) is treated classically. Hence, R is considered an external parameter. In the particular case θ = 0 • , the component of the angular momentum along the z-axis is conserved and the term linear in B becomes 1 2 L • B for ξ = 1 2 . The spin Zeeman term S • B with the total electron spin S = S 1 + S 2 is included in the Hamiltonian. However, for the spin-singlet states with S = 0 this term does not contribute to the total energy and can be excluded.

Finally, the nuclear motion can be treated as vibrations and rotations following the BO approximation with the electronic energy taken as the potential in the nuclear Hamiltonian.

The trial function

Following physical relevance arguments (see, e.g. [START_REF] Turbiner | The eigenvalue spectrum in quantum mechanics and the nonlinearization procedure[END_REF]) we designed a spatial trial function which is a product of Landau orbitals, Coulomb orbitals and a correlation term in exponential form:

ψ(r 1 , r 2 ) = 2 k=1 (e -α ka r ka -α kb r kb -Bβ kx 4 x 2 k - Bβ ky 4 y 2 k )e α 12 r 12 (3)
where α ka,kb , β kx , β ky with k = 1, 2 as well as α 12 are variational parameters. In (3) the variational parameters α ka , α kb (k = 1, 2) have the meaning of screening (or anti-screening) factors (charges) for the nucleus a, b respectively, as it is seen from the k-th electron. The variational parameters β kx , β ky account for the screening (or anti-screening) factors for the magnetic field seen from k-th electron in x, y direction respectively, and the parameter α 12 "measures" the screening (or anti-screening) of the electron correlation interaction. This spatial function reproduces adequately the behavior of the electrons near the Coulomb singularities and the harmonic oscillator at long distances arising from the magnetic field. In a certain way the trial function (3) is a generalization of the trial function presented in [START_REF] Turbiner | A note about the ground state of the hydrogen molecule[END_REF] for the field free case. It reproduces two physical situations: for small internuclear distances the trial function (3) mimics the interaction H + 2 + e (if α 1a = α 1b and α 2a = α 2b ) while for large internuclear distances it mimics the interaction H -H (if α 1a = α 2b and α 1b = α 2a ).

We consider a trial function which is a superposition of three Ansätze: a general Ansatz of the type (3), a H -H type Ansatz and a H

+ 2 + e type Ansatz Ψ = A 1 ψ + A 2 ψ H+H + A 3 ψ H + 2 +e , (4)
where A 1,2,3 are linear variational parameters. Each Ansatz has its own set of variational parameters. Without loss of generality A 1 may be set equal to the unity, therefore the total number of variational parameters is 27 including the internuclear distance R and ξ as variational parameters.

In the singlet state (S = 0) the trial function (4) must be symmetric with respect to the exchange of the electrons and in the gerade (g) state the trial function (4) must be symmetric with respect to the exchange of nuclei. Therefore the operator

(1 + Pab )(1 + P12 ) , (5)
where Pab is the operator of symmetrization of nuclei and P12 is the operator of symmetrization of the electrons, must be applied to the trial function [START_REF] Goldreich | Pulsar electrodynamics[END_REF]. The calculation of the variational energy using the trial function (3) involves two major parts: (i) 6-dimensional numerical integrations which were implemented by an adaptive multidimensional integration C-language routine (cubature) [START_REF] Genz | Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N -dimensional rectangular region[END_REF], and (ii) a minimizer which was implemented with the Fortran minimization package MINUIT from CERN-LIB. Our C-Fortran hybrid program was parallelized using MPI. The 6-dimensional integrations were carried out using a dynamical partitioning procedure: the domain of integration is manually divided into sub-domains following the profile of the integrand. Then each sub-domain is integrated on separated processors using the routine CUBATURE. In total, we have a division into 960 subregions for the numerator and ∼ 1000 for the denominator of the variational en-ergy. With a maximal number of sampling points ∼ 10 8 for the numerical integrations for each subregion, the time needed for one evaluation of the variational energy (two integrations) is 2 × 10 3 seconds (∼ 37 min) with 96 processors at the cluster KAREN (ICN-UNAM, Mexico). It was checked that this procedure stabilizes the estimated accuracy and is reliable in the first three to four decimal digits. However, in order to localize the domain of minimal parameters, a minimization procedure with much less sample points was used in each sub-domain, and a single evaluation of the energy usually took ∼ 15 -20 mins. Once a domain is roughly localized, the number of sample points is increased by a factor of ∼ 10 2 . Typically, a minimization process required several hundreds of evaluations. As a general strategy, the variational energy corresponding to the general Ansatz only is calculated in first place. Then, either the H -H type Ansatz or the H + 2 + e type Ansatz is added as a first correction, depending on which configuration yields a better variational result, and the energy is minimized using the superposition of two Ansätze. Eventually, the remaining configuration is included in the final trial function and a final minimization is carried out. The whole process is very lengthy and cumbersome due to the absence of a fast minimization procedure. Computations were mainly performed in parallel on 96 processors on the cluster ROMEO at the University of Reims, France, and on the cluster KAREN at ICN-UNAM, Mexico.

Results

The electronic energies and the equilibrium distances of H 2 in the 1 g state are presented in Table 1 for magnetic fields B = 0, 0.1, 0.175 and 0.2 a.u. Variational energies indicate that for B ≤ B cr = 0.178 a.u. the lowest energy state of H 2 is realized by the 1 g state in parallel configuration. For B cr = 0.178 a.u. the energy of the 1 g state at the equilibrium minimum coincides with the energy of two hydrogen atoms infinitely separated and having both electron spins antiparallel to the magnetic field direction. Thus, for B = 0.2 a.u. the state 1 g is, in fact, a m-eta-stable state. We studied the geometrical configurations with angles θ = 0 • , 45 • and 90 • between the magnetic field direction and molecular axis in great detail, while some sample calculations were carried out for the intermediate angles θ = 15 • , 30 • and 60 • , 75 • to check the smoothness of the angular dependence of both, the energy and the equilibrium distance (see below). For all inclinations the potential energy curve E vs R exhibits a well pronounced minimum at a finite internuclear distance R. As the magnetic field increases, for any given inclination the system becomes more strongly bound (the binding energy increases) and more compact (the internuclear equilibrium distance R eq is reduced), see Table 1. Note that for the field-free case B = 0 our energy is in agreement with one of the most accurate results [START_REF] Sims | High precision variational calculations for the Born-Oppenheimer energies of the ground state of the hydrogen molecule[END_REF] in ∼ 2 × 10 -4 a.u. in spite of the very simple form of the trial function that we used. We must emphasize that for parallel configuration θ = 0 our energies are systematically better than the ones from [START_REF] Detmer | Hydrogen molecule in magnetic fields: The ground states of the Σmanifold of the parallel configuration[END_REF] in 3 decimal digits (d.d.), which leads to a more accurate value of the critical magnetic field strength B cr . For a given magnetic field, the total energy increases while the equilibrium distance R eq shows a small decrease with growth of the inclination angle from θ = 0 to 90 • , see Figures 1 and2. Such an increase in E, and decrease in R eq , are more pronounced as the magnetic field increases. Thus, for all magnetic fields studied, the optimal configuration corresponds to the parallel configuration as it is expected. The angular dependence of the variational energy E(B, θ) and the equilibrium distance R eq (B, θ) for a fixed magnetic field strength B is very simple and is well-described by the hindered rotator model, see Eq. ( 15) and captions of Figs. 1 and2. This observation is in agreement with the test calculations for angles θ = 15 • , 30 • , 60 • and 75 • .

Potential Energy Curves

Potential energy curves E vs R of the state 1 g of the H 2 molecule in magnetic fields B = 0, 0.1, 0.2 a.u. and inclinations θ = 0, 45 • , 90 • are built from variational results obtained in the domain R ∈ [1, 2] a.u. and extended beyond following the procedure discussed in [START_REF] Olivares-Pilón | H 2 + , HeH and H 2 : Approximating potential curves, calculating rovibrational states[END_REF] for approximating potential curves in diatomic molecules (see also references therein). It is evident that the asymptotic behavior of the electronic energy of H 2 at small distances R → 0 is given by

E ≈ 1 R + E He (B) + c 1 R + O(R 2 ) , (6)
where E He (B) is the ground state energy of the helium atom in a magnetic field (B) (the so-called united atom limit), and the coefficient in front of R depends on the magnetic field and the inclination θ, c 1 = c 1 (B, θ); at B = 0 this coefficient vanishes c 1 = 0 (see [START_REF] Olivares-Pilón | H 2 + , HeH and H 2 : Approximating potential curves, calculating rovibrational states[END_REF] and references therein). As for the asymptotic limit R → ∞, the expansion of the energy E is given by where E 2H (B) is the energy of two (infinitely separated) hydrogen atoms in their ground state in the magnetic field of strength B (however, with opposite electron spin projections so that S • B = 0), the term ∝ 1/R 5 corresponds to the quadrupole-quadrupole interaction (repulsive for 0, 90 • and attractive for 45 • ) between two separated hydrogen atoms in the magnetic field (which is the leading order interaction at R → ∞). The term ∝ 1/R 6 corresponds to the induced dipole-dipole interaction (in second order perturbation theory in 1/R for B = 0) between two separated hydrogen atoms (see [START_REF] Turbiner | Hydrogen molecule in a strong magnetic field[END_REF] and [START_REF] Kadomtsev | Matter in a superstrong magnetic field[END_REF]). The coefficients c 5,6,7 can depend on the magnetic field strength and inclination c 5,6,7 = c 5,6,7 (B, θ). In absence of a magnetic field c 5,7 = 0. In general, the quadrupole-quadrupole interaction energy (in a.u.) is given by

E ≈ E 2H (B) + c 5 R 5 - c 6 R 6 + c 7 R 7 + O 1 R 8 , (7)
E Q = 3 4 Q 2 zz (B) P 4 (cos θ) R 5 , (8)
where Q zz is the quadrupole moment of the hydrogen atom in a magnetic field of strength B (see [START_REF] Turbiner | Hydrogen molecule in a strong magnetic field[END_REF]), P 4 is 4th Legendre polynomial. Thus, the coefficient c 5 is known. For weak magnetic fields B we use the approximation the quadrupole moment in perturbation theory (see [START_REF] Potekhin | Hydrogen atom in a magnetic field: The quadrupole moment[END_REF])

Q zz = - 5 2 B 2 + 615 32 B 4 + . . . . (9)
Now we interpolate both asymptotic expansions ( 6) and ( 7) via the twopoint Padé approximant P ade[N/N + 4](R) with N = 2 as the minimal degree, which guarantees that the expansions ( 6) and ( 7) are described functionally correct,

E(R) = 1 R a 0 + a 1 R + a 2 R 2 (b 0 + b 1 R + b 2 R 2 + b 3 R 3 + b 4 R 4 + b 5 R 5 + b 6 R 6 ) + E 2H (B) , (10)
where the constraints

b 0 = a 0 , b 1 = a 0 (E 2H (B) -E He (B)) + a 1 ,

Magnetic Susceptibility

The trial function (4), in spite of its simplicity, incorporates accurately the many physics features of the H 2 molecule in a magnetic field. In order to verify this assertion for weak magnetic fields we calculated the magnetic susceptibility (magnetizability). To make this calculation we follow the recipe proposed in our work on H + 2 [START_REF] Medel Cobaxin | H 2 + in a weak magnetic field[END_REF]. It is well known that the response of the molecule to an external magnetic field falls into two parts: diamagnetic and paramagnetic. Correspondingly, there are two contributions to the magnetic susceptibility: a paramagnetic, χ p , originating from the linear Zeeman term in the Hamiltonian (2) when treated within second order perturbation theory in powers B, and a diamagnetic, χ d , coming from the quadratic Zeeman term ∼ B 2 in the first order of perturbation theory. Thus, the total magnetic susceptibility is the sum of the two terms χ = χ d + χ p .

In general, the magnetic susceptibility tensor χ αβ is defined by the coef-

H = - 1 2 α,β χ αβ B α B β ,
with B α , B β being the components of the magnetic field. As for the diamagnetic susceptibility it is given in first order PT in B 2 as

χ d αβ = - 1 4 2 i=1 r 2 i δ αβ -r i,α r i,β , (11)
where both the expectation value r 2 i of the position vector squared of the i-th electron and the 2nd order tensor r i,α r i,β , α, β = x, y, z are taken with respect to the field-free wavefunction at equilibrium distance R eq . If the magnetic field direction is chosen along the z-axis, B = B ẑ the tensor χ d αβ appears in diagonal form and contains a single non-zero component, χ d zz ≡ χ d ,

χ d = - 1 4 2 i=1 ( x 2 i + y 2 i ) , (12)
where the symmetric gauge is assumed to be taken. On the other hand, the paramagnetic contribution to the susceptibility is much more difficult to calculate, since it occurs in second order PT. In general, the paramagnetic susceptibility is much smaller than the diamagnetic one. In principle, this contribution to the susceptibility can be easily evaluated as the difference χ p = χ-χ d , where χ is the total magnetic susceptibility at a given inclination.

As for the ground state, the total magnetic susceptibility can be calculated in a straightforward way as the coefficient in front of the B 2 term in the energy expansion

E(B, θ) = E(0) - 1 2 χ(θ)|B| 2 + . . . , (13)
at R = R eq . The results for the susceptibilities are presented in Table 3 for inclinations θ = 0 • , 45 • , 90 • , they are compared with the experimental data from Ramsey [START_REF] Ramsey | Molecular Beams, International series of monographs on physics[END_REF], and with other calculations, when available. In general, all susceptibilities grow with the inclination angle. For θ = 0 our χ d are larger than the values obtained in the past in [START_REF] Riley | The octopole magnetizability of the hydrogen molecule[END_REF][START_REF] Kolos | Potential-energy curves for the X 1 Σ + g , b 3 Σ + u , and C 1 Π u states of the hydrogen molecule[END_REF]. They are closer to experimental data being different from experimental data in one portion ×10 -3 in spite of the fact that our trial function is much simpler than the ones used in [START_REF] Riley | The octopole magnetizability of the hydrogen molecule[END_REF] and [START_REF] Kolos | Potential-energy curves for the X 1 Σ + g , b 3 Σ + u , and C 1 Π u states of the hydrogen molecule[END_REF]. As for θ = 45 • , the susceptibilities are calculated for the first time to the best of our knowledge, while for θ = 90 • our χ d agrees in 2 d.d. with [START_REF] Kolos | Potential-energy curves for the X 1 Σ + g , b 3 Σ + u , and C 1 Π u states of the hydrogen molecule[END_REF] and differs from experimental data in 2 × 10 -2 . Concerning χ p , it is superior to all nine values calculated previously and collected in Table XII of [START_REF] Rychlewski | Electric and magnetic properties for the ground and excited states of molecular hydrogen[END_REF], however, it still differs from experimental data in ∼ 20%. In general, our results for the susceptibility agree very well with the experimental data and with other calculations.

Rovibrational levels

The lowest rovibrational states of H 2 and D 2 were calculated for the field strengths B = 0.1 B 0 and B = 0.2 B 0 , where B 0 = 2.35 × 10 9 Gauss = 2.35×10 5 T, as described in [START_REF] Medel Cobaxin | H 2 + in a weak magnetic field[END_REF]. To keep the present paper self-contained, the method is briefly summarized below. Starting point is the nuclear Hamiltonian expressed in spherical coordinates,

Ĥnuc = - 2 M s 1 R ∂ 2 ∂R 2 R + 2 M s R 2 L2 R - 1 M s B Lz + 1 8M s B 2 R 2 sin 2 θ + Ṽ (R, θ). (14)
Here, M s denotes the total mass of the nuclei, Lz is the projection of angular momentum along z-axis and θ the angle between the molecular and the zaxis. The two-dimensional potential, Ṽ (R, θ), is parametrized as a hindered rotator, where only the lowest expansion term is maintained, to yield (in a.u.) are also included for B = 0 at the equilibrium distance R eq = 1.40 a.u., they were obtained using the trial function (4). exp Experimental results from [START_REF] Ramsey | Molecular Beams, International series of monographs on physics[END_REF], see also Table I in [START_REF] Pflug | Calculation of the magnetic susceptibility of H 2 by the distinguishable electron method[END_REF]. Results marked as a are from [START_REF] Riley | The octopole magnetizability of the hydrogen molecule[END_REF]

Ṽ (R, θ) = Ṽ (R, 0) + n V 90,n (R) 2 [1 -cos(2nθ)]
≈ Ṽ (R, 0) + V 90 (R) sin 2 θ ( 15
)
V 90 (R) = Ṽ (R, 90) -Ṽ (R, 0) is the barrier height for a given value of R.

The rovibrational wave function can be expanded in terms of vibrational and rotational basis functions as

Ψ(R, θ, φ) = v,L c v,L ξ v (R; θ ) R Y M L (θ, φ) (16)
where ξ v (R; θ ) are solutions of the vibrational part of Eq. ( 14) at the reference orientation θ , chosen as θ = 0. These are obtained numerically using the renormalized Numerov algorithm. The Y M L (θ, φ) in the above equation are spherical harmonics.

In this basis, the matrix elements of the Hamiltonian in Eq. ( 14) are given by

v L M | Ĥnuc |vLM = E v δ L L δ v v + 2 M s v | 1 R 2 |v L(L + 1)δ L L - BM M s δ L L δ v v + B 2 12M s v |R 2 |v + 2 3 v |V 90 (R)|v δ L L - B 2 12M s v |R 2 |v + 2 3 v |V 90 (R)|v × (-1) M (2L + 1)(2L + 1) × L 2 L 0 0 0 L 2 L M 0 -M (17)
The terms in parentheses are Wigner 3j-symbols. The matrix Eq. ( 17) is diagonal in M as expected, since M is an exact quantum number. L-functions are coupled in steps of 2, conserving z-parity, π = (-1) L+M . Diagonalization of the Hamiltonian matrix, Eq. ( 17), yields the eigenvalues and eigenvectors of the rovibrational problem.

We have computed the lowest rovibrational states for H 2 and D 2 . Allowed rovibrational states must obey the permutational symmetry of the two identical nuclei. In the case of H 2 , with two fermions, the symmetry of the vibrational and rotational parts of the rovibrational wavefunction must be opposite, while in the case of D 2 , with two bosons, it must be the same if we consider ortho nuclear spins. For a rovibrational state of given vibrational quantum number, v, and projection of the angular momentum on the magnetic field axis, M , the z-parities are thus π = (-1) M +v+1 = -(-1) M for v even (-1) M for v odd [START_REF] Turbiner | The eigenvalue spectrum in quantum mechanics and the nonlinearization procedure[END_REF] for H 2 , and

π = (-1) M +v = (-1) M for v even -(-1) M for v odd (19) for D 2 .
The results of our calculations for the lowest vibrational states, v = 0, 1, 2, 3 and M ≤ 5 are presented in Tables 4-7 for H 2 and in Tables 8-11 for D 2 , for the magnetic field strengths B = 0.1 B 0 and B = 0.2 B 0 . We note that at B = 0.2 B 0 the molecule is meta-stable. As in our previous work on H + 2 , two models have been considered: the approximate model 1, in which off-diagonal terms in v are omitted when setting up the rovibrational matrix, Eq. ( 17), and model 2, in which they are included. The closeness of the two sets of results demonstrate that a simple expansion, with just one vibrational function, yields a good approximation of the final rovibrational wavefunction, at least for the lowest vibrational states. Therefore, in the full expansion of model 2, the coefficients c v,L allow easy identification of the vibrational quantum number of each computed eigenstate. Five digits after the decimal point are given, which reflect the relative accuracy of the rovibrational calculations. The underlying potential curves are accurate to about four digits in absolute terms.

All states are located above the rotational barrier, which is at E barrier = -1.16972 E h for B = 0.1 B 0 and E barrier = -1.15713 E h for B = 0.2 B 0 , respectively, and hence, L, which is an exact quantum number in the field-free case, can still be considered a "good" quantum number. It is interesting to analyse the orientation with respect to the magnetic field axis of the lowest rovibrational state. The lowest state of H 2 , at the field strength of B = 0.1 B 0 , is located 0.0101 E h , or 2222 cm -1 , above the barrier. Yet only one of the basis functions of the expansion in Eq. ( 16) contributes effectively to its eigenvector, with coefficient c 0,1 = 0.997. The eigenfunction of the lowest state is thus

Ψ(R, θ, φ) ∼ [ξ v (R; θ = 0)/R] Y 0 1 (θ, φ) ∼ [ξ v (R; θ = 0)/R]
cos θ, which shows that molecule essentially vibrates in the direction of the magnetic field.

In general, within each L-layer, the rotational energy of a vibrational state increases with |M |. Figures 5 and7 show some exceptions for the states v = 1, 3 of H 2 and v = 0, 2 of D 2 , where the M = 0 state corresponding to L = 2 is above |M | = 1. A similar effect has been observed in the case of H + 2 and D + 2 . It is due to strong coupling of the L = 0 and L = 2 basis functions, a kind of Fermi resonance of the zero-order states with well-defined L. The effect scales as B 2 and is not visible for the lower field strength, B = 0.1 B 0 . No strong effect can be seen for the states v = 0, 2 of H 2 and v = 1, 3 of D 2 , which have L = 1, 3 . . ., where the L = 1 and L = 3 layers are sufficiently separated in energy.

Table 4: Rotational energy levels of H 2 in presence of a uniform magnetic field B for the vibrational state v = 0. The pure vibrational state (L = 0 in the field-free case) is forbidden but shown here nevertheless as it corresponds to the origin of the rotational band. In the simple model 1, terms off-diagonal in v are neglected. In model 2, the full matrix 17 is diagonalized. Values in parentheses are from Ref. [START_REF] Olivares-Pilón | H 2 + , HeH and H 2 : Approximating potential curves, calculating rovibrational states[END_REF].

L

Energy/E h M π Energy/E h B = 0.0

B = 0.1 B = 0.2 model 1 model 2 model 1 model 2 -5 1 -1.15140 -1.15156 -1.13842 -1.13852 5 1 -1.15167 -1.15183 -1.13896 -1.13907 -4 -1 -1.15155 -1.15171 -1.13886 -1.13898 4 -1 -1.15177 -1.15193 -1.13929 -1.13941 -3 1 -1.15167 -1.15184 -1.13919 -1.13932 L = 5 -1.15627 3 1 -1.15184 -1.15200 -1.13952 -1.13965 (-1.15660) -2 -1 -1.15177 -1.15193 -1.13944 -1.

Conclusions

We have investigated the problem of the hydrogen molecule vibrating and rotating in the presence of an external magnetic field for the field strengths of B = 0.1, 0.15, 0.175 a.u. and B = 0.2 a.u. (4.7×10 4 T). It was shown that for B > B cr = 0.178 a.u. H 2 exists in the form of two isolated hydrogen atoms with anti-parallel electron spins to the magnetic field direction. For magnetic fields larger than 12 a.u. the molecule gets bound in parallel configuration with 3 Π u as the ground state, see e.g. [START_REF] Detmer | Hydrogen molecule in magnetic fields: The ground states of the Σmanifold of the parallel configuration[END_REF][START_REF] Detmer | Hydrogen molecule in a magnetic field: The lowest states of the Π-manifold and the global ground state of the parallel configuration[END_REF].

Highly accurate variational calculations, based on a few-parameter physically adequate trial function, are carried out for inclined configurations, where the molecular axis forms an angle θ with respect to the direction of a uniform constant magnetic field. We calculated diamagnetic and paramagnetic susceptibilities (for θ = 45 • for the first time), they closely describe experimental data and agree very well with other calculations or are superior. The two-dimensional potential energy surfaces were built for magnetic fields for B = 0.1 and 0.2 a.u. The parallel orientation of the H 2 molecule with respect to the magnetic field is the most stable one even though the molecule becomes metastable for B = 0.2 a.u. being in domain B > B cr . This holds true also if the vibrational zero-point energy is taken into account. Though the rovibrational ground state is located well above the barrier to perpendicular orientation, the vibrating molecule remains in its parallel orientation. The lowest rovibrational states have then been calculated for the first time. Their energy values are reported for the four lowest vibrational states and rotational excitation up to M = 5, for both the H 2 and D 2 isotopologues.

Figure 1 :

 1 Figure 1: Total energy E of H 2 , 1 g state vs inclination θ for B = 0 and B = 0.1, 0.175, and 0.2 a.u. The solid lines correspond to the hindered rotator model E(B, θ) = E(B, 0 • ) + A sin 2 (θ), where A = (E(B, 90 • ) -E(B, 0 • )), see Eq. (15).

Figure 2 : 2 = c 5 b 6 .

 226 Figure2: Equilibrium distance R eq of H 2 in the 1 g state vs inclination angle θ for B = 0, 0.1, 0.175, 0.2 a.u. The solid lines correspond to R eq (B, θ) = R eq (B, 0 • ) + C sin 2 (θ), where C = R eq (B, 90 • ) -R eq (B, 0 • ) (c.f. Eq[START_REF] Medel Cobaxin | H 2 + in a weak magnetic field[END_REF]).

Figure 3 :

 3 Figure 3: Potential Energy curves of the 1 g ground state of the H 2 molecule for B = 0, 0.1, 0.2 a.u. and θ = 0, 45 • , 90 • . The insets show amplified energy curves for B = 0.1 a.u. (bottom), and B = 0.2 a.u. (top) around the equilibrium distance.For each given magnetic field, the lowest energy curve around the minimum always corresponds to the parallel configuration θ = 0 • , while the highest energy curve corresponds to the perpendicular configuration θ = 90 • . For B = 0.2 a.u. the minimal energy of H 2 corresponds to the asymptotic energy of the repulsive 3 u triplet state (the energy of two hydrogen atoms infinitely separated with spins antiparallel to the magnetic field direction) and lies below the minimum for the 1 g state.

  1.15655 -1.15657 -1.14393 -1.14393 -2 -1 -1.15659 -1.15662 -1.14422 -1.14424 L = 3 -1.16099 2 -1 -1.15670 -1.15672 -1.14444 -1.14445 (-1.16130) -1 1 -1.15672 -1.15674 -1.14450 -1.14452 1 1 -1.15677 -1.15679 -1.14461 -1

LL

  Energy/E h M π Energy/E h B = 0.0 B = 0.1 B = 0.2 model 1 model 2 model 1 model 2 -5 1 -1.11475 -1.11491 -1.10065 -1.10073 5 1 -1.11502 -1.11518 -1.10119 -1.10127 -4 -1 -1.11493 -1.11509 -1.10123 -1.10133 4 -1 -1.11515 -1.11531 -1.10167 -1.10176 -3 1 -1.11508 -1.11524 -1.10165 -11.11963 -1.11966 -1.10645 -1.10646 (-1.12475) -1 1 -1.11966 -1.11969 -1.10649 -1.10651 1 1 -1.11972 -1.11974 -1.10660 -1Energy/E h M π Energy/E h B = 0.0 B = 0.1 B = 0.2 model 1 model 2 model 1 model 2 -4 1 -1.10041 -1.10047 -1.08581 -1.08583 4 1 -1.10063 -1.10069 -1.08624 -1.08626 -3 -1 -1.10063 -1.10069 -1.08659 -1.08662 3 -1 -1.10079 -1.10085 -1.08692 -1

Figure 4 :

 4 Figure 4: Rotational structure, up to L = 5, of the four lowest vibrational states of H 2 in the presence of an external magnetic field B = 0.1 B 0 .

Figure 5 :

 5 Figure 5: Rotational structure, up to L = 5, of the four lowest vibrational states of H 2 in the presence of an external magnetic field of B = 0.2 B 0 .

Figure 6 :

 6 Figure 6: Rotational structure, up to L = 5, of the four lowest vibrational states of D 2 in the presence of an external magnetic field B = 0.1 B 0 .

Figure 7 :

 7 Figure 7: Rotational structure, up to L = 5, of the four lowest vibrational states of D 2 in the presence of an external magnetic field of B = 0.2 B 0 .

Table 1 :

 1 Total electronic energy and equilibrium distance of H 2 in the state 1 g vs magnetic field B and inclination θ based on trial function[START_REF] Goldreich | Pulsar electrodynamics[END_REF], see text. Energies E and equilibrium distances R eq rounded to 5th and 3rd d.d., respectively. * For B = 0.2 a.u. the 1 g state is no longer the ground state. Results marked † are from Ref.[START_REF] Detmer | Hydrogen molecule in magnetic fields: The ground states of the Σmanifold of the parallel configuration[END_REF], those marked ‡ from[START_REF] Sims | High precision variational calculations for the Born-Oppenheimer energies of the ground state of the hydrogen molecule[END_REF] (rounded). The binding energy E bind ≡ 2E(H) -E(H 2 ) with respect to dissociation to H + H is shown in the last column, where the energies for the H atom in ground state are taken from[START_REF] Kravchenko | Exact solution for a hydrogen atom in a magnetic field of arbitrary strength[END_REF].

	B (a.u.) θ (degrees)	E (a.u.)	R eq (a.u.) E bind (a.u.)
	0.0	-	-1.17420	1.40	0.17420
			-1.174476 ‡	1.40
	0.1	0	-1.17047	1.397	0.17542
			-1.16965 †	1.39 †
		45	-1.17014	1.396	0.17508
		90	-1.16983	1.394	0.17477
	0.175	0	-1.16282	1.390	0.17768
		45	-1.16187	1.387	0.17673
		90	-1.16107	1.384	0.17592
	0.2 *	0	-1.15941	1.385	0.17864
			-1.15877 †	1.39 †
		45	-1.15816	1.382	0.17740
		90	-1.15713	1.379	0.17636

Table 2 :

 2 Fitted parameters (rounded to 5 d.d.) in the Padé approximant[START_REF] Turbiner | One-electron molecular systems in a strong magnetic field[END_REF] for the H 2 potential energy curves E vs R for B = 0, 0.1, 0.2 a.u., see Fig.3.

	B(a.u.) θ	a 1	b 2	b 3	b 4	b 5	b 6
	0	0 • -1.28814 1.09317 -0.94343 0.79009 -0.28161 0.04758
	0.1	0 • -1.29188 0.97020 -0.64470 0.51708 -0.17165 0.03136
		45 • -1.29133 1.16648 -1.12137 0.94652 -0.34039 0.05644
		90 • -1.29330 1.30768 -1.38159 1.13784 -0.40315 0.06395
	0.2	0 • -1.30393 1.04298 -0.80653 0.65880 -0.22496 0.03886
		45 • -1.30284 1.09918 -0.99912 0.90103 -0.35196 0.06408
		90 • -1.30207 1.13486 -1.06132 0.93558 -0.35573 0.06400
	ficients in the operator

Table 3 :

 3 Diamagnetic χ d , paramagnetic χ p and total χ susceptibilities of H 2 in the state 1 g for different inclinations θ at R = R eq . The paramagnetic susceptibility χ p obtained as the difference χ p = χ -χ d (see text) is included for completeness. The expectation values of the squared components of the position vector of each electron x 2

	1,2 , y 2 1,2 and z 2 1,2

  ], b from[START_REF] Kolos | Potential-energy curves for the X 1 Σ + g , b 3 Σ + u , and C 1 Π u states of the hydrogen molecule[END_REF], c from[START_REF] Rychlewski | Electric and magnetic properties for the ground and excited states of molecular hydrogen[END_REF].

	θ	x 2 1,2	y 2 1,2	z 2 1,2	χ d	χ p	χ
	0 • 0.76465 0.76465	1.00929	-0.7647	0.0	-0.7647
		0.7608 a 0.7608 a	0.9730 a -0.7608 a	0.0 a
		0.76169 b 0.76169 b 1.02297 b -0.7617 b	0.0 b
					-0.766 exp
	45 • 0.88697 0.76465	0.88697	-0.8258	0.0240	-0.8046
	90 • 1.00929 0.76465	0.76465	-0.8870	0.0258	-0.8612
		0.9730 a 0.7608 a	0.7608 a	-0.8669 a
		1.02297 b 0.76169 b 0.76169 b -0.8923 b
						(0.027 -0.082) c
					-0.913 exp	0.022 exp

Table 5 :

 5 Rotational energy levels of H 2 in presence of a uniform magnetic field B for the vibrational state v = 1. See Caption of Table4for explications.

	L	Energy/E h M π	Energy/E h
		B = 0.0	B = 0.1	B = 0.2
			model 1 model 2 model 1 model 2
			-4 1 -1.13514 -1.13521 -1.12176 -1.12179
			4 1 -1.13536 -1.13543 -1.12220 -1.12222
			-3 -1 -1.13534 -1.13541 -1.12235 -1.12239
			3 -1 -1.13551 -1.13558 -1.12268 -1.12272
	L = 4	-1.14010	-2 1 -1.13549 -1.13556 -1.12272 -1.12277
		(-1.14050)	2 1 -1.13560 -1.13567 -1.12294 -1.12299
			-1 -1 -1.13558 -1.13566 -1.12297 -1.12302
			1 -1 -1.13564 -1.13571 -1.12308 -1.12313
			0 1 -1.13563 -1.13570 -1.12308 -1.12313
			-2 1 -1.13888 -1.13888 -1.12575 -1.12575
			2 1 -1.13898 -1.13899 -1.12597 -1.12597
	L = 2	-1.14364	-1 -1 -1.13918 -1.13918 -1.12666 -1.12666
		(-1.14405)	1 -1 -1.13923 -1.13924 -1.12677 -1.12677
			0 1 -1.13924 -1.13924 -1.12651 -1.12652
	L = 0	-1.14517	0 1 -1.14071 -1.14071 -1.12840 -1.12840
		(-1.14555)

Table 6 :

 6 Rotational energy levels of H 2 in presence of a uniform magnetic field B for the vibrational state v = 2. The pure vibrational state (L = 0 in the field-free case) is forbidden but shown here nevertheless as it corresponds to the origin of the rotational band. See Caption of Table4for explications.

Table 8 :

 8 Rotational energy levels of D 2 in presence of a uniform magnetic field B for the vibrational state v = 0. See Caption of Table4for explications.

	L	Energy/E h M π	Energy/E h
		B = 0.0	B = 0.1	B = 0.2
			model 1 model 2 model 1 model 2
			-4 1 -1.15985 -1.15987 -1.14715 -1.14715
			4 1 -1.15996 -1.15998 -1.14737 -1.14737
			-3 -1 -1.16001 -1.16003 -1.14764 -1.14765
			3 -1 -1.16009 -1.16011 -1.14780 -1.14781
	L = 4 -1.158877	-2 1 -1.16012 -1.16013 -1.14790 -1.14791
			2 1 -1.16017 -1.16019 -1.14800 -1.14802
			-1 -1 -1.16019 -1.16020 -1.14808 -1.14809
			1 -1 -1.16021 -1.16023 -1.14813 -1.14815
			0 1 -1.16022 -1.16024 -1.14814 -1.14816
			-2 1 -1.16185 -1.16185 -1.14936 -1.14936
			2 1 -1.16190 -1.16190 -1.14947 -1.14947
	L = 2 -1.162594	-1 -1 -1.16209 -1.16209 -1.15012 -1.15012
			1 -1 -1.16212 -1.16212 -1.15017 -1.15017
			0 1 -1.16211 -1.16211 -1.14981 -1.14982
	L = 0 -1.164212	0 1 -1.16291 -1.16291 -1.15118 -1.15119

Table 9 :

 9 Rotational energy levels of D 2 in presence of a uniform magnetic field B for the vibrational state v = 1. The pure vibrational state (L = 0 in the field-free case) is forbidden but shown here nevertheless as it corresponds to the origin of the rotational band. See Caption of Table4for explications.

	L	Energy/E h M π	Energy/E h
		B = 0.0	B = 0.1	B = 0.2
			model 1 model 2 model 1 model 2
			-5 1 -1.14472 -1.14475 -1.13162 -1.13163
			5 1 -1.14485 -1.14489 -1.13189 -1.13190
			-4 -1 -1.14487 -1.14491 -1.13208 -1.13210
			4 -1 -1.14498 -1.14501 -1.13230 -1.13232
			-3 1 -1.14499 -1.14502 -1.13238 -1.13241
	L = 5 -1.137626	3 1 -1.14507 -1.14511 -1.13255 -1.13257
			-2 -1 -1.14507 -1.14511 -1.13260 -1.13263
			2 -1 -1.14513 -1.14516 -1.13271 -1.13274
			-1 1 -1.14513 -1.14517 -1.13274 -1.13277
			1 1 -1.14515 -1.14519 -1.13279 -1.13282
			0 -1 -1.14515 -1.14519 -1.13280 -1.13283
			-3 1 -1.14716 -1.14716 -1.13423 -1.13423
			3 1 -1.14724 -1.14725 -1.13439 -1.13439
			-2 -1 -1.14738 -1.14738 -1.13490 -1.13490
	L = 3 -1.142112	2 -1 -1.14743 -1.14744 -1.13501 -1.13501
			-1 1 -1.14748 -1.14749 -1.13496 -1.13497
			1 1 -1.14751 -1.14752 -1.13502 -1.13502
			0 -1 -1.14753 -1.14754 -1.13521 -1.13522
			-1 1 -1.14862 -1.14862 -1.13616 -1.13616
	L = 1 -1.144658	1 1 -1.14864 -1.14864 -1.13621 -1.13622
			0 -1 -1.14899 -1.14899 -1.13715 -1.13715
	L = 0 -1.145172	0 -1 -1.14908 -1.14908 -1.13716 -1.13716

Table 10 :

 10 Rotational energy levels of D 2 in presence of a uniform magnetic field B for the vibrational state v = 2. See Caption of Table4for explications.

	L	Energy/E h M π	Energy/E h
		B = 0.0	B = 0.1	B = 0.2
			model 1 model 2 model 1 model 2
			-4 1 -1.13286 -1.13288 -1.11944 -1.11944
			4 1 -1.13297 -1.13299 -1.11965 -1.11966
			-3 -1 -1.13306 -1.13307 -1.12005 -1.12006
			3 -1 -1.13314 -1.13315 -1.12022 -1.12022
	L = 4 -1.122394	-2 1 -1.13318 -1.13320 -1.12030 -1.12032
			2 1 -1.13324 -1.13325 -1.12041 -1.12042
			-1 -1 -1.13326 -1.13328 -1.12054 -1.12055
			1 -1 -1.13329 -1.13331 -1.12059 -1.12060
			0 1 -1.13330 -1.13332 -1.12058 -1.12060
			-2 1 -1.13474 -1.13474 -1.12161 -1.12162
			2 1 -1.13480 -1.13480 -1.12172 -1.12173
	L = 2 -1.125757	-1 -1 -1.13504 -1.13504 -1.12254 -1.12254
			1 -1 -1.13507 -1.13507 -1.12259 -1.12259
			0 1 -1.13503 -1.13503 -1.12204 -1.12206
	L = 0 -1.127217	0 1 -1.13582 -1.13582 -1.12367 -1.12367

Table 11 :

 11 Rotational energy levels of D 2 in presence of a uniform magnetic field B for the vibrational state v = 3. The pure vibrational state (L = 0 in the field-free case) is forbidden but shown here nevertheless as it corresponds to the origin of the rotational band. See Caption of Table4for explications.

	L	Energy/E h M π	Energy/E h
		B = 0.0	B = 0.1	B = 0.2
			model 1 model 2 model 1 model 2
			-5 1 -1.11895 -1.11898 -1.10501 -1.10502
			5 1 -1.11908 -1.11912 -1.10529 -1.10529
			-4 -1 -1.11912 -1.11915 -1.10559 -1.10560
			4 -1 -1.11923 -1.11926 -1.10581 -1.10582
			-3 1 -1.11924 -1.11928 -1.10593 -1.10595
	L = 5 -1.102915	3 1 -1.11933 -1.11936 -1.10609 -1.10611
			-2 -1 -1.11934 -1.11938 -1.10619 -1.10621
			2 -1 -1.11939 -1.11943 -1.10630 -1.10632
			-1 1 -1.11940 -1.11944 -1.10634 -1.10636
			1 1 -1.11942 -1.11946 -1.10639 -1.10642
			0 -1 -1.11943 -1.11947 -1.10641 -1.10644
			-3 1 -1.12123 -1.12124 -1.10752 -1.10752
			3 1 -1.12131 -1.12132 -1.10768 -1.10769
			-2 -1 -1.12148 -1.12148 -1.10836 -1.10836
	L = 3 -1.107489	2 -1 -1.12153 -1.12154 -1.10847 -1.10847
			-1 1 -1.12158 -1.12159 -1.10828 -1.10829
			1 1 -1.12161 -1.12162 -1.10833 -1.10835
			0 -1 -1.12164 -1.12164 -1.10867 -1.10868
			-1 1 -1.12262 -1.12262 -1.10960 -1.10960
	L = 1 -1.109872	1 1 -1.12264 -1.12264 -1.10966 -1.10966
			0 -1 -1.12303 -1.12303 -1.11073 -1.11073
	L = 0 -1.110336	0 -1 -1.12310 -1.12310 -1.10896 -1.10898

Acknowledgements

The authors thank the high-performance computer centres ROMEO of the University of Reims Champagne-Ardenne, CRIANN of the Region of Normandy, France and cluster KAREN (ICN-UNAM, Mexico) for generous allowance of super-computer time. The research by J.C.L.V., D.J.N., A.V.T. is partially supported by CONACyT grant A1-S-17364 and DGAPA grant IN108815 (Mexico). This work was also supported by the Programme National de Planétologie (PNP) of CNRS/INSU, co-funded by CNES. Two of the authors A.A. and A.V.T. have the honor and the privilege to know closely Vladimir Tyuterev to whom this paper is dedicated.