Ionization and dissociation equilibrium in strongly-magnetized helium atmosphere

Recent observations and theoretical investigations of neutron stars indicate that their atmospheres consist not of hydrogen or iron but possibly other elements such as helium. We calculate the ionization and dissociation equilibrium of helium in the conditions found in the atmospheres of magnetized neutron stars. For the first time this investigation includes the internal degrees of freedom of the helium molecule. We found that at the temperatures and densities of neutron star atmospheres the rotovibrational excitations of helium molecules are populated. Including these excitations increases the expected abundance of molecules by up to two orders of magnitude relative to calculations that ignore the internal states of the molecule; therefore, if the atmospheres of neutron stars indeed consist of helium, helium molecules and possibly polymers will make the bulk of the atmosphere and leave signatures on the observed spectra from neutron stars. We applied our calculation to nearby radio-quiet neutron stars with B_{dipole}~10^{13}-10^{14} G. If helium comprises their atmospheres, our study indicates that isolated neutron stars with T_{BB}~10^6 K such as RXJ0720.4-3125 and RXJ1605.3+3249 will have He^+ ions predominantly, while isolated neutron stars with lower temperature (T_{BB}~5x10^5 K) such as RXJ1856.5-3754 and RXJ0420.0-5022 will have some fraction of helium molecules. If helium molecules are abundant, their spectroscopic signatures may be detected in the optical, UV and X-ray band. (Abridged)


INTRODUCTION
Hydrogen has been considered as the surface composition of isolated neutron stars (INSs) because gravitational stratification forces the lightest element to the top of the atmosphere (Alcock & Illarionov 1980). Only a tiny amount of material is required to constitute an optically thick layer on the surface (Romani 1987). However, recent studies of  and  have shown that the NS surface may be composed of helium or heavier elements since hydrogen may be quickly depleted by diffusive nuclear burning. Observationally, helium and heavier element atmospheres have been proposed for interpreting the spectral features observed in several INS partially because the existing hydrogen atmosphere models do not reproduce the observed spectra (Sanwal et al. 2002;van Kerkwijk & Kaplan 2006). However, atomic and molecular data in the strong magnetic field regime are scarce for non-⋆ e-mail: kaya@cita.utoronto.ca † e-mail: heyl@physics.ubc.ca hydrogenic elements. Accurate atomic and molecular data are available mostly for the He + ion (Pavlov & Bezchastnov 2005), the helium atom (Neuhauser et al. 1987;Demeur et al. 1994), He 3+ 2 (Turbiner & López Vieyra 2004) and He 2+ 2 (Turbiner & Guevara 2006). Helium molecular binding energies have been crudely calculated by density functional theory (Medin & Lai 2006a,b) (hereafter ML06). Unlike hydrogen atmospheres (Lai & Salpeter 1997), the ionization and dissociation balance in strongly-magnetized helium atmosphere has not been investigated yet.
In this paper, we extend our Hartree-Fock type calculation  to helium molecules in the Born-Oppenheimer approximation. For molecular ions that exist in strong magnetic fields (B = 10 12 -10 14 G), we achieved < ∼ 1% and < ∼ 10% agreement in binding energies and vibrational energies in comparison with other more accurate studies mainly on hydrogen molecules. Including numerous electronic, vibrational and rotational states, we studied ionization and dissociation equilibrium in helium atmospheres at B = 10 12 -10 14 G. We also applied our calculations to several INSs which may have helium atmospheres on their surfaces. c 2006 RAS

MOLECULAR BINDING AND VIBRATIONAL ENERGY
At first we adopt the Born-Oppenheimer approximation and neglected any effects associated with motion of atoms and molecules in a magnetic field. Later, we will discuss rotovibrational states ( §3) and how the finite nuclear mass modifies results ( §5). In the Landau regime (β Z > 1 where β Z ≡ B/(4.7 × 10 9 Z 2 G) and Z is the atomic number), bound electrons in an atom and molecule are well specified by two quantum numbers (m, ν). m is the absolute value of a magnetic quantum number (which is negative to lower the total energy in strong magnetic fields) and ν is a longitudinal quantum number along the field line. We consider only tightly-bound states with ν = 0. Electronic excited states with ν > 0 have small binding energies, therefore their population in the atmosphere is tiny due to small Boltzmann factors and pressure ionization. Hereafter we denote atomic and molecular energy states as (m 1 , m 2 , m 3 , m 4 , ...).
We computed molecular binding energies with a simple modification to our Hartree-Fock type calculation for atoms . We replaced the nuclear Coulomb term V m (z) in the Schrödinger equation by (Lai et al. 1992) The function Φ m is the ground state Landau wavefunction, and a is the separation between two nuclei. We added Z 2 e 2 a as the Coulomb repulsion energy between the two nuclei. We computed binding energies with a grid size ∆a ∼ 0. 1 [a.u.] up to a ∼ 1 [a.u.] and ∆a ∼ 0.01 [a.u.] near the energy minimum. Figure 1 shows the binding energy curve of He 2 at B = 10 12 G fitted with the Morse function defined as (Morse 1929) where E m (> 0) is the molecular binding energy for an electronic state m (in this paper this usually denotes the magnetic quantum number of the outermost electron), a 0 is the separation between two nuclei at the minimum energy. We defined two different dissociation energies: The calculated binding energy values are not smooth near the energy minimum ( Figure 1). This is due to our numerical errors. Binding energy does not change by more than 0.1% for ∆a = 0.01 [a.u.] near the energy minimum. We determined E m from the fitting procedure using the function given by equation (3) since we found that it provides more accurate results than the minimum energies from our grid calculation. However, in most cases, E m from our grid calculation and the fitted E m do not differ by more than 1%. We computed D m andD m using the atomic data we calculated numerically.
Our results for H + 2 and H 2 are in good agreement  with Turbiner & López Vieyra (2003) (hereafter TL03) and Lai & Salpeter (1996) (LS96) with less than 1% deviation in total binding energy (Tab. 1). TL03 performed highly accurate variational studies mainly on one-electron molecular systems (e.g., H + 2 , H 2+ 3 , He 3+ 2 ). LS96 studied hydrogen molecular structure similarly by a Hartree-Fock calculation in the adiabatic approximation. While our calculation takes into account higher Landau levels using perturbation theory, the difference in binding energies by including higher Landau levels is tiny for helium atoms and molecules at B > 10 12 G.

Electronic excitation
The electronic excited states in molecules occupy higher (m, ν) states than those in atoms. Since excitation energies from the ground state to (0, 1, 2, m) states with m > 3 are small, there may be numerous tightly-bound electronic excited states until they dissociate into atoms and ions at large m. We did not consider excited states with ν > 0 because their binding energies are small therefore they are likely to be dissociated. We calculated binding energies for the (0, 1, 2, m) state up to m = 9 and estimated binding energies for higher m states using the well-known m dependence of the energy spacing ∆E m ∼ ln 2m+3 2m+1 (Lai et al. 1992). We found that the difference between the exact solutions and those from the scaling law is tiny at m & 9. Figure  2 shows the He 2 binding energy of (0, 1, 2, m) states at B = 10 12 G. Note that the excited states of He 2 , (0, 1, 2, m) with m > 5, are unbound with respect to two atoms in the ground state at B = 10 12 G.

ROTOVIBRATIONAL EXCITATION
We consider molecular excitation levels associated with vibrational and rotational motion of molecules in a magnetic field. In contrast to the field-free case, the strong magnetic field induces molecular oscillations with respect to the field line similar to a twodimensional harmonic oscillator. Accordingly, there are three types of molecular motion; vibration along and transverse to the magnetic field and rotation around the magnetic field. Hereafter we briefly describe energy levels of rotovibrational states.
Strictly speaking, the aligned and transverse vibrations are coupled (Khersonskii 1984;Lai & Salpeter 1996). However, using perturbation theory, Khersonskii (1985) has shown that the coupling energy is tiny (less than 1%) compared to the total binding energy. Neglecting the coupling, the rotovibrational energy levels are well approximated by (4) ǫ V is the aligned vibrational energy given by The integer V(> 0) is the quantum number for the aligned vibration and~ω is the aligned vibrational energy quanta (Morse 1929;Khersonskii 1985).
On the other hand, the transverse rotovibration energy (ǫ NΛ ) consists of transverse vibration and rotation around the magnetic field axis and it is given as (Khersonskii 1985) The integer N(> 0) is the quantum number for the transverse vibration, while the integer Λ is the projection of angular momentum in the B-field direction (Khersonskii 1985).~Ω B is the nuclear cyclotron energy (= Z~eB/(Am p c) = 6.3(Z/A)B 12 [eV] where A is the atomic mass) and ω t = 4ω 2 ⊥ + Ω 2 B . ω ⊥ is the transverse vibrational energy quanta. The nuclear cyclotron energy term takes into account the magnetic restoring force on the nuclei (Lai & Salpeter 1996). In the following subsections, we calculate vibrational energy quanta in the Born-Oppenheimer approximation.

Aligned vibrational excitation
In the Born-Oppenheimer approximation, the motion of two nuclei along the magnetic field is governed by the binding energy curves determined in section 2. It is therefore straightforward to calculate aligned vibrational energy quanta using the results from section 2. We fit the Morse function given by equation (3) to molecular binding energies as a function of the nuclear separation a. Once ζ is determined, the aligned vibrational energy quanta is given by (LS96) where µ is the reduced mass of the two nuclei in units of the electron mass (918 for H 2 and 3675 for He 2 ). At large a, another electron configuration is mixed with tightly-bound states (Lai et al. 1992). Since configuration interaction is neglected in our calculation, our ζ values (therefore~ω ) are overestimated by 10-30% in comparison with LS96 (table 3). We found that ζ is nearly identical for different (tightly-bound) electronic excited states. Therefore we computed~ω for electronic excited states with large m (for which we did not perform grid calculation) using ζ from the lower excited states. In most cases, our aligned vibrational energy quanta agree with other more accurate results within ∼ 10% (table 3). Table 4 compares our results for He 2+ 2 with Turbiner & Guevara (2006) at various magnetic fields. Table 5 shows aligned vibrational energy quanta for helium molecular ions along with some results on He 3+ 2 from Turbiner & López Vieyra (2004).
It is apparent that the discrepancies with other calculations are larger for helium molecules than for hydrogen molecules. There are two effects both of which reduce the accuracy of our results particularly for highly-ionized molecules at low magnetic fields. First, we did not take into account configuration interaction, while the work Notes: We multiplied the results of TL03 by a factor of two to correct the discrepancy due to different definitions of aligned vibrational energy quanta. by Turbiner et al. employed the full 2-dimensional variation energy calculation. The degree of configuration interaction is larger at small β Z since the increasing effects of the nuclear Coulomb field mix different electron configurations. Second, highly-ionized molecular ions are either unbound or weakly-bound at low magnetic fields, therefore the numerical errors in binding energies significantly affect determination of ζ since the binding energy curve is shallow. Nevertheless, the accuracy of vibrational energies of highly-ionized molecules such as He 3+ 2 is irrelevant for dissociation balance since they are unbound or their abundance is negligible (section 6). For the abundant molecular ions in 10 12 -10 14 G (e.g., He + 2 and He 2 ), we expect the accuracy of aligned vibrational energy is < ∼ 10%.

Transverse vibrational excitation
We calculated the energy curve as a function of transverse position of nuclei R following Ansatz A described in section IIIB of LS96. We fixed a to the equilibrium separation a 0 and supposed that the two nuclei are located at (±R/2, ±a 0 /2). As LS96 pointed out, this method is appropriate for small R (.ρ whereρ is the cyclotron radius) and gives only an upper limit to transverse vibrational en-   (2004) are 11 eV and 51 eV at B 12 = 1 and 10 respectively. ergy quanta~ω ⊥ . We replaced the nuclear Coulomb term in the Schrödinger equation V(z) by We added Z 2 e 2 (a 2 0 +R 2 ) 1/2 as the Coulomb repulsion energy between two nuclei. Once we calculated the molecular binding energy at different R grid points, we fit a parabolic form 1 2 µω 2 ⊥ R 2 to binding energies at R .ρ. We found that~ω ⊥ is nearly identical for different electronic excitation levels. Therefore, we adopted~ω ⊥ of the ground state for electronic excited states. Our transverse vibrational energy quanta agree with those of LS96 and TL03 within 10% (table 6). Table 7 shows the results for helium molecular ions in comparison with those for He 3+ 2 from Turbiner & López Vieyra (2004). Compared to aligned vibrational energy, our results for transverse vibrational energy are in better agreement with other more accurate results. The better agreement is well-understood because the transverse potential well is deeper than in the aligned direction in magnetic field, therefore our results are less subject to the numerical errors as discussed in section 3.1.

Perturbative approach
It is also possible to estimate the transverse vibrational energy perturbatively, by calculating the lowest order perturbation to the energy of the molecule induced by tilting it. We assume that the energy of the tilted molecule is almost the same as the energy required to displacing the electron cloud relative to the molecule by an electric field (E) : The first-order change to the wavefunction is where ǫ (0) κ is the unperturbed energy of the state κ and κ ′ denotes the other states of the system.
The expectation value of x for this situation is where r is half the distance between the nuclei. We can solve for the value of eE that we need to apply to give a particular displacement and substitute it into the expression for the energy of the state to second order The only states that contribute to the sum have m ′ = m ± 1 for the electronic states of the molecule (the rotational states do not count because that is what we are examining).
States that have been excited along the field (ν > 0) or by increasing the Landau number do not have much overlap in the integral in the numerator and also a large energy difference in the denominator.
The frequency for low amplitude oscillations is given by where I is the moment of inertia of the molecule, 2Mr 2 (M is the nuclear mass). The size of the molecule has cancelled out. For a single electron system in the ground state if we assume that the bulk of the contribution to the sum in the equation is given by the m = 1, ν = 0 state we have where the final term is the overlap of the longitudinal wavefunction of the two states. For a multielectron system, evaluating equation (13) and (14) is somewhat more complicated. For clarity of nomenclature we shall write the wavefunctions of the various electronic states of the molecule as K ′ . The symbol K denotes the state that we are focused upon. The change in the energy of the system due to an applied electric field is where N is the number of electrons and where j counts over the electrons in the molecule. We have assumed that the multielectron wavefunctions are normalized such that K|K = 1 The first-order change to the wavefunction is where ǫ (0) K ′ is the unperturbed energy of the state K ′ . Now let us calculate the expectation value ofx for this situation, where r is half the distance between the nuclei. We can solve for the value of NeE that we need to apply to give a particular displacement and substitute it into the expression for the energy of the state to second order In strongly magnetized atoms or molecules, it is natural to expand the wavefunctions in terms of the ground Landau level using various values of m. We assume that the wavefunction for each electron is written as In this case we have otherwise it vanishes. Combining these results yields an estimate for the frequency of low amplitude oscillations of where the subscript on energy labels the value of m in the outermost shell. The number of electrons appears in this equation because the expectation value ofx is a factor N smaller than the expectation value of x for the single shifted electron.
For the ground state of the molecule we have m + 1 = N yielding a simpler expression, Table 6 and 7 compare the perturbative estimates with the numerical calculations. Within the approximations made the agreement is encouraging.

ROTOVIBRATIONAL SPECTRUM
Given the aligned and transverse vibrational energies, we construct the rotovibrational spectra of helium molecules. From the equations (5) and (6), the molecular system has a finite zero-point energy associated with aligned and transverse vibration when N = Λ = V = 0. Therefore, the actual dissociation energy will be reduced from D m by ǫ 000 (LS96). Figure 3 shows the rotovibrational energy spectrum of He 2 in the ground state at B = 10 12 G and 10 13 G respectively.  with excitation energies that are smaller than or similar to the thermal energy have large statistical weights in the partition function ( §6).

EFFECTS OF FINITE NUCLEAR MASS
The separation of the center-of-mass motion is non-trivial when a magnetic field is present ). The total pseudomomentum ( K) is often used to take into account motional effects in a magnetic field since K is a constant of motion (Lai 2001). In the following subsections, we discuss two effects associated with finite nuclear mass in strong magnetic fields. We denote the binding energy in the assumption of fixed nuclear location (e.g., Born-Oppenheimer approximation) as ǫ (0) κ (< 0) for an electronic state κ. Since we consider states with n = 0 and ν = 0, the relevant quantum numbers are κ = {m i } (i denotes each bound electron in multi-electron atoms and molecules).

Finite nuclear mass correction
The assumption of zero transverse pseudomomentum introduces an additional term s κ~ΩB in the binding energy ;  Wunner et al. 1981). Ω B is the nuclear cyclotron energy and s κ = i m i is the sum of magnetic quantum numbers for a given electronic state κ (e.g., s κ = 4 for He 2 molecule in the ground state). However, the scheme assuming the zero transverse pseudomomentum does not necessarily give the lowest binding energies at B > ∼ B Q . Instead, LS95 and LS96 estimated lower binding energies at B > ∼ B Q using another scheme which relaxed the assumption of the zero transverse pseudomomentum. A more rigorous calculation was performed for He + ion by Bezchastnov et al. (1998) and Pavlov & Bezchastnov (2005). An application of such schemes to multi-electron systems is beyond the scope of this paper. We will discuss the limitation of our models at very high magnetic field in §6.

Motional Stark effects
When an atom or molecule moves across the magnetic field B, a motional Stark electric field E MS = K× B MC is induced in the center-ofmass frame. M is the mass of atom or molecule. The Hamiltonian for the motional Stark field is given by where Ω B0 ≡ eB Mc . For a given pseudomomentum K, the motional Stark field separates the guiding center of the nucleus and that of the electron by Since the motional Stark field breaks the cylindrical symmetry preserved in magnetic field, it is non-trivial to evaluate motional Stark field effects. A non-perturbative (therefore more rigorous) approach has been applied only for one-electron systems (Vincke et al. 1992;Potekhin 1994;Pavlov & Bezchastnov 2005). However, such an approach is quite complicated and timeconsuming especially for multi-electron atoms and molecules. Therefore, following LS95, we considered two limiting cases: (1) R K ≪ρ and (2) R K ≫ρ and determined general formula which can be applied to a wide range of B and K ⊥ . For diatomic molecules, we can apply a nearly identical scheme used for calculating transverse vibrational energy to the both cases.

Decentered states
When R K ≫ρ, it is convenient to utilize the so-called decentered formalism (LS95). We replace the nuclear Coulomb term by and compute binding energies at different R K grid points. For diatomic molecules, we replace the nuclear Coulomb term by The grid calculation for R K is identical to the one for transverse vibrational energy in §3.2 except that the Coulomb repulsion term between two nuclei is Z 2 e 2 a 0 (instead of Z 2 e 2 (a 2 0 +R 2 K ) 1/2 ) since motional Stark field shifts the guiding center of the two nuclei by R K in the transverse direction but the separation between the two nuclei is still a 0 .
LS95 found that the binding energy curves are well fit by the following formula. where A 1 , A 2 and A 3 are the fit parameters. E ⊥ (R K ) = E ⊥ (K ⊥ ) is the transverse energy and ǫ (0) κ is the binding energy in the infinite nuclear mass approximation. Figure 4 shows the binding energy curve of He 2 molecule as a function of R K at B = 10 12 G. At small R K , the fitted function is well matched with the results from the perturbation approach. Although mixing between different m states is ignored, a comparison with Potekhin (1998) and Potekhin (1994) for hydrogen atoms indicates this approach gives better than 30% accuracy over a large range of K ⊥ (Lai & Salpeter 1995). This is adequate for our purpose of investigating ionization and dissociation balance.
Along with the finite nuclear mass term discussed in §5.1, the electronic energy of an atom or molecule moving with transverse pseudomomentum K ⊥ is given by Note that ǫ (0) κ < 0 and both the second and third term decrease the binding energy. Later, we will discuss the validity of this approach at B > ∼ B Q .

IONIZATION AND DISSOCIATION EQUILIBRIUM
We have investigated the ionization and dissociation balance of magnetized helium atmospheres including the following chemical reaction channels. Table 9 and 10 list ionization and dissociation energies of various helium ions and molecular ions in the assumption with fixed nuclear location. We did not take into account the He − ion since its ionization energy is . 20 [eV] and He − is not abundant at all in the temperature range considered here (T & 10 5 K).
• Ionization (1) He + ↔ α + e (2) He ↔ He + + e (3) He 2+ 2 ↔ He 3+ 2 + e (4) He + 2 ↔He 2+ 2 + e (5) He 2 ↔He + 2 + e • Dissociation (6) He 3+ 2 ↔ α + He + (7) He 2+ 2 ↔ α + He (8) He 2+ 2 ↔ He + + He + (9) He + 2 ↔ He + + He (10) He 2 ↔ He + He The Saha-Boltzmann equation for ionization balance is given as, n i n i+1 n e = z i z i+1 z e where z e is the partition function for an electron, η e =~ω B /2kT and ω B = eB/m e c is the electron cyclotron frequency (Gnedin et al. 1974). The quantity z i is the partition function for ionization state i, where E i , Z i and M i are the ground state energy, the charge and the mass of an ion i. η i =~Ω Bi /2kT where Ω Bi = Z i eB/M i c is the ion cyclotron frequency. ǫ i,κ (K ⊥ )(< 0) is the binding energy of a bound state κ given in equation (36). w i,κ is the occupation probability. In general, w i,κ (K ⊥ ) is a function of Z i , Z p (the charge of perturbing ions, usually the effective charge of the plasma), ǫ i,κ (K ⊥ ) and ρ (the plasma density). w depends not only on the type of ion i, the electronic state κ and the transverse pseudomomentum K ⊥ , but also it is different for each bound electron. Obviously, electrons in the outer shells are subject to a stronger electric field from neighboring ions than those in the inner shells. We explicitly computed w i,κ for all the bound electrons using the electronic microfield distribution of Potekhin et al. (2002). The Saha-Boltzmann equation for dissociation balance is given as wherez j is the partition function for molecular ionization state j, where E j , Z j and M j are the ground state energy, the charge and the mass of a molecular ion j.
where Ω B j = Z j eB/M j c is the molecular ion cyclotron frequency. ǫ j,κ (K ⊥ )(< 0) and w j,κ (K ⊥ ) are the binding energy and occupation probability of an electronic bound state κ. ǫ NΛM (> 0) is the excitation energy of a rotovibrational state (N, Λ, V). We took the summation N,Λ,V until ǫ NΛV exceeds the dissociation energy D κ . The set of the Saha-Boltzmann equations along with the condition for the baryon number conservation and charge neutrality are iteratively solved until we reached sufficient convergence in ∆n e /n e (< 10 −6 ) where n e is the density of free electrons. The convergence was achieved rapidly in most cases (less than 10 iterations were required). Figure 5, 6 and 7 show the fraction of helium ions and molecular ions at B = 10 12 , 10 13 and 10 14 G. At B = 10 12 G, He 3+ 2 and He 2+ 2 are not present because they are not bound with respect to their dissociated atoms and ions (table 10). In all the cases, the He 3+ 2 fraction is negligible because helium molecular ions with more electrons have much larger binding energies. At B = 10 14 G, He 2 is not present because even the ground state becomes auto-ionized to He + 2 ion due to the finite nuclear mass effect although He 2 remains bound with respect to two helium atoms. Note that the ionization energy of He 2 at 10 14 G is 643 eV (table 9) while the difference in the energy due to the finite nuclear mass term between He 2 and He + 2 is 945 eV. As discussed in §5.1, our scheme using the zero pseudomomentum becomes invalid at B > ∼ B Q and in reality He 2 will have lower binding energy than He + 2 . Therefore, our results at B = 10 14 G should be taken with caution. Figure 8 shows the temperature dependence of the helium atomic and molecular fractions at different B-field strengths. We fixed the plasma density to the typical density of the X-mode photosphere (∼ 10, 10 2 and 10 3 g/cm 3 for B = 10 12 , 10 13 and 10 14 G (Lai & Salpeter 1997)). The transition from molecules to helium atoms and ions takes place rapidly at T ∼ 3 × 10 5 K and ∼ 6 × 10 5 K for B = 10 12 G and 10 13 G.
In order to illustrate different physical effects, we show the fraction of helium atoms and He 2 molecules as a function of temperature in Figure 9. When motional Stark effects are ignored (dotted line), the molecular fraction is underestimated because helium ions and atoms are subject to a larger motional Stark field due to their smaller masses and binding energies. As expected, molecules become more abundant when rotovibrational states are included. The discrepancy from the results without rotovibrational states (dashed line) increases toward higher temperature since more rotovibrational states have larger statistical weight as their excitation energies are an order of 100 eV at B = 10 13 G (see figure 3). When we take into account only the ground states (dot-dashed line), the results are close to the case without motional Stark effects because the ground states are less affected by motional Stark effects than the excited states.

He 3 molecule and larger helium molecular chains
When He 2 is abundant, larger helium molecules such as He 3 may become abundant. We roughly estimated the fraction of He 3 molecules by neglecting the finite nuclear mass effects and zero- point vibrational energy correction both of which will reduce the dissociation energy. Similarly to He 2 , we computed the binding energy of He 3 molecules by Hartree-Fock calculation. The dissociation energy of He 3 →He 2 +He is 289 and 1458 eV at B = 10 13 and 10 14 G. Our results are close to those of the density functional calculation by ML06 (384 and 1647 [eV] at B = 10 13 and 10 14 G). Note that the density function calculation overestimates binding energies by ∼ 10% compared with more accurate Hartree-Fock calculation (ML06). Assuming the internal degrees of freedom are similar between He 2 and He 3 molecule and neglecting the bending degree of freedom for He 3 , we computed dissociation equilibrium between He 2 and He 3 following Lai & Salpeter (1997). Figure 10 shows contours of the He 3 molecule fraction with respect to He 2 molecule at B = 10 13 G. In the dotted area, the fraction of diatomic helium molecular ions is larger than 10%. It is seen that He 3 molecule will be abundant at T . 3 × 10 5 K at B = 10 13 G. The results are consistent with the estimates of molecular chain formation by Lai (2001) and ML06. ML06 provided binding energy of infinite He molecular chain as E ∞ ≃ 1.25B 0.38 13 keV at B & 10 13 G. From this fitting formula, the cohesive energy of helium molecular chains is given as E co ∼ 0.36, 1.5 and 5.1 keV at B = 10 13 , 10 14 and 10 15 G. When kT . 0.1E co , molecular chains are likely to be formed (Lai 2001); consequently, we expect helium molecular chains to form at T . 3 × 10 5 , 1 × 10 6 and 5 × 10 6 K and B & 10 13 , 10 14 and 10 15 G.

APPLICATION
We investigated the ionization and dissociation balance of helium atmospheres for two classes of INS whose X-ray spectra show absorption features.

Radio-quiet neutron stars
A class of INS called radio-quiet neutron stars (RQNS) is characterized by their X-ray thermal spectra with T . 10 6 K and spin-down dipole B-field strength B & 10 13 G. A single or multiple absorption features have been detected from six radio-quiet NS (Haberl et al. 2003(Haberl et al. , 2004van Kerkwijk et al. 2004). Although the interpretation of these features is still in debate, helium is certainly one of the candidates for the surface composition of RQNS (van Kerkwijk & Kaplan 2006). Timing analysis suggests that a couple of RQNS have dipole magnetic field strengths in the range of B = 10 13 -10 14 G (Kaplan & van Kerkwijk 2005). Therefore, we investigated the ionization/dissociation balance of a helium atmosphere at B = 3 × 10 13 G. Figure 11 shows the contours of He + and He molecular fractions at 3 × 10 13 G. He + is dominant at T ∼ 10 6 K and molecules become largely populated at T . 5 × 10 5 K.
Suppose all the RQNS have helium atmospheres on the surface. RQNS with higher temperatures (T ∼ 10 6 K) such as RXJ0720.4-3125 and RXJ1605.3+3249 will have He + ions predominantly with a small fraction of He atoms. Indeed, several bound-bound transition lines of He + ion have energies  (Pavlov & Bezchastnov 2005). On the other hand, RQNS with lower temperatures (T ∼ 5 × 10 5 K) such as RXJ0420.0-5022 may have some fraction of helium molecules, although the He + ion is still predominant at low densities where absorption lines are likely formed.
However, our study shows that the fraction of He 3+ 2 molecular ions is negligible at any B-field and temperature because more Figure 8. Temperature dependence of helium atomic and molecular fraction at B = 10 12 G (ρ = 10 g/cm 3 ), B = 10 13 G (ρ = 10 2 g/cm 3 ) and B = 10 14 G (ρ = 10 3 g/cm 3 ). Figure 9. Fraction of helium atoms (left) and helium molecules (right) ρ = 10 2 g/cm 3 at B = 10 13 G. The solid lines show the results by taking into account all the physical effects discussed in this paper. The dotted and dashed lines ignored motional Stark effects and rotovibrational states respectively. The dot-dashed lines include only the ground states. neutral molecular ions have significantly larger binding energies. We also investigated ionization/dissociation balance of helium atmospheres at 2 × 10 14 G (figure 12). At B > ∼ B Q , our scheme of treating finite nuclear mass effects becomes progressively inaccurate. As a result, the ionization energy of helium atoms becomes significantly smaller and the He 2 molecule becomes auto-ionized due to the finite nuclear mass correction discussed in §5.1. Therefore, the molecular fraction will be underestimated in this case (the left panel in figure 12). On the other hand, when we ignore the finite nuclear mass effects (the right panel in figure 12), the molecular fraction is likely overestimated. We expect that a realistic ionic and molecular fraction will be somewhere between the two cases. It is premature to conclude whether He + is dominant for the case of 1E1207 since we do not have a self-consistent temperature and density profile. The study of ML06 and Lai (2001) also suggests the critical temperature below which helium molecule chains form is ∼ 2 × 10 6 K at B = 2 × 10 14 G. This is close to the blackbody temperature of 1E1207. Further detailed studies are necessary to Figure 10. Ratio of He 3 with respect to He 2 at B = 10 13 G (solid curve contours). The dotted area shows the (ρ, T ) regime where a fraction of diatomic helium molecular ions is larger than 10%. Figure 11. Ionization and dissociation balance of helium at B = 3 × 10 13 G. The dotted area is where the fraction of helium molecules and molecular ions (including He 2 , He + 2 , He 2+ 2 and He 3+ 2 ) becomes more than 50%. The black area is where the fraction of He + becomes more than 50%. The white region on the right is where the fraction of bare helium ions becomes more than 50%. In the narrow white region betweeen the dotted and black area, atomic helium is somewhat abundant, so neither molecular nor singly ionized helium dominates.
conclude the composition of helium atmospheres at B = 2 × 10 14 G.

DISCUSSION
We have examined the ionization-dissociation balance of the helium atmospheres of strongly magnetized neutron stars. As the observational data on isolated neutron stars has improved over the past decade (Sanwal et al. 2002;van Kerkwijk & Kaplan 2006), both hydrogen and iron atmospheres have lost favour, and atmospheres composed of other elements have been invoked to interpret both the continuum and spectral properties of neutron stars. On the other hand the theoretical investigations of  and  have shown that diffusive nuclear burning may quickly deplete the hydrogen by so the NS surface may be composed of helium or heavier elements.
The results found here indicate several avenues for future work. Although the treatment of finite nuclear mass effects is problematic in the strong field limit, only by including this accurately can we know the physical state of matter in super-critical magnetic fields -solid, liquid or gas. The state can have significant effects of the outgoing radiation from these objects. It is also important to repeat a similar calculation to that presented here for molecular chains, including all of the relevant degrees of freedom of the chains, i.e. including the bending modes which were neglected here. Regardless, the high abundance of molecules under the conditions of observed neutron star atmospheres is bound to spur additional research into the statistical mechanics of highly magnetized molecules and polymers.
At the temperatures and densities of neutron star atmospheres the rotovibrational excitations of helium molecules are populated. Including these excitations increases the expected abundance of molecules by up to two orders of magnitude relative to calculations that ignore the internal states of the molecule. Ionization, dissociation and electric excitation energies of helium molecules are larger than 100 eV at B > ∼ 10 13 G. On the other hand, rotovibrational excitation energies are in the range of 10-100 eV at B = 10 12 -10 14 G.
If helium molecules are abundannt, their spectroscopic signatures may be detected in the optical, UV and X-ray band. If helium comprises the atmospheres of isolated neutron stars, clearly it is crucial to understand the structure of helium molecules and molecular chains in order to interpret the spectra from neutron stars.

ACKNOWLEDGMENTS
We would like to thank the anonymous referee for comments that helped to improve the paper. KM acknowledges useful discussions with George Pavlov and Marten van Kerkwijk. This work was supported in part by a Discovery Grant from NSERC (JSH). This work made use of NASA's Astrophysics Data System. The authors were visitors at the Pacific Institute of Theoretical Physics during the nascent stages of this research.