format: text

text:
# ISO spectroscopy of circumstellar dust in the Herbig Ae systems AB Aur and HD 163296 Based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with the participation of ISAS and NASA and on observations collected at the European Southern Observatory, La Silla, Chile.

## 1 Introduction

Herbig Ae/Be stars are intermediate-mass pre-main sequence stars surrounded by disks of gas and dust which might be the site of on-going planet formation (see Waters & Waelkens 1998 for a recent review). AB Aur (A0Ve) and HD 163296 (A3Ve) belong to the best studied Herbig stars and are sometimes considered prototypical for the entire class. As early as 1933, Merrill & Burwell remarked upon the similarity of both systems, which was confirmed by numerous subsequent authors. Apart from the fact that the stellar mass, effective temperature and age of AB Aur and HD 163296 are nearly identical (van den Ancker et al. 1997, 1998), the similarity between the two systems also extends to their circumstellar environment: both AB Aur and HD 163296 are examples of relatively isolated star formation, and are not hindered by confusion with other sources (Henning et al. 1998; Di Francesco et al. 1998).

Both AB Aur and HD 163296 show a rich, variable, emission-line spectrum from the ultraviolet to the optical. With a few exceptions, most noticeably the observed \[O i\] emission, these lines have been successfully modelled as arising in an inhomogeneous stellar wind (Catala et al. 1989; Böhm et al. 1996; Bouret et al. 1997). Infall of material has been detected in both HD 163296 and AB Aur through monitoring of UV and optical absorption and emission lines and can be explained by the presence of infalling evaporating exocomets (Grady et al. 1996, 1999). In the infrared, both stars are among the sources with the strongest 10 $`\mu `$m silicate feature in emission (Cohen 1980; Sitko 1981; Sorrell 1990). Sitko et al. (1999) compared the 10 $`\mu `$m silicate feature in HD 163296 with that of solar-system comets and found a striking resemblance with that of comet Hale-Bopp. Neither star shows the 3.29 $`\mu `$m UIR emission band present toward many Herbig Ae/Be stars (Brooke et al. 1993; Sitko et al. 1999). Basic astrophysical parameters of both stars are listed in Table 1.

There is strong evidence for the presence of a circumstellar disk in both AB Aur and HD 163296. Bjorkman et al. (1995) detected a 90° flip of the polarization angle between the optical and the ultraviolet in HD 163296, which they interpreted as evidence for a flattened, disk-like structure. Mannings & Sargent (1997) resolved the gaseous disks surrounding AB Aur and HD 163296 using CO millimeter wave aperture synthesis imaging. Using continuum measurements at 1.3 mm, the same authors also resolved the circumstellar dust disk of HD 163296. The AB Aur dust disk has also been resolved in the infrared, and shows a surprisingly strong dependence of disk diameter on wavelength, ranging from 0$`\stackrel{}{.}`$0065 (0.94 AU) at 2.2 $`\mu `$m (Millan-Gabet et al. 1999), through 0$`\stackrel{}{.}`$24 (35 AU) at 11.7 $`\mu `$m to 0$`\stackrel{}{.}`$49 (70 AU) at 17.9 $`\mu `$m (Marsh et al. 1995).

In this paper we present new infrared spectra of AB Aur and HD 163296 obtained with the Short- and Long Wavelength Spectrometers on board the Infrared Space Observatory (ISO; Kessler et al. 1996). We will discuss these spectra and their implications for the evolution of dust in Herbig systems. In a subsequent paper (Bouwman et al. 2000), we will describe a model for the circumstellar dust disks of AB Aur and HD 163296 and apply it to these data.

## 2 Observations

ISO Short Wavelength (2.4–45 $`\mu `$m) Spectrometer (SWS; de Graauw et al. 1996) and Long Wavelength (43–197 $`\mu `$m) Spectrometer (LWS; Clegg et al. 1996) full grating scans of AB Aur were obtained in ISO revolutions 680 (at JD 2450717.747) and 835 (JD 2450872.380), respectively. An SWS full grating scan of HD 163296 was made in revolution 329 (JD 2450367.398). Observing times were 3666 seconds for the SWS and 2741 seconds for the LWS observations. Data were reduced in a standard fashion using calibration files corresponding to OLP version 7.0 (SWS) or 6.0 (LWS), after which they were corrected for remaining fringing and glitches. To increase the S/N in the final spectra, statistical outliers were removed and the detectors were aligned, after which the spectra were rebinned to a lower spectral resolution. The resulting spectra are shown in Fig. 1.

$`N`$-band (10.1 $`\mu `$m) images of HD 163296 were obtained on July 25, 1997 (at JD 2450654.608) using TIMMI on the ESO 3.6m telescope at La Silla. Total integration time was 65 minutes. The pixel size was 0$`\stackrel{}{.}`$336, with a total field of view of 21$`\stackrel{}{.}`$$`\times `$ 21$`\stackrel{}{.}`$5. After a standard reduction procedure, the HD 163296 image was indistinguishable from that of the standard star $`\eta `$ Sgr. After deconvolution of the HD 163296 image with that of $`\eta `$ Sgr, the resulting stellar image stretched across 2 pixels. We conclude that the bulk of the 10 $`\mu `$m flux of HD 163296 comes from an area less than 0$`\stackrel{}{.}`$7 (90 AU at 122 pc) in diameter.

## 3 Contents of spectra

The infrared spectra of AB Aur and HD 163296 (Fig. 1) show big differences: whereas AB Aur shows the cool, strong continuum expected for a Herbig star, for HD 163296 the continuum appears to be so weak that the entire SWS spectrum is dominated by solid-state emission features. The IRAS fluxes also plotted in Fig. 1 suggest that an underlying continuum in HD 163296 is present, but peaks longward of 100 $`\mu `$m and therefore is much cooler than the $``$40 K continuum in AB Aur.

Both AB Aur and HD 163296 show a strong 9.7 $`\mu `$m amorphous silicate feature in emission together with a broad emission complex ranging from 14 to 38 $`\mu `$m. The emissivities for various dust components in the spectra are also included in Fig. 1. In the HD 163296 spectrum the broad emission complex ranging from 14 to 38 $`\mu `$m is too broad and too intense to be solely attributed to the 19 $`\mu `$m feature due to the O–Si–O bending mode. We tentatively attribute this feature to a blend of silicates and iron oxide. Lab spectra of FeO also show a strongly rising emissivity in the short-wavelength range of the SWS (Henning et al. 1995). When folded with a $``$ 800 K blackbody, this naturally produces a broad emission feature peaking around 3 microns, which is present in both stars.

The large degree of redundancy in the SWS data makes it possible to assess the reality of weak spectral features which at first glance may appear to be lost in the noise. Each part of the spectrum was scanned twice by twelve detectors, so by checking whether a particular feature is seen in all detectors and in both scan directions, it is possible to disentangle real features from noise. The features identified in this way are listed in Table 2. AB Aur clearly shows the familiar 6.2 and 7.7 and possibly also the 8.6 and 11.2 $`\mu `$m UIR bands usually attributed to emission by polycyclic aromatic hydrocarbons (PAHs), as well as a new UIR band at 15.9 $`\mu `$m. The 3.29 $`\mu `$m UIR band is absent.

The HD 163296 spectrum shows a number of small emission features at wavelengths corresponding to those of crystalline olivines ((MgxFe1-x)2SiO4). Remarkably, the amorphous silicate feature has a higher band-strength relative to the continuum in our ISO data than in the ground-based 8–13 $`\mu `$m spectrum of HD 163296 by Sitko et al. (1999). In contrast, the 11.2 $`\mu `$m shoulder, due to crystalline silicates, appears weaker in our spectrum, suggesting a significant time variability of both components.

In addition to this, HD 163296 shows emission from the 44 $`\mu `$m H2O ice feature. The relative location of the IRAS 60 $`\mu `$m measurement in comparison to the SWS spectrum suggests that the long-wavelength H2O ice feature around 69 $`\mu `$m as well as the broad unidentified feature longward of 100 $`\mu `$m, observed in HD 100546 and HD 142527 (Malfait et al. 1998, 1999), might also be very prominent in HD 163296. PAHs are very weak or absent in HD 163296.

In addition to the solid-state features, both AB Aur and HD 163296 also contain a number of H i recombination lines at shorter wavelengths. All lines from the Bracket and Pfund included in the SWS wavelength range series are present, while the higher H i series are not detected, possibly due to the combined effect of a lower instrumental sensitivity and a higher background in this part of the spectrum. These recombination line data will be discussed in more detail in a forthcoming paper.

The LWS spectrum of AB Aur is relatively smooth and featureless. The only line that is clearly visible in the spectrum is the \[C ii\] line at 157.7 $`\mu `$m. The strength of this line (8.6 $`\times `$ 10-16 W m-2) is compatible with it originating in the background rather than being circumstellar. As can be seen from Fig. 1, there is a $``$25% difference in the flux scales between the AB Aur SWS and LWS spectra in the overlapping region. Although within the formal errors of the absolute flux calibration for SWS and LWS, this discrepancy is larger than that found in other sources. The difference cannot be attributed to the different aperture sizes (33″ $`\times `$ 20″ for SWS versus a circular 80″ FWHM for LWS), and confusion with extended emission, since then the LWS spectrum would have to have a higher flux level than the SWS spectrum. It is interesting to note that in the time interval between the SWS and LWS measurements, AB Aur did show an optical photometric event which could have also affected the infrared brightness (van den Ancker et al. 1999 and references therein), so the difference in flux between SWS and LWS might in fact be due to real variability. In Fig. 1 we also plotted the IRAS fluxes of AB Aur. Although compatible with both spectra, the IRAS 60 $`\mu `$m flux agrees better with the SWS spectrum, so we rather arbitrarily choose to adopt the SWS flux calibration for the region around 45 $`\mu `$m.

## 4 Spectral energy distributions

Spectral Energy Distributions (SEDs) of AB Aur and HD 163296 were constructed from literature data as well as our new ISO spectra and newly obtained VLA photometry for HD 163296 and are shown in Fig. 2. All the submm fluxes in these SEDs refer to single-dish measurements. As can be seen from Fig. 2, the SED can be naturally decomposed in three parts: the optical wavelength range, where the total system flux is dominated by the stellar photosphere, the infrared to submm, where emission originates from the circumstellar dust disk, and the radio, where free-free emission from the stellar wind becomes dominant.

The difference in behaviour of the dust component in HD 163296 and AB Aur is striking: after a nearly flat energy distribution in the infrared, the sub-mm and mm fluxes of AB Aur drop rapidly ($`\lambda F_\lambda \lambda ^{4.3}`$), indicative of the dust becoming optically thin at these wavelengths, whereas toward HD 163296 the slope of the sub-mm to mm fluxes ($`\lambda ^{2.9}`$) is within errors equal to that of the Rayleigh-Jeans tail of a black body ($`\lambda ^3`$). The new radio points at 1.3, 3.6 and 6 cm for HD 163296 do not follow the simple power-law dependence expected if these were solely due to free-free radiation. This demonstrates that even at wavelengths as long as 1.3 cm, a significant fraction of the system flux is due to circumstellar dust. The 3.6 and 6 cm fluxes are probably dominated by free-free emission.

The energy distribution of a circumstellar dust disk is governed by its temperature profile, the density distribution and the dust properties (chemical composition and size distribution). Since the circumstellar disks of AB Aur and HD 163296 stars are expected to be passive (Waters & Waelkens 1998) and the properties of the central stars are nearly identical, the temperature profiles in the disks are expected to be similar as well. One possibility to explain the different sub-mm to mm slope for AB Aur and HD 163296 could be a much flatter density distribution for HD 163296. However, with a standard sub-mm dust emissivity ($`\beta `$ = 2) the inferred dust mass for HD 163296 would become implausibly large. A better explanation may be that the dust properties of AB Aur and HD 163296 are different, a fact already concluded independently from the ISO spectra. To be able to radiate efficiently, dust particles must have a size similar to (or larger than) the wavelength, pointing to the existence of a population of mm- to cm-sized cold dust grains in the circumstellar environment of HD 163296, whereas those in AB Aur must be micron-sized. The ISO spectrum of HD 163296 also contains warm ($``$ 800 K) dust, suggesting a significant lack of emission from dust of intermediate temperatures.

## 5 Discussion and conclusions

We have shown that the main difference between the AB Aur and HD 163296 systems is that HD 163296 contains a population of very large (mm to cm-sized), cold, partially crystalline, dust grains, which is absent in AB Aur. AB Aur contains a population of small dust grains (PAHs), which is absent in HD 163296. In view of the fact that the stellar parameters are nearly identical except for stellar rotation, these differences are remarkable. This must mean that either the evolution of the dust composition in protoplanetary disks happens within the error in the age determination of both systems (2$`{}_{1}{}^{}{}_{}{}^{+2}`$ Myr for AB Aur vs. 4$`{}_{2}{}^{}{}_{}{}^{+4}`$ Myr for HD 163296; van den Ancker et al. 1998), or that the evolution of the dust is dominated by external factors.

We have also shown that in AB Aur and HD 163296 iron oxide is a constituent of the circumstellar dust mixture and can also be responsible for the observed excess emission near 3 $`\mu `$m. Since the near-infrared excess exhibited by our programme stars is by no means unusual for a Herbig Ae/Be star, this means that the same applies for the entire group of Herbig stars. Therefore the results of models attributing this near infrared excess emission to very hot dust from an actively accreting disk (e.g. Hillenbrand et al. 1992) must be regarded with some caution. The case of HD 163296 demonstrates that infrared broad-band photometry can be completely dominated by emission from solid-state features, which must also be taken into account in any future modelling of the energy distributions of Herbig stars.

The detection of PAHs in AB Aur shows that ground-based surveys (e.g. Brooke et al. 1993) have underestimated the fraction of Herbig stars containing PAHs. The case of HD 163296 shows that the fraction of Herbig stars showing PAH emission will not go up to 100%, so models depending on the presence of very small dust grains to explain the observed near-infrared excess in Herbig Ae/Be stars (Natta et al. 1993; Natta & Krügel 1995) will not be successful in all cases. We believe iron oxide to be a more plausible explanation for this near-infrared excess.

It is difficult to attribute the absence of the 3.29 $`\mu `$m feature in AB Aur to an unusual temperature of the PAHs. Since the 6.2 and 7.7 $`\mu `$m C–C stretches are strong, while the bands due to C–H bonds are weak or absent, a more promising possibility seems a very low hydrogen covering factor of the PAHs in AB Aur or the presence of a population of large ($`>`$ 100 C atoms) PAHs (Schutte et al. 1993). If the new 15.9 $`\mu `$m UIR band in AB Aur is also caused by PAHs, this suggests that it is also due to a C–C bond.

To gain more insight in the evolution of dust in protoplanetary disks, it is useful to compare the spectra presented here to those of other Herbig Ae stars. Except for the absence of crystalline material, the AB Aur spectrum and energy distribution are nearly identical to those of its older counterpart HD 100546 (Malfait et al. 1998). In the case of AB Aur (2.5 M), any possible crystallization of circumstellar dust must therefore occur at a stellar age older than 2 $`\times `$ 106 years. The differences between HD 100546 and HD 163296 are larger: the crystalline dust in HD 100546 is much more prominent than that in HD 163296 and the population of large cold dust grains seen in HD 163296 is absent in HD 100546.

As these cases of AB Aur, HD 163296 and HD 100546 demonstrate, the age of the central star and the degree of crystallization do not show a one to one correspondence, but the processes of grain growth and crystallization in protoplanetary disks are also not necessarily coupled. Studying a larger sample of Herbig Ae/Be stars might shed more light on what causes this large observed diversity in dust properties in systems which appear very similar in other aspects.

###### Acknowledgements.

This paper is based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with the participation of ISAS and NASA and on observations collected at the European Southern Observatory, La Silla, Chile. The authors would like to thank the SWS IDT for their help with the SWS observations and Norman Trams and Michelle Creech-Eakman for their help with the LWS observations. Koen Malfait and Xander Tielens are gratefully acknowledged for reading of the manuscript of this paper prior to publication.


text:
# Topological approach to Luttinger’s theorem and the Fermi surface of a Kondo lattice

## Abstract

A non-perturbative proof of Luttinger’s theorem, based on a topological argument, is given for Fermi liquids in arbitrary dimensions. Application to the Kondo lattice shows that even the completely localized spins do contribute to the Fermi sea volume as electrons, whenever the system can be described as a Fermi liquid.

Landau’s Fermi liquid theory is among the most important theories in quantum many-body problem. At zero temperature, a Fermi liquid has a Fermi surface, similarly to the noninteracting fermions. One of the most fundamental results on the Fermi liquid is Luttinger’s theorem, which states that the volume inside the Fermi surface is invariant by the interaction, if the number of particles is held fixed. Luttinger argued, in his 1960 paper , the correction to the volume vanishes order by order in the perturbation expansion.

Recently, there have been renewed interests in the Luttinger’s theorem. Since Luttinger’s original proof was based on the perturbation theory, Luttinger’s theorem could be violated by non-perturbative effects. In fact, several claims of possible breakdown of Luttinger’s theorem have been reported recently . On the other hand, such non-perturbative effects can violate the Fermi liquid theory itself. In fact, the Fermi liquid theory is known to be invalid generally in one dimension, where Tomonaga-Luttinger (TL) liquid is the generic behavior. Although not being a Fermi liquid, a TL liquid in one dimension has a well-defined Fermi surface (actually Fermi points in one dimension). Thus the question of the validity of Luttinger’s “theorem” still exists in this case, where Luttinger’s original proof certainly does not apply. This question was answered recently by a perturbative proof  and a more general non-perturbative proof , which can be applied to one-dimensional TL liquids. However, the question on higher (especially two) dimensions remains unanswered. In fact, it is not clear whether a Fermi liquid which violates Luttinger’s theorem can exist.

Another interesting problem, which is not answered by the Luttinger’s perturbative proof, is the Fermi surface of the Kondo lattice. The Kondo lattice contains a periodic array of localized spins which are coupled to conduction electrons. The Kondo lattice is believed to belong to the Fermi liquid (or TL liquid, in one dimension) in some region of the phase diagram. Even if we assume Luttinger’s theorem to be valid, there is a problem in how to count number of particles. It is rather difficult, by conventional methods, to clarify whether a localized spin should be counted as an electron (“large Fermi surface” picture) or not (“small Fermi surface” picture). In one dimension, the non-perturbative proof of the Luttinger’s theorem was also applied to the Kondo lattice , to show that the localized spins do participate in the Fermi sea. (See also  for numerical evidences.) On the other hand, there has been no definite answer for higher dimensions, although there are several results  supporting the “large Fermi surface” picture.

The argument in Ref.  is a generalization of Lieb-Schultz-Mattis (LSM) theorem , which was given at about the same time as the apparently unrelated Luttinger’s theorem . Since the LSM argument itself cannot be applied to higher dimensions, the discussion in Ref.  was restricted to one dimension. However, very recently the LSM argument was combined with Laughlin’s gauge invariance argument  on the Quantum Hall Effect (QHE) and extended to higher dimensions . Inspired by this observation, we will extend the non-perturbative proof of the Luttinger’s theorem to arbitrary dimensions in the present letter.

We consider an interacting fermion system on a $`D`$-dimensional lattice with periodic boundary conditions. We start from a finite system of size $`L_x\times L_y\times \mathrm{}\times L_D`$, where the length is defined so that the unit cell has the size $`1\times 1\times \mathrm{}\times 1`$. The number of fermions is assumed to be conserved. If the system satisfies a commensurability condition, it can have a finite excitation gap . In this letter, we will rather focus on the gapless case, which is expected for general incommensurate particle density. For simplicity, let us first start with the case of spinless fermion of single species. We introduce a fictitious electric charge $`e`$ for each particle, and a coupling to an externally controlled fictitious electromagnetic field. Because of the periodic boundary condition, the system is topologically equivalent to torus. Following Refs. , we consider an adiabatic increase of a (fictitious) magnetic flux $`\mathrm{\Phi }`$ piercing through the “hole” of the torus so that the uniform electric field is induced, say, in the $`x`$-direction.

While in general the Hamiltonian of the system $`H(\mathrm{\Phi })`$ depends on the flux $`\mathrm{\Phi }`$ reflecting the Aharanov-Bohm (AB) effect, the AB effect is absent when the flux reaches the unit flux quantum $`\mathrm{\Phi }_0=hc/e`$. We consider the adiabatic increase of the flux from $`\mathrm{\Phi }=0`$ to $`\mathrm{\Phi }=\mathrm{\Phi }_0`$. In the following, we will consider how the total momentum of the system is changed during the adiabatic process in two different ways, and compare those results. In the remainder of this letter, we take the units in which $`\mathrm{}=1`$, for simplicity.

First, we analyze the momentum change in a system of interacting fermions for general. We remind the reader that the momentum itself is a gauge dependent quantity in the presence of the gauge field; a meaningful comparison between momenta can only be made under the same gauge choice. In a simplest gauge choice, the AB flux $`\mathrm{\Phi }`$ is represented by the uniform vector potential $`A_x=\mathrm{\Phi }/L_x`$ in the $`x`$-direction. In this gauge, the Hamiltonian always commutes with the translation operator $`T_x`$ to the $`x`$-direction. We further assume that the translation symmetry is not spontaneously broken, as it should not in a Fermi liquid. Thus the ground state is an eigenstate of the total momentum $`P_x`$: $`P_x|\mathrm{\Psi }_0=P_x^0|\mathrm{\Psi }_0`$ with the eigenvalue $`P_x^0`$. The $`x`$-component of the total momentum $`P_x`$ is related to $`T_x`$ as $`T_x=e^{iP_x}`$. After the adiabatic process, the original groundstate $`|\mathrm{\Psi }_0`$ evolves into some state $`|\mathrm{\Psi }_0^{}`$. While the state $`|\mathrm{\Psi }_0^{}`$ could be different from $`|\mathrm{\Psi }_0`$, it belongs to the same eigenvalue $`P_x^0`$ of $`P_x`$, because the Hamiltonian always commutes with $`T_x`$ (and thus $`P_x`$) in the uniform gauge  during the adiabatic process. Although it naively means that the momentum is unchanged after the adiabatic process, it is not true. The Hamiltonian $`H(\mathrm{\Phi }_0)`$ with the unit flux quantum in the uniform gauge is different from the original one $`H(0)`$, although the spectrum should be identical. Namely, they correspond to different choices of the gauge for the same physics. In order to get back to the original gauge, we must perform a large gauge transformation

$$U=\mathrm{exp}[\frac{2\pi i}{L_x}\underset{\stackrel{}{r}}{}xn_\stackrel{}{r}],$$
(1)
where $`n_\stackrel{}{r}`$ is the particle number operator at site $`\stackrel{}{r}`$, and $`x`$ is the $`x`$-coordinate of $`\stackrel{}{r}`$. This transforms the Hamiltonian $`H(\mathrm{\Phi }_0)`$ back to the original one: $`UH(\mathrm{\Phi }_0)U^1=H(0)`$. After this gauge transformation, the adiabatic evolution of the groundstate becomes $`U|\mathrm{\Psi }_0^{}`$.

Now we can examine the total momentum $`P_x`$ of this state, and compare it with the original one $`P_x^0`$. Here we can employ the arguments used in the LSM theorem and its generalizations . By using the identity

$$U^1T_xU=T_x\mathrm{exp}[2\pi i\underset{\stackrel{}{r}}{}\frac{n_\stackrel{}{r}}{L_x}]$$
(2)
we see that $`U|\mathrm{\Psi }_0^{}`$ is an eigenstate of $`P_x`$ with

$$P_x=P_x^0+2\pi \nu L_yL_z\mathrm{}L_D,$$
(3)
where $`\nu `$ is the particle density (number of particles per unit cell). This result is valid regardless of the interaction strength.

Next, we analyze the momentum change assuming that the system is a Fermi liquid. The Fermi liquid is described in terms of quasiparticles, which are almost non-interacting. More precisely, the low-energy effective Hamiltonian of a Fermi liquid is given by

$$\underset{\stackrel{}{k}}{}ϵ(\stackrel{}{k})\stackrel{~}{n}_\stackrel{}{k}+\underset{\stackrel{}{k},\stackrel{}{k}^{}}{}f(\stackrel{}{k},\stackrel{}{k}^{})\stackrel{~}{n}_\stackrel{}{k}\stackrel{~}{n}_\stackrel{}{k}^{},$$
(4)
where $`\stackrel{~}{n}_\stackrel{}{k}`$ is the quasiparticle number operator of momentum $`\stackrel{}{k}`$. Namely, there is an interaction energy due to the second term but no scattering between the quasiparticles. Thus the eigenstates of $`\stackrel{~}{n}_\stackrel{}{k}`$ are also eigenstates of Hamiltonian. In the groundstate, the Fermi sea (region inside the Fermi surface) is completely filled with quasiparticles, while the outside is empty in terms of quasiparticles. Excitations on the groundstate are given by quasiparticles outside the Fermi sea and/or quasiholes inside the Fermi sea. In fact, the quasiparticle (or quasiholes) are free from scattering only in the vicinity of the Fermi surface; the very notion of quasiparticle/hole is useful only in this case. The Fermi liquid theory is valid for the low-energy phenomena, in which the relevant excitations consist only of quasiparticles (quasiholes) near the Fermi surface.

Let us define the Fermi sea volume $`V_F^{(L)}`$ in the finite size system $`L_x\times L_y\times \mathrm{}\times L_D`$. The quasiparticles are scattering free, and their momenta are discretized as in the case of free particles. Thus we can define the Fermi sea volume $`V_F^{(L)}`$ by an integer “occupation number” of the quasiparticles $`N_F^{(L)}`$:

$$V_F^{(L)}=\frac{(2\pi )^DN_F^{(L)}}{L_xL_y\mathrm{}L_D}.$$
(5)
Although the quasiparticles are not free from scattering (and thus are not meaningful) away from the Fermi surface, this expression is still valid because the Fermi sea volume is uniquely determined by its surface. The $`V_F^{(L)}`$ should approach the true volume of the Fermi sea $`V_F`$, in the thermodynamic limit $`L_j\mathrm{}`$.

The adiabatic evolution is determined by the low-energy effective Hamiltonian (4). In the Fermi liquid theory, the charge of the quasiparticle is identical to that of the original particle $`e`$. The coupling of the quasiparticles to the uniform vector potential $`A_x`$ is thus given by the substitution of the momentum $`k_xk_x+eA_x/c`$ in the Hamiltonian. After the adiabatic insertion of the unit flux quantum, and getting back to the original Hamiltonian by the gauge transformation, each quasiparticle gets a momentum shift: $`k_x`$ is increased by $`2\pi /L_x`$. This produces quasiparticles on one side of the Fermi surface, and quasiholes on the opposite side.

Since the result of the adiabatic process is equivalent to the shift of the whole Fermi sea by $`2\pi /L_x`$, the change of the $`x`$-component of total momentum $`P_x`$ of the system during the adiabatic process is given by

$$\mathrm{\Delta }P_x=\frac{2\pi }{L_x}N_F^{(L)}$$
(6)
We note that the only changes after the adiabatic process involve the quasiparticles and quasiholes near the Fermi surface, so that the Fermi liquid theory is still valid. To violate eq. (6), the system must break some of the properties of Fermi liquid used in the present argument. For example, if a quasiparticle had a charge $`e^{}`$ which is different from the charge $`e`$ of the original particle, we would obtain a different result.

Now, comparing the two results eqs. (3) and (6) obtained with different arguments, we obtain $`N_F^{(L)}/L_x\nu L_yL_z\mathrm{}L_D=\text{(integer)}`$, where we have used the fact that each component of momenta is defined modulo $`2\pi `$. Let us choose the system size so that $`L_x,L_y,\mathrm{}`$ and $`L_D`$ are mutually prime with the others. We also assume $`L_x=ql_x`$ where $`l_x`$ is an integer. (It should be recalled that system size should be an integral multiple of $`q`$, to allow the filling factor $`\nu =p/q`$.) Then, from eq. (5) we obtain $`N_F^{(L)}pl_xL_yL_z\mathrm{}L_D=L_x\times \text{(integer)}`$.

Furthermore, we can consider other adiabatic processes, in which the gauge field is induced in one of the other directions $`y,z,\mathrm{}`$, instead of $`x`$. Similar calculations for these cases lead to $`N_F^{(L)}=L_\alpha \times \text{(integer)}`$, where $`\alpha =y,z,\mathrm{},D`$. Because we have chosen the lengths $`L_j`$’s mutually prime, we conclude that $`N_F^{(L)}pl_xL_yL_z\mathrm{}L_D=nL_xL_yL_z\mathrm{}L_D`$, where $`n`$ is an integer. Writing this in terms of Fermi sea volume, we arrive at

$$\frac{V_F}{(2\pi )^D}\nu =n,$$
(7)
where we have replaced the Fermi sea volume $`V_F^{(L)}`$ for the finite size system by its thermodynamic limit $`V_F`$, because this relation is exact already for the finite system. The thermodynamic limit $`V_F`$ should be independent of our special (mutually prime) choice of $`L_j`$’s, if $`V_F`$ is well-defined.

The relation (7) is nothing but the statement of Luttinger’s theorem. The integer $`n`$ corresponds to the number of completely filled bands. It is valid also when the Fermi sea consists of several disjoint regions, if $`V_F`$ is understood as the sum of volumes of all regions. Our proof is much simpler than the original one . Moreover, in contrast to Ref. , our argument is non-perturbative and relies only on some of the basic properties of Fermi liquid.

It is straightforward to extend our argument to spinful electrons. When the numbers of up-spin electrons and down-spin electrons are conserved separately, we consider the fictitious electromagnetic field coupled to only up-spin (or down-spin) electrons. Assuming the spinful Fermi liquids, the volume of the Fermi sea $`V_F^\sigma `$ for spin $`\sigma `$ is given by $`V_F^\sigma =(2\pi )^D\nu _\sigma `$, where $`\nu _\sigma `$ is the number of particles with spin $`\sigma `$ per unit cell. For the spin-symmetric case $`\nu _{}=\nu _{}`$, it reads $`V_F=V_F^{}=V_F^{}=(2\pi )^D\nu /2`$ where $`\nu `$ is the total particle density $`\nu _{}+\nu _{}`$.

As a nontrivial application, let us consider the Kondo lattice. Luttinger’s original perturbative proof does not apply to this case, and the question on the volume of the Fermi sea has remained. For the sake of clarity, we consider the Kondo lattice model given by the Hamiltonian

$$H=\underset{j,k}{}t_{jk}c_{j\sigma }^{}c_{k\sigma }+\text{h.c.}+\underset{j}{}U_jc_j^{}c_jc_j^{}c_j+\underset{l}{}J_l\stackrel{}{s}_l\stackrel{}{S}_l,$$
(8)
where $`c_{j\sigma }^{}`$ and $`c_{j\sigma }`$ are standard Fermion creation/annihilation operators at site $`j`$ with spin $`\sigma `$, $`\stackrel{}{s}_l=c_{l\alpha }^{}\stackrel{}{\sigma }^{\alpha \beta }c_{l\beta }/2`$ is the spin operator of the conduction electron, and $`\stackrel{}{S}_l`$ is the localized spin at site $`l`$. As in the previous case, we couple the fictitious electromagnetic field only to the up-spin electrons. After the adiabatic insertion of the AB flux of unit flux quantum, we make the gauge transformation as in the previous cases. However, the naive one $`U_{}^e=\mathrm{exp}[\frac{2\pi i}{L_x}_\stackrel{}{r}xn_{\stackrel{}{r}\sigma }]`$, does not bring the Hamiltonian back to the original one, because it changes the Kondo coupling. In order to recover the original Hamiltonian, we must also twist the localized spins. The transformation

$$U_{}=\mathrm{exp}[\frac{2\pi i}{L_x}\underset{\stackrel{}{r}}{}x(n_{\stackrel{}{r}\sigma }+S_\stackrel{}{r}^z)]$$
(9)
does the required job. We obtain the total momentum after the adiabatic process as

$$P_x=P_x^0+2\pi [\nu _{}+N_s(S+m)]L_yL_z\mathrm{}L_D,$$
(10)
where $`N_s`$ is the number of localized spins per unit cell and $`m`$ is the magnetization per single localized spin. The special contribution proportional to $`S`$ comes from the boundary term $`\mathrm{exp}(2\pi iN_sS_1^z)`$ appearing in $`U_{}^1T_xU_{}`$, similarly to the one dimensional case .

Thus, provided that the system belongs to a spinful Fermi liquid, the volume of the Fermi sea is given by $`V_F^\sigma =(2\pi )^D[\nu _\sigma +N_s(S\pm m)]`$, where $`\pm `$ takes $`+`$ for $`\sigma =`$ and $``$ for $`\sigma =`$. For the spin-symmetric case $`\nu _{}=\nu _{}`$ and $`m=0`$, we obtain

$$V_F=V_F^{}=V_F^{}=\frac{(2\pi )^D}{2}[\nu +2N_sS],$$
(11)
for the total particle density $`\nu =\nu _{}+\nu _{}=2\nu _{}`$. This is exactly what we obtain if we apply the Luttinger’s theorem to the Anderson-type model in which the localized spins are represented by electrons. It means that the localized spin $`S`$ does contribute to the Fermi sea volume as $`2S`$ electrons, even though it is completely immobile. This is the picture conventionally called as the “large Fermi surface”.

It should be noted that we did not answer the non-trivial question whether (or when) the Kondo lattice belongs to the Fermi liquid. We have only proved that, if the Kondo lattice is a Fermi liquid (as it is believed to be true in some region of the phase diagram), the localized spins participate in the Fermi sea.

Finally, let us comment on claims  of the violation of the Luttinger’s theorem. There are several possibilities regarding the apparent contradiction to our non-perturbative proof. Of course, it should be checked whether our argument applies to the model under consideration. However, our argument does apply to a very wide range of lattice models, including the Hubbard and $`t`$-$`J`$ models, for which the violation of the Luttinger’s theorem has been also proposed. A possibility is that the system is not a Fermi liquid in these cases. In other words, a violation of the Luttinger’s theorem requires the system to be a non-Fermi liquid. We note that, however, merely not being a Fermi liquid is insufficient, as the TL liquid in one dimension does satisfy the Luttinger’s theorem . Our approach could be extended to a non-Fermi liquid which has an appropriately defined Fermi surface, if such a liquid does exist. Our argument reveals a rigid relationship between the structure of low-energy excitations and the Fermi sea volume.

Another possibility is that the claimed violation of Luttinger’s theorem is actually incorrect. In particular, numerical results are only available for restricted system size and/or temperature, and can miss the possibly small singularity at the true Fermi surface. On the other hand, even if they are incorrect in identifying the true Fermi surface, they might still be of physical relevance because the actual experiments are done also at finite energy scale; the experimentally measured “Fermi surface” could be different from the true Fermi surface defined in the low-energy limit, to which our argument applies. In any case, our definite result on the Fermi surface of the Fermi liquid in the low-energy limit would be useful as a guideline. Claims of the violation of the Luttinger’s theorem should be examined in the light of the present result.

During the forty years after the Luttinger’s paper , several examples of “quantization” of a physical quantity have been found in many-body physics. Namely, despite the complexity of the interacting many-body states, some physical quantity takes a special value which is stable against various perturbations such as interaction strength. Presumably the most natural understanding of such a quantization is given by a topological argument. Indeed, typical examples of the quantization – QHE and the quantized magnetization plateaus have been related to topological mechanisms .

Luttinger’s theorem perhaps does not look like a quantization, because the volume of the Fermi sea takes continuous values depending on the particle density. However, the insensitivity to the interaction resembles to other quantization phenomena, and may well be regarded as a certain kind of quantization, especially when written as in eq. (7). In fact, we have revealed a close theoretical relationship among Luttinger’s theorem, QHE and magnetization plateaus. In addition, our argument can be related also to the chiral anomaly in quantum field theory . Luttinger’s theorem might be actually the first example of the topological quantization discovered in quantum many-body problem, although the topological understanding has been missing for a long time.

I would like to thank Ian Affleck, Hal Tasaki and Masanori Yamanaka for stimulating discussions which were essential for the present work. I am also grateful to an anonymous referee for pointing out a logical gap in the original manuscript. This work is supported by Grant-in-Aid from Ministry of Education, Science, Sports and Culture of Japan.


text:
# Phase transition in the ground state of a particle in a double-well potential and a magnetic field

## I Introduction

It is a well-known fact  that in the absence of a magnetic field the ground state of bosons is non degenerate, and therefore has the symmetry of the hamiltonian. Mathematically this result from the fact that the kernel of the operator $`e^{tH}`$ is positive. This last property no more holds in the presence of a magnetic field so that degeneracy of the ground state may be expected, as well as symmetry breaking in it. One-body systems already may show this phenomenon. Indeed Lavine and O’Carrol  proved the existence of spherically symmetric potentials for which, in the presence of a magnetic field, the ground state has a non-vanishing value for the $`z`$ component of angular momentum, so that the rotational symmetry is broken.

Further examples were provided by Avron, Herbst and Simon <sup>,</sup> . On the opposite side, these last authors were able to prove that for the hydrogen atom the symmetry is not broken, as well as in the case where the potential is monotonically increasing with the distance. These authors, however, mainly concerned with problems of atomic physics, did not discuss the degeneracy and the physical significance of it.

On the other hand, two of us analysing the problem of a particle confined to a disc or an annulus in the presence of a magnetic field found that the ground state was degenerate in the case of an annulus and for a disc with Neumann boundary conditions (with Dirichtlet boundary conditions in the disc case the degeneracy disappears.) The degeneracy appears each time the magnetic field reaches a critical value and the magnetisation jumps at these critical values, which form a discrete set .

Motivated by these results we consider in this article a class of systems for which similar phenomena occur. Namely we analyse the ground state of a particle in three dimensions moving in a double-well type potential, cylindrically symmetric, and submitted to a constant magnetic field in the $`z`$ direction.

We find that the ground state has an azimuthal momentum $`\mathrm{}m`$ taking increasing values $`m=0,1,2,\mathrm{}`$ when we increase the magnetic field $`B`$. At critical values of $`B`$ $`(B_m)`$ the ground state is twice degenerate between the $`m`$ and the $`m+1`$ state. Moreover the magnetisation jumps at these critical values and shows in general an oscillatory behaviour reminiscent of the well known de Haas von Halphen oscillations in solid state physics.

We show that this phenomenon can be understood by an analysis of the minima of the potential energy, fixing however the angular momentum to its quantised value $`\mathrm{}m`$. In the two-dimensional case we can use the WKB method and obtain bounds on the energy in order to estimate the critical fields. But in general, we had to compute numerically the energies and compare them to estimates based on trial wave functions. The agreement is quite good in general.

Concerning possible experimental verifications of these effects, which require basically to have a potential which has a minimum sufficiently far from the origin, we could think of two cases. The first one would be in some molecules where proton dynamics could be described by such an effective potential. The second one, more thrilling, would be the case of charged bosons undertaking a Bose-Einstein condensation. Our results suggest that in this case, the bosons would undertake a phase transition in their condensate, when we apply an increasing magnetic field. This phase transition would manifest itself by appearance of oscillations in the magnetisation, which would jump at certain critical values of the magnetic field.

## II The Model

We will consider the case of a particle of mass $`\mu `$, charge $`q`$, in a potential $`V`$ with a cylindrical symmetry, submitted to a magnetic field $`\stackrel{~}{B}`$ in the $`z`$ direction. We do not consider the effect of the spin of the particle. We choose for a unit of energy $`V_0`$, and length $`r_0`$, both being characteristic of the potential. The dimensionless hamiltonian reads if $`r=\sqrt{x^2+y^2}`$

$$(iϵ\stackrel{}{}\stackrel{}{A})^2+V(r,z)$$
(1)
where

$$ϵ=\frac{\mathrm{}}{r_0\sqrt{2\mu V_0}}$$
(2)
measures the importance of the quantum effects and the vector potential in the symmetric gauge is given by

$$\stackrel{}{A}=(\frac{By}{2},\frac{Bx}{2},0)$$
(3)
$`B={\displaystyle \frac{q}{c}}{\displaystyle \frac{r_0}{\sqrt{2\mu V_0}}}\stackrel{~}{B}`$ being the dimensionless magnetic field.

Thanks to the cylindrical symmetry, we can replace the $`z`$ component of the angular momentum $`L_z`$ by its eigenvalue $`ϵm`$ so that the reduced hamiltonian reads

$$H_m=ϵ^2\left[\frac{1}{r}\frac{}{r}r\frac{}{r}+\frac{^2}{z^2}\right]+\left(\frac{ϵm}{r}\frac{rB}{2}\right)^2+V(r,z)$$
(4)
The ground state energy of this hamiltonian and the corresponding eigenfunction will be denoted $`E_m`$ and $`\psi _m`$.

It remains to specify $`V`$. We will basically consider a double-well potential of the form:

$$V(r,z)=r^4+z^42(r^2+z^2)+vr^2z^2$$
(5)
with $`v`$ satisfying $`v2`$, so that $`V`$ is bounded from below. If $`v`$ is equal to $`0`$ we can decouple the motion in the $`z`$ direction form the one in the plane perpendicular to the magnetic field. This is what we will call the two-dimensional case. If $`v=2`$, we have in three dimensions a potential with spherical symmetry.

We have chosen this double-well form because if we had taken the simple well $`V=r^4+z^4+2(r^2+z^2)+vr^2z^2`$ with $`v0`$ it follows from  that the ground state is not degenerate and corresponds to $`m=0`$.

A physical quantity of interest is the magnetisation in the ground state

$$M=\frac{E}{B}$$
(6)
in units $`\frac{q}{c}r_0\sqrt{\frac{V_0}{2\mu }}`$

We will denote by $`e_m`$ the ground state energy of the hamiltonian

$$h_m=ϵ^2\left[\frac{1}{r}\frac{}{r}r\frac{}{r}+\frac{^2}{z^2}\right]+V_m(r,z)$$
(7)
with

$$V_m=\frac{(ϵm)^2}{r^2}+\frac{B^2}{4}r^2+V$$
(8)
and by

$$E_m=e_mϵmB$$
(9)
the ground state energy of $`H_m`$ given in (4), so that the real ground state energy is given by

$$E=\underset{m0}{inf}E_m$$
(10)
since obviously negative $`m`$ give a larger energy.

Finally we will use the following useful scaling property of the energy $`e_m`$

$$e_m(ϵ,\lambda ,v)=s^2e_m(\frac{ϵ}{s^{3/2}},\frac{\lambda }{s},v)s0$$
(11)
where

$$\lambda =\frac{B^2}{4}2$$
(12)
is the parameter multiplying $`r^2`$ in the potential. Equation (11) follows simply from the scaling transformation : $`r^2sr^2`$ and $`z^2sz^2`$. This relation shows that we have effectively a two parameter dependence of the energy $`e_m`$ in general and a one parameter dependence in the two dimensional case.

The choice $`s=|\lambda |`$ or $`s=m^{\frac{1}{3}}(m1)`$ shows that large magnetic field or large angular momenta correspond to the semi-classical limit. In fact we shall see that in the classical limit $`ϵ0`$ ground state with $`m0`$ are favoured inducing ground state degeneracies at some values of the magnetic field. It thus appears that the tendency to have a ground state with the same symmetry as the hamiltonian and therefore non degenerate is an effect due to quantum mechanics.

## III The classical limit

One can gain some qualitative understanding of the problem by looking at the classical limit of it. This means that we neglect the quantum kinetic energy and define the ground state energy as

$$E=\underset{m0}{inf}\underset{(r,z)}{inf}[V_m(r,z)ϵmB]$$
(13)
where

$$V_m=\frac{(ϵm)^2}{r^2}+r^4+z^4+\frac{B^2}{4}r^22r^2+vr^2z^2$$
(14)
and consider that $`m`$ is an integer.

Two cases need to be considered separately: $`|v|<2`$ and $`v2`$. If $`|v|<2`$ we denote by $`x_m`$ and $`t_m`$ respectively the values of $`r^2`$ and $`z^2`$ which minimise the potential $`V_m`$, and we find

$$t_m=\mathrm{\hspace{0.17em}1}\frac{vx_m}{2}$$
(15)
$$\left(2\frac{v^2}{2}\right)x_m+\left(v2+\frac{B^2}{4}\right)=\frac{(ϵm)^2}{x_m^2}$$
(16)
On the other hand, considering for a while $`m`$ as a continuous variable, the absolute minimum of $`V_mϵmB`$ is given by

$$ϵ\widehat{m}=\frac{B}{2}x_{\widehat{m}}$$
(17)
From (15) this gives an absolute minimum of $`V_mϵmB`$ given by

$$x_{\widehat{m}}=t_{\widehat{m}}=\frac{1}{1+\frac{v}{2}}$$
(18)
and therefore

$$ϵ\widehat{m}=\frac{B}{2}\frac{1}{1+\frac{v}{2}}$$
(19)
In considering the variable $`m`$ as a continuous one we have treated the problem purely classically and the corresponding ”ground state” energy is

$$E^{cl}=\frac{2}{1+\frac{v}{2}}$$
(20)
We know that $`m`$ is a discrete variable but for consistency we must consider $`ϵ`$ as a small number. Then if $`m`$ designates the integer part of $`\widehat{m}`$, we have $`\widehat{m}=m+\theta `$ and if $`0\theta <\frac{1}{2}`$, the ground state has the quantum number $`m`$, whereas if $`\frac{1}{2}<\theta 1`$ it has $`m+1`$.

From this analysis we conclude that if $`B_{m1}<B<B_m`$ where

$$B_m=ϵ(1+\frac{v}{2})(2m+1)$$
(21)
the ground state has the quantum number $`m`$. Hence we see that by increasing the magnetic field, we find in increasing order the values of $`m=0,1,2,\mathrm{}`$ and an infinite set of critical values of the magnetic field exist, $`B_m`$ for which the ground state is twice degenerate, being both $`m`$ and $`m+1`$.
This picture is entirely confirmed by the numerical results in the quantum case. It is also quite interesting to look at the magnetisation. In the state whose quantum number is $`m`$, we have

$$M_m=ϵm\frac{B}{2}x_m$$
(22)
so that using (15)

$$M_m=[ϵm\frac{B}{2}\frac{1}{1+\frac{v}{2}}][\frac{1\frac{v}{2}}{1\frac{v}{2}+\frac{B^2}{4}}]$$
(23)
when $`B_{m1}<B<B_m`$.

This shows that the magnetisation has an ”oscillatory” type of behaviour reminiscent of the familiar de Haas von Halphen one in solid state physics and that the magnetisation jumps at the critical values of the magnetic field, the jump being given by

$$\mathrm{\Delta }M_m=ϵ\frac{1\frac{v}{2}}{1\frac{v}{2}+\frac{B^2}{4}}$$
(24)
Once again this general behaviour is reproduced by the numerical results in the quantum case and the spacing between the values of the critical field is rather well represented by formula (21) when $`m1`$. In the two-dimensional case, i.e. $`v=0`$ and neglecting the trivial $`z`$ dependence, we can proceed further and look at a really semi-classical approximation namely WKB, for the ground state energy

$$_r_{}^{r_+}dr\sqrt{e_mV_m(r)}=\frac{ϵ\pi }{2}$$
(25)
where

$$V_m(r)=\frac{(ϵm)^2}{r^2}+r^4+\left(2+\frac{B^2}{4}\right)r^2$$
(26)
and the ground state energy is $`E_m=e_mϵmB`$.

In fact this WKB approximation will give the best analytical results, apart from the variational estimates for the energy, which give unfortunately only exact upper bounds on the energy.

When the potential has spherical symmetry $`v=2`$, quantum effects are much more important and the classical analysis gives only that the ground state has $`m=0`$ if $`B<2ϵ`$, is degenerate between $`m=0`$ and $`m=1`$ when $`2ϵB<4ϵ`$, has possibly $`m=0,1,2`$ for $`4ϵB<6ϵ`$ and so on. This only suggests that we have again the increasing sequence of $`m`$, when we increase the magnetic field and that critical values appear near $`2ϵm`$.

When $`v>2`$, we find that $`m=0`$ is the ground state except when $`B=2ϵm`$, where it is degenerate between $`m`$ and $`0`$. We may note however that the classical ground state correspond to points $`(r=0,z=\pm 1)`$ in configuration space for $`m=0`$, whereas it corresponds to two circles $`(r=\frac{ϵ}{2B},z=\pm \sqrt{1\frac{ϵ}{2B}})`$ for $`m=1`$ and $`2ϵ<B<4ϵ`$, so that the wave function can be more spread in the $`m=1`$ state than the in the $`m=0`$ state, and that the kinetic energy of the $`m=1`$ state is lower, favouring the $`m=1`$ state. Hence we should expect, at least when $`ϵ`$ is small, a ground state with $`m=0`$ for small fields and a ground state with $`m=1`$, when $`2ϵ<B<4ϵ`$. A similar argument can be given for the higher values of $`m`$.

Finally, it is worth noticing that if we had taken a simple well type potential

$$V(r,z)=r^4+z^4+2(r^2+z^2)+vr^2z^2$$
(27)
the classical analysis gives a ground state with $`m=0`$, at least when $`v1`$. This is a correct result when $`v0`$ at the quantum level.

## IV Numerical results and variational bounds

It is quite useful to undertake a numerical analysis of this problem. We have used a finite element method, choosing for the basis a product of two triangles functions. We discuss separately the two-dimensional problem and the three dimensional ones.

### A Two dimensions

We first give pictures of the ground state energy for two typical values of $`ϵ`$, a small ($`ϵ=0.03`$) and a large one ($`ϵ=0.5`$) as a function of the magnetic field $`B`$. (figure 1). The cusps at the critical values of $`B`$ indicate a jump of the corresponding magnetisation.

This last quantity shows first a diamagnetic behaviour at small field, but then a paramagnetic - diamagnetic oscillation at least when $`ϵ0.3`$. Beyond this value the magnetisation is entirely negative (figure 1 bottom right). We can also note that when $`B`$ becomes large the magnetisation tends to $`ϵ`$, its value in the Landau regime.

The results clearly indicate that we go progressively from the states with $`m=0,1,2\mathrm{}`$ by increasing the magnetic field and that the magnetisation jumps at the critical values. The effect is more pronounced in the classical regime. All these results are in qualitative argument with the classical picture presented before and the agreement is even quantitative when $`ϵ=0.03`$ for example.

The jumps of the magnetisation given by formula (23) are reproduced (figure 2) with a precision of less than $`1`$ percent when $`ϵ=0.03`$, and the spacing between the critical values of the magnetic field

$$\frac{B_{m+1}B_m}{ϵ}=\mathrm{\hspace{0.17em}2}+\mathrm{\Delta }_m$$
(28)
is given by $`\mathrm{\Delta }_m0.04`$ if $`m1`$ and $`ϵ=0.1`$. $`\mathrm{\Delta }_m`$ decreases when $`m`$ increases in agreement with the scaling relation $`B_m=(2m+1)ϵ`$, so that the simple classical formula reproduces rather well the results. By contrast, the jump between the $`m=0`$ and the $`m=1`$ state is largely of quantum mechanical origin, as well as the precise values of the critical fields.

Figure 3 describes the various regions in the $`ϵB`$ plan. We can note that even when $`ϵ>0.25`$ a linear relation exists between $`B_m`$ and $`ϵ`$, as in the classical regime, which is a bit surprising.

It is also interesting to look at the eigenfunctions when the magnetic field reaches its critical value. In figure 4 we give pictures of them at the critical value between the state $`m=0`$ and $`m=1`$ when $`ϵ=0.2`$. We see that their maxima are located very near the minimum of the potential.

Finally we compare the results with two theoretical estimates: first of all the WKB one, and a variational one. This last estimate is based on the following two parameters trial wave function

$$\psi _m=r^me^{\alpha r^2\beta (r1)^2}$$
(29)
The variational upper bound on the energy can be expressed in terms of Weber cylindrical functions, but we directly computed the corresponding integrals.

Tables I,II,III and IV give a comparison of the results for two values of the parameter $`ϵ`$, and for the critical fields. Excellent agreement is found for the variational method (maximal error of the order of 2 % when $`ϵ=0.5`$). WKB works quite well when $`ϵ`$ is small ($`ϵ=0.03`$) as expected, but even better on the energies when $`ϵ=0.5`$ and the error does not exceed 1%.

### B Three dimensions

For the spherically symmetric potential $`(v=2)`$, figure 5 gives the ground energies a well as the corresponding magnetisation for two different values of $`ϵ`$ : 0.03, 0.5.

Once again we see that the values of $`m`$ in the ground state increases with $`B`$, and that the magnetisation jumps at critical values $`B_m`$ of the magnetic field, where the ground state is doubly degenerate. These results are in qualitative agreement with the classical analysis. Figure 6 summaries the results in the $`ϵ`$ \- $`B`$ plane. Notice that in this

case, when $`ϵ0.1`$ already the relation between $`B_m`$ and $`ϵ`$ is no more linear. On the other hand the spacing between the critical values of $`B`$ predicted by the crude classical estimate:

$$\mathrm{\Delta }B_m=B_{m+1}B_m\mathrm{\hspace{0.17em}2}ϵ$$
(30)
is satisfied with a precision of $`25\%`$ at $`m=1`$ and becomes more accurate when $`m`$ increases, at least in the range $`ϵ0.1`$.

Our best variational estimate for the energy was made with a three parameter trial wave function

$$\psi _{\alpha ,\beta ,\zeta }=r^me^{\alpha r^2\beta (\sqrt{r^2+z^2}\zeta )^2}$$
(31)
Table V gives the values of the critical field $`B_m`$ and Table VI the corresponding ground state energies, when $`ϵ=0.05`$ estimated by the variational method and computed with the simulation.

Obviously there is a very good agreement, since the largest error for $`B_m`$ is less than $`2\%`$ and for $`E_m`$ less than $`0.7\%`$. Table VI gives the same but for $`ϵ=0.5`$. Again we see a good agreement (error less than 5%). When $`ϵ`$ increases we found that $`\alpha `$ increases and $`\beta `$ decreases as well as $`\zeta `$ and our trial wave function becomes less accurate, because the double-well nature of the potential is less important compared to the kinetic energy.

Figure 7 describes the situation in the $`v`$ \- $`B`$ plane for $`m=0,1,\mathrm{},10`$ and different values of $`ϵ`$. We notice that when $`v`$ is less than $`2`$ and $`ϵ`$ is not too large ($`ϵ0.2`$), the situation is similar to the one already discussed, but that there is an abrupt change at $`v=2`$ when $`ϵ`$ is small in agreement with the classical analysis. However when $`ϵ>0.2`$ the ground state $`m=0`$ is definitely favoured as $`v`$ increases.

Figure 8 shows the energies for the first five $`m`$ values computed with three different $`v`$: two-dimensional ($`v=0`$), spherical potential ($`v=2`$), and $`v=3`$. We can see a new crossing between the $`m=0`$ and the other $`m`$ levels when $`v`$ becomes larger than $`2`$, although this does not concern the ground state.

## V Bounds on the critical field in the two dimensional case

One might desire to get rigorous upper and lower bounds on the critical fields. One possible approach would consist in getting upper and lower bounds on the ground state energies $`E_m`$. Whereas we have seen that one can obtain very good variational upper bounds, it is rather difficult to get good lower ones. In order to test these results, we analysed only the two-dimensional case.

First we want to obtain conditions under which $`m=0`$ is the ground state. Using the inequality

$$\frac{l^2}{x}+x^2x^2\frac{l^2}{a^2}x+\frac{2l^2}{a}$$
(32)
valid for any $`x`$ and $`a`$ positive, we deduce that

$$e_0[\lambda ]\frac{2(ϵm)^2}{a}+e_0\left[\lambda \left(\frac{ϵm}{a}\right)^2\right]$$
(33)
On the other hand

$$e_0[\lambda ]e_0\left[\lambda \left(\frac{ϵm}{a}\right)^2\right]=_{\lambda \left(\frac{ϵm}{a}\right)^2}^\lambda d\lambda ^{}r^2_0(\lambda ^{})$$
(34)
$$\left(\frac{ϵm}{a}\right)^2r^2_0\left[\lambda \left(\frac{ϵm}{a}\right)^2\right]$$
(35)
since $`r^2_0(\lambda )`$ is decreasing in $`\lambda `$.

But

$$\left|r^2_0[\lambda ]\frac{\lambda }{2}\right|\left[e_0[\lambda ]+\left(\frac{\lambda }{2}\right)^2\right]^{\frac{1}{2}}$$
(36)
The scaling relation and the fact that $`e_0`$ is increasing in $`ϵ`$ imply that when $`\frac{\lambda }{2}1`$

$$e_0[\lambda ]+\left(\frac{\lambda }{2}\right)^2\left(\frac{\lambda }{2}\right)^2(e_0[2]+1)$$
(37)
Taking now $`a`$ such that $`\frac{ϵm}{a}\frac{B}{2}`$ $`(m1)`$ we get combining these inequalities that

$$E_0E_mm1$$
(38)
if we can find $`t>\frac{B}{2}`$ such that

$$t^2\left\{1+\frac{1}{2}\left(t^2\frac{B^2}{4}\right)\right\}\delta 2ϵ\left(t\frac{B}{2}\right)<0$$
(39)
where $`\delta =1+\sqrt{e_o[2]+1}`$

In the estimate for $`\delta `$ we can use our best variational upper bound. Inequality (39) will be satisfied if $`B`$ is less than some value $`B_0`$, so that in this range $`m=0`$ is the ground state. In order to see when $`m0`$ is a ground state, we use the following trial wave function $`\psi (r)`$ for a state with angular momentum $`m^{}`$.

$$\psi (r)=r^{m^{}m}\psi _m(r)m^{}m$$
(40)
where $`\psi _m(r)`$ is the exact ground state wave function for the state with angular momentum $`m`$. An integration by parts shows that

$$_0^{\mathrm{}}drr\left[\psi _{m}^{}{}_{}{}^{2}r^{2(m^{}m)}+2(m^{}m)r^{2(m^{}m)1}\psi _m^{}\psi _m\right]$$
(41)
$$=_0^{\mathrm{}}drr^{2(m^{}m)}\psi _m(r\psi _m)^{}$$
(42)
Therefore if we use the fact that

$$\frac{ϵ^2}{r}(r\psi _m^{})^{}=[V_m(r)e_m]\psi _m$$
(43)
We see that

$`{\displaystyle _0^{\mathrm{}}}`$ $`\mathrm{d}rr\left[ϵ^2\psi _{}^{}{}_{}{}^{2}+V_m^{}(r)\psi ^2\right]=e_m{\displaystyle _0^{\mathrm{}}}drr\psi ^2`$ (44)
$`+ϵ^22m^{}(m^{}m){\displaystyle _0^{\mathrm{}}}drr^{2(m^{}m)1}\psi _m^2`$ (45)
and we conclude that

$$e_m^{}e_m+ϵ^22m^{}(m^{}m)\frac{r^{2(m^{}m1)}_m}{r^{2(m^{}m)}_m}$$
(46)
In particular

$$e_1e_0+2ϵ^2\frac{1}{r^2_0}$$
(47)
If we have a lower bound $`c`$ on $`r^2_0`$ then we see that

$$E_1<E_0$$
(48)
if

$$B>\frac{2ϵ}{c}$$
(49)
We can use for the lower bound $`c`$ the one given in equation (36)

$$c=\frac{\lambda }{2}\sqrt{e_0[\lambda ]+(\frac{\lambda }{2})^2}$$
(50)
which is satisfactory when $`B`$ is not too large, but which becomes negative for large $`B`$. We can repair this by using the fact  that if $`f`$ is an increasing function of $`r`$, its expectation value in the ground state is lowered by adding to the potential a new increasing potential. We can find a useful comparison potential

$$W=a_1r^2+a_2r^4+a_3r^6V$$
(51)
which has a ground state wave function of the form

$$\psi =e^{b_1r^2b_2r^4}b_2>0$$
(52)
so that $`r^2_W`$ can be computed explicitly for this potential and we can take $`c=r^2_W`$ in equation (49), which gives a more satisfactory result for large $`B`$.

In any case we see that the state $`m=1`$ if favoured over the state $`m=0`$ if $`B`$ is larger than some value, and by continuity there must exist a field for which both states have equal energy. But in order to prove that the ground state is $`m=1`$ when $`B`$ is in some range requires to show that $`E_m>E_1`$ $`m2`$. For this purpose let us consider $`m`$ as a continuous parameter. Then

$$\frac{E_m}{m}=\mathrm{\hspace{0.17em}2}ϵ^2m\frac{1}{r^2}_mϵB$$
(53)
If we can show that $`\frac{E_m}{m}0`$ for all $`m1`$, then we will have shown that $`E_m>E_1`$. When $`m1`$ we have

$$\frac{1}{r^2}_m\frac{1}{r^2_m}$$
(54)
and

$$\frac{(ϵm)^2}{r^2_m}+r^2_m^2+\lambda r^2_me_m$$
(55)
In order to get a variational bound on $`e_m`$ we can use the trial wave function $`\psi =r^me^{ar^2}`$, which gives

$$e_m\frac{m+1}{m+2}V_{m+2}(x_{m+2})$$
(56)
where $`x_m`$ is the value of $`x`$ which minimises

$$V_m(x)=\frac{(ϵm)^2}{x}+x^2+\lambda x$$
(57)
Noting that equation (55) implies that

$$r^2_mx_m+\sqrt{e_m+V_m(x_m)}$$
(58)
one can see by combining equations (53), (54), (56) and (58) that $`E_mE_1`$ for all $`m2`$ if

$$\frac{B^2}{8}<\frac{1}{1+c^2}$$
(59)
with

$$c^2=\frac{1}{ϵ^2}|\lambda |\left[x_1+\sqrt{V_{1+2m}(x_{1+2m})V_1(x_1)}\right]^2$$
(60)
which implies that $`B`$ should be less than some value.

We give in the table VII some numerical values for the bounds obtained by these methods.

They show that whereas the range of values of $`B`$ for which $`E_0<E_1`$ and $`E_1<E_0`$ is reasonably estimated for $`ϵ0.1`$, there is no range of values of $`B`$ for which our bounds show that $`m=1`$ is the ground state except when $`ϵ`$ is very small (0.01) But in this range WKB works perfectly well. Obviously we have too poorly estimated the effect of the kinetic energy and that of the centrifugal barrier. Numerical computations for example show that the replacement of $`\frac{1}{r^2}_1`$ by $`\frac{1}{r^2_1}`$ is not appropriate when $`ϵ`$ or $`B`$ are too large.

In conclusion, even in two dimensions improved rigorous bounds on the critical values of the magnetic field are needed, and the WKB method for which we have no estimate of the error gives the best analytic results.

## VI conclusion

It could be of course quite interesting to see an experimental verification of these surprising effects of the magnetic field. Even though we have found them in the case of a double-well, we think that the details of the potential do not matter too much. What is needed is a potential whose minimum is taken sufficiently far from the origin.

We have thought of two possible fields where one could observe such effects. The first one is molecular physics where often the dynamics of electrons or protons is modelled by the motion of a quantum particle in a double-well (although admittedly often a one-dimensional one.) If we consider the case of the electron in the rotationally symmetric double-well, the smallest value of the critical field where the $`m=1`$ and $`m=0`$ states are degenerate is about $`15`$ Tesla if we take for the depth of the potential $`1`$ eV and for the distance to the origin of the minimum 2 $`\AA `$. For protons the situation is more favourable since a field of $`5`$ Tesla can create a degeneracy when the depth is kept to $`1`$ eV and the minimum is at a distance of $`1.5`$ $`\AA `$. Obviously a more detailed investigation is needed if one wants to see these unusual effects (like a change from diamagnietism to paramagnetism) in molecules.

The other field is that of Bose-Einstein condensates of very cold atoms, which recently has made spectacular progress. If we consider free charged bosons in a magnetic field and in a potential $`V(\frac{r}{r_0})`$ one can show that there is a Bose-Einstein condensation in the ground state in three dimensions, in the limit $`r_0`$ going to infinity, for all potentials which have a quadratic dependence of $`r`$ near the origin. Our result supports therefore that free charged bosons in their condensate would show a phase transition when one varies the magnetic field. This transition would manifest itself by jumps of the magnetisation at some critical values of the magnetic field. The phenomenon would probably persist in a dilute gas of charged bosons in a neutralising background. It is however probably quite difficult to create such a jellium in the laboratory and this remains a challenging task.

## VII Acknowledgements

We thank Ph. Martin and N. Datta for some useful discussions on the Bose-Einstein condensation in the presence of a magnetic field.


text:
# A classical Over Barrier Model to compute charge exchange between ions and one–optical–electron atoms

## I Introduction

The electron capture process in collisions of slow, highly charged ions with neutral atoms and molecules is of great importance not only in basic atomic physics but also in applied fields such as fusion plasmas and astrophysics. The process under study can be written as:

$$A^{+q}+BA^{(qj)+}+B^{j+}.$$
(1)
Theoretical models are regularly developed and/or improved to solve (1) from first principles for a variety of choices of target $`A`$ and the projectile $`B`$, and their predictions are compared with the results of ever more refined experiments.
In principle, one could compute all the quantities of interest by writing the time-dependent Schrödinger equation for the system (1) and programming a computer to solve it. This task can be performed on present–days supercomputers for moderately complicated systems. Notwithstanding this, simple approximate models are still valuable: (i) they allow to get analytical estimates which are easy to adapt to particular cases; (ii) allow to get physical insight on the features of the problem by looking at the analytical formulas; (iii) finally, they can be the only tools available when the complexity of the problem overcomes the capabilities of the computers. For this reason new models are being still developed .

The present author has presented in a recent paper  a study attempting to develop a more accurate OBM by adding some quantal features. The model so developed was therefore called a semi–classical OBM. Its results showed somewhat an improvement with respect to other OBMs, but not a dramatic one.
In this paper we aim to present an OBM for dealing with one of the simplest processes (1): that between an ion and a target provided with a single active electron. Unlike the former one , this model is entirely developed within the framework of a classical model, previously studied in  (see also ), but with some important amendments and improvements which, as we shall see, allow a quite good accordance with experiments.

The paper is organized as follows: a first version of the model is presented and discussed in section II. In section III we will test our model against a first test case. From the comparison a further improvement to the model is proposed (section IV) and tested against the same case, as well as other data in section V. It will be shown that predictions with this correction are in much better agreement.

## II The model: first picture

We consider the standard scattering experiment and label T, P, and e respectively the target ion, the projectile and the electron. The system T \+ e is the initial neutral atom. Let r be the electron vector relative to T and R the internuclear vector between T and P. In the spirit of classical OBM models, all particles are considered as classical objects.
Let us consider the plane $`𝒫`$ containing all the three particles and use cylindrical polar coordinates $`(\rho ,z,\varphi )`$ to describe the position of the electron within this plane. We can arbitrarily choose to set the angle $`\varphi =0`$, and assign the $`z`$ axis to the direction along the internuclear axis.
The total energy of the electron is (atomic units will be used unless otherwise stated):

$$E=\frac{p^2}{2}+U=\frac{p^2}{2}\frac{Z_t}{\sqrt{\rho ^2+z^2}}\frac{Z_p}{\sqrt{\rho ^2+(Rz)^2}}.$$
(2)
$`Z_p`$ and $`Z_t`$ are the effective charge of the projectile and of the target seen by the electron, respectively. Notice that we are considering hydrogenlike approximations for both the target and the projectile. We assigne an effective charge $`Z_t=1`$ to the target and an effective quantum number $`n`$ to label the binding energy of the electron:$`E_n=Z_t^2/2n^2=1/2n^2`$.

As long as the electron is bound to T, we can also approximate $`E`$ as

$$E(R)=E_n\frac{Z_p}{R}.$$
(3)
This expression is used throughout all calculations in (I); however, we notice that it is asimptotically correct as long as as $`R\mathrm{}`$. In the limit of small $`R`$, instead, $`E(R)`$ must converge to a finite limit:

$$E(R)(Z_p+1)^2E_n$$
(4)
(united atom limit). For the moment we will assume that $`R`$ is sufficiently large so that eq . (3) holds, but later we will consider the limit (4), too.

On the plane $`𝒫`$ we can draw a section of the equipotential surface

$$U(z,\rho ,R)=E_n\frac{Z_p}{R}.$$
(5)
This represents the limit of the region classically allowed to the electron. When $`R\mathrm{}`$ this region is divided into two disconnected circles centered around each of the two nuclei. Initial conditions determine which of the two regions actually the electron lives in. As $`R`$ diminishes there can be eventually an instant where the two regions become connected. In fig. 1 we give an example for this.
In the spirit of OBMs it is the opening of the equipotential curve between P and T which leads to a leakage of electrons from one nucleus to another, and therefore to charge exchange. We make here the no-return hypothesis: once crossed the barrier, the electron does not return to the target. It is well justified if $`Z_p>>1`$. As we shall see just below, this hypothesis has important consequences.

It is easy to solve eq. (5) for $`R`$ by imposing a vanishing width of the opening ($`\rho _m=0`$); furthermore, by imposing also that there be an unique solution for $`z`$ in the range $`0<z<R`$:

$$R_m=\frac{(1+\sqrt{Z_p})^2Z_p}{E_n}.$$
(6)
In the region of the opening the potential $`U`$ has a saddle structure: along the internuclear axis it has a maximum at

$$z=z_0=R\frac{1}{\sqrt{Z_p}+1}$$
(7)
while this is a minimum along the orthogonal direction.

Charge exchange occurs provided the electron is able to cross this potential barrier. Let $`N_\mathrm{\Omega }`$ be the fraction of trajectories which lead to electron loss at the time $`t`$. It is clear from the discussion above that it must be function of the solid opening angle angle $`\mathrm{\Omega }`$, whose projection on the plane is the $`\pm \theta _m`$ angle. The exact expression for $`N_\mathrm{\Omega }`$ will be given below. Further, be $`W(t)`$ the probability for the electron to be still bound to the target, always at time $`t`$. Its rate of change is given by

$$dW(t)=N_\mathrm{\Omega }dt\frac{2}{T_{em}}W(t),$$
(8)
with $`T_{em}`$ the period of the electron motion along its orbit.
It is important to discuss the factor $`dt(2/T_{em})`$ since it is an important difference with (I), where just half of this value was used. The meaning of this factor is to account for the fraction of electrons which, within the time interval $`[t,t+dt]`$ reach and cross the potential saddle. In (I) it was guessed that it should be equal to $`dt/T_{em}`$, on the basis of an uniform distribution of the classical phases of the electrons. However, let us read again what the rhs of eq. (8) does mean: it says that the probability of loss is given by the total number of available electrons within the loss cone ($`W(t)\times N_\mathrm{\Omega }`$), multiplied by the fraction of electrons which reach the potential saddle. However, on the basis of the no–return hypothesis, only outgoing electrons can contribute to this term: an electron which is within the loss cone and is returning to the target from the projectile is not allowed, it should already have been captured and therefore would not be in the set $`W`$. It is clear, therefore, that the effective period is $`T_{em}/2`$, corresponding to the outgoing part of the trajectory.

A simple integration yields the leakage probability

$$\begin{array}{cc}\hfill P_l=P(+\mathrm{})& =1W(+\mathrm{})=\hfill \\ & =1\mathrm{exp}\left(\frac{2}{T_{em}}_{t_m}^{+t_m}N_\mathrm{\Omega }𝑑t\right).\hfill \end{array}$$
(9)
In order to actually integrate Eq. (9) we need to know the collision trajectory; an unperturbed straight line with $`b`$ impact parameter is assumed:

$$R=\sqrt{b^2+(vt)^2}.$$
(10)
The extrema $`\pm t_m`$ in the integral (9) are the maximal values of $`t`$ at which charge exchange can occur. If we identify this instant with the birth of the opening, using eq. (6) and (10), we find

$$t_m=\frac{\sqrt{R_m^2b^2}}{v}.$$
(11)
At this point it is necessary to give an explicit expression for $`N_\mathrm{\Omega }`$. To this end, we will consider first the case of an electron with zero angular momentum ($`l=0`$), and then will extend to nonzero values.

In absence of the projectile, the classical electron trajectories, with zero angular momentum, are ellipses squeezed onto the target nucleus. We are thus considering an electron moving essentially in one dimension. Its hamiltonian can be written as

$$\frac{p^2}{2}\frac{1}{r}=E_n.$$
(12)
The electron has a turning point at

$$r_c=\frac{1}{E_n}.$$
(13)
Obviously the approaching of the projectile modifies these trajectories. However, in order to make computations feasible, we make the following hypothesis: electron trajectories are considered as essentially unperturbed in the region between the target and the saddle point. The only trajectories which are thus allowed to escape are those whose aphelia are directed towards the opening within the solid angle whose projection on the $`𝒫`$ plane is $`\pm \theta _m`$ (see fig. 1) provided that the turning point of the electron is greater than the saddle-point distance: $`r_cz_0`$. The validity of these approximations can be questionable, particularly if we are studying the collision with highly–charged ions, which could deeply affect the electron trajectory. We limit to observe that it is necessary in order to make analytical calculations. A posteriori, we shall check the amount of error introduced by such an approximation.

The angular integration is now easily done, supposing a uniform distribution for the directions of the electrons:

$$N_\mathrm{\Omega }=\frac{1}{2}(1\mathrm{cos}\theta _m).$$
(14)
In order to give an expression for $`\theta _m`$ we notice that $`\mathrm{cos}\theta _m=z_0/(\rho _m^2+z_0^2)^{1/2}`$, with $`\rho _m`$ root of

$$E(R)=\left(\rho _m^2+\frac{R^2}{(\sqrt{Z_p}+1)^2}\right)^{1/2}+Z_p\left(\rho _m^2+\frac{Z_pR^2}{(\sqrt{Z_p}+1)^2}\right)^{1/2}.$$
(15)
It is easy to recognize that, in the right-hand side, the first term is the potential due to the electron–target interaction, and the second is the electron–projectile contribution. Eq. (15) cannot be solved analytically for $`\rho _m`$ except for the particular case $`Z_p=1`$, for which case:

$$\rho _m^2=\left(\frac{2}{E(R)}\right)^2\left(\frac{R}{2}\right)^2.$$
(16)
The form of $`E(R)`$ function of $`R`$ cannot be given analytically, even though can be quite easily computed numerically . In order to deal with expressions amenable to algebraic manipulations, we do therefore the approximation: first of all, divide the space in the two regions $`R<R_u,R>R_u`$, where $`R_u`$ is the internuclear distance at which the energy given by eq. (3) becomes comparable with its united–atom form:

$$E_n+\frac{Z_p}{R_u}=(Z_p+1)^2E_nR_u=\frac{Z_p}{(Z_p+1)^21}\frac{1}{E_n}.$$
(17)
We use then for $`E(R)`$ the united–atom form for $`R<R_u`$, and the asymptotic form otherwise:

$$\begin{array}{cc}\hfill E(R)& =E_n+\frac{Z_p}{R},R>R_u\hfill \\ & =(Z_p+1)^2E_n,R<R_u\hfill \end{array}$$
(18)
It is worthwhile explicitly rewriting eq. (16) for the two cases:

$$\begin{array}{cc}\hfill \rho _m^2& =R^2\left(\frac{4}{(E_nR+1)^2}\frac{1}{4}\right),R>R_u\hfill \\ & =\frac{1}{4}\left(\frac{1}{E_n^2}R^2\right),R<R_u\hfill \end{array}$$
(19)
and the corresponding expressions for $`N_\mathrm{\Omega }`$ are:

$$\begin{array}{cc}\hfill N_\mathrm{\Omega }=\frac{1\mathrm{cos}\theta _m}{2}& =\frac{1}{8}(3E_nR),R>R_u\hfill \\ & =\frac{1}{2}(1E_nR),R<R_u.\hfill \end{array}$$
(20)
Note that $`N_\mathrm{\Omega }=1/2`$ for $`R=0`$. This is a check on the correctness of the model, since, for symmetrical scattering at low velocity and small distances we expect the electrons to be equally shared between the two nuclei.

When $`Z_p>1`$ we have to consider two distinct limits: when $`R\mathrm{}`$ we know that eventually $`\rho _m0`$ (eq. 6). It is reasonable therefore to expand (15) in series of powers of $`\rho _m/R`$ and, retaining only terms up to second order:

$$\rho _m^2\frac{2\sqrt{Z_p}}{\left(\sqrt{Z_p}+1\right)^4}R^2\left[\left(\sqrt{Z_p}+1\right)^2Z_pE_nR\right].$$
(21)
Consistently with the limit $`R\mathrm{}`$, we have used the large–$`R`$ expression for $`E(R)`$.
The limit $`R0`$ is quite delicate to deal with: a straightforward solution of eq. (15) would give

$$\rho _m\frac{1}{(Z_p+1)E_n}+𝒪(R),$$
(22)
but calculating $`\mathrm{cos}\theta _m`$ and eventually $`N_\mathrm{\Omega }`$ from this expression gives wrong results: it is easy to work out the result $`N_\mathrm{\Omega }=1/2,R0`$. This is wrong because, obviously, the limit $`N_\mathrm{\Omega }1,Z_p\mathrm{}`$ must hold. The reason of the failure lies in the coupling of eq. (15) with the united–atom form for $`E(R)`$: one can notice that the expression thus written is perfectly simmetrical with respect to the interchange projectile–target. Because of this symmetry, electrons are forced to be equally shared between the two nuclei. This is good when dealing with symmetrical collisions, $`Z_p=Z_t=1`$, and is actually an improvement with respect to (I), where eq. (21) was used even for small $`R`$’s and one recovered the erroneous value $`N_\mathrm{\Omega }(R=0)=3/8`$. But when $`Z_p>1`$ the asymmetry must be retained in the equations. The only way we have to do this is to extend eq. (21) to small $`R`$, obtaining

$$1\mathrm{cos}\theta _m\frac{\sqrt{Z_p}}{(\sqrt{Z_p}+1)^2}\left[(\sqrt{Z_p}+1)^2Z_pE_nR\right].$$
(23)
It is straightforward to evaluate eq. (23) in the limit $`Z_p\mathrm{},R0`$, and find the sought result, 2.
We notice that, from the numerical point of view, it is not a great error using eq. (21) everywhere: the approximation it is based upon breaks down when $`R`$ is of the order of $`R_u`$ or lesser, which is quite a small range with respect to all other lengths involved when $`Z_p>1`$, while even for the case $`Z_p=1`$ it is easy to recover (see equations below) that the relative error thus introduced on $`P_l`$ is $`\mathrm{\Delta }P_l/P_l=1/24`$ for small $`b`$ (and–obviously–it is exactly null for large $`b`$). Therefore, eq. (21) could be used safely in all situations. However, we think that the rigorous altough quite lengthy derivation given above was needed since it is not satisfactory working with a model which does not comply with the very basic requirements required by the symmetries of the problem at hand.

We have now to take into account that the maximum escursion for the electron is finite. If we put $`r_c=z_0`$ and use for $`z_0`$, $`r_c`$ respectively the expressions given by (7) and (13), we obtain an equation which can be easily solved for $`R`$:

$$R=R_{}^{}{}_{m}{}^{}=(\sqrt{Z_p}+1)r_c.$$
(24)
The $`R_{}^{}{}_{m}{}^{}`$ thus computed is the maximum internuclear distance at which charge exchange is allowed under the present assumptions. Since $`R_{}^{}{}_{m}{}^{}<R_m`$ (compare the previous result with that of eq. 6 ) we have to reduce accordingly the limits in the integration in eq. (9): it must be performed between $`\pm t_{}^{}{}_{m}{}^{}`$, with the definition of $`t_{}^{}{}_{m}{}^{}`$ the same as $`t_m`$ but for the replacement $`R_mR_{}^{}{}_{m}{}^{}`$.
The result for the leakage probability is:

$$P_l=1\mathrm{exp}\left(2\frac{F(u_m)+G_Z}{T_{em}}\right),$$
(25)
where we have defined

$$\begin{array}{cc}\hfill F(u)& =\frac{\sqrt{Z_p}}{(\sqrt{Z_p}+1)^2}\left[\left((\sqrt{Z_p}+1)^2Z_p\right)\frac{b}{v}u\left(\frac{E_nb^2}{2v}\right)\left(u\sqrt{1+u^2}+\mathrm{arcsinh}(u)\right)\right],\hfill \\ \hfill G_Z& =(3F(u_u)2t_u)(Z_p=1)\hfill \\ & =0(Z_p>1),\hfill \\ \hfill u_m& =vt_m^{}/b,\hfill \\ \hfill u_u& =vt_u/b,\hfill \\ \hfill t_u& =\frac{\sqrt{R_u^2b^2}}{v}.\hfill \end{array}$$
(26)
The period can be easily computed by

$$T_{em}=2_0^{1/E_n}\frac{dr}{p}=\sqrt{2}_0^{1/E_n}\frac{dr}{\sqrt{\frac{1}{r}E_n}}=2\pi n^3$$
(27)
(this result could be found also in ).
The cross section can be finally obtained after integrating over the impact parameter (this last integration must be done numerically):

$$\sigma =2\pi _0^{b_m}bP_l(b)𝑑b.$$
(28)
Again, we have used the fact that the range of interaction is finite: the maximum allowable impact parameter $`b_m`$ is set equal to $`R_{}^{}{}_{m}{}^{}`$.

Finally, we consider the case when the angular momentum is different from zero. Now, orbits are ellipses whose minor semiaxis has finite length. We can still write the hamiltonian as function of just $`(r,p)`$:

$$\frac{p^2}{2}\frac{1}{r}+\frac{L^2}{2r^2}=E_n.$$
(29)
$`L`$ is the usual term: $`L^2=l(l+1)`$. The turning points are now

$$r_c^\pm =\frac{1\pm \sqrt{12E_nL^2}}{2E_n}.$$
(30)
and $`R_m^{}=(\sqrt{Z_p}+1)r_c^+`$.

Now the fraction of trajectories entering the loss cone is much more difficult to estimate. In principle, it can still be determined: it is equal to the fraction of ellipses which have intersection with the opening. Actual computations can be rather cumbersome. Thus, we use the following approximation, which holds for low angular momenta $`l<<n`$ (with $`n`$ principal quantum number): ellipses are approximated as straight lines (as for the $`l=0`$ case), but their turning point is correctly estimated using eq. (30). Note that also the period is modified: its correct expression is

$$T_{em}=\sqrt{2}_r^{}^{r^+}\frac{dr}{\sqrt{\frac{1}{r}E_n\frac{l(l+1)}{2r^2}}}.$$
(31)
## III A test case

As a first test case we consider the inelastic scattering $`\mathrm{Na}^++\mathrm{Na}(28\mathrm{d},29\mathrm{s})`$. We investigate this sytem since: (i) it has been studied experimentally in ; (ii) some numerical simulations using the Classical Trajectory Monte Carlo (CTMC) method have also been done on it , allowing to have detailed informations about the capture probability $`P_l`$ function of the impact parameter, and not simply integrated cross sections; (iii) finally, it has been used as test case in (I), thus allowing to assess the relative quality of the fits.

In fig. (2) we plot the normalized cross section $`\stackrel{~}{\sigma }=\sigma /n^4`$ versus  the normalized impact velocity $`\stackrel{~}{v}=vn`$ for both collisions $`nl=`$ 28d and $`nl=`$ 29s (solid line). The two curves are very close to each other, reflecting the fact that the two orbits have very similar properties: the energies of the two states differ by a very small amount, and in both cases $`E_nL^2<<1`$. The two curves show reversed with respect to experiment: $`\sigma `$(28d) it is greater than $`\sigma `$(29s). The reason is that the parameter $`r_c`$ is larger in the former case than in the latter.
We can distinguish three regions: the first is at reduced velocity around 0.2, where a steep increase of cross section appears while going towards lower velocities. Over–barrier models do not appear to fully account for this trend: they have a behaviour at low speed which is ruled approximately by the $`1/v`$ law, consequence of the straight-line impact trajectory approximation: it is well possible that this approximation too becomes unadequate in this region.
The second region covers roughly the range 0.3 $`÷`$ 1.0. Here the $`nl=`$ 29s data are rather well simulated while the present model overestimates the data for $`nl=`$ 28d. The bad agreement for $`nl=`$ 28d was already clear to Ostrovsky which attributed it to a deficiency of the model to modelize $`l`$-changing processes. It seems clear that neither our treatment of the angular momentum is sufficient to cure this defect.

Finally, there is the region at $`\stackrel{~}{v}>1`$, where again the OBM, as it stands, is not able to correctly reproduce the data. The reason for this discrepancy can be traced back to the finite velocity of the electron: the classical electron velocity is $`v_e=1/n`$, so $`\stackrel{~}{v}`$ can be given the meaning of the ratio between the projectile and the electron velocity. When $`\stackrel{~}{v}1`$ the projectile is less effective at collecting electrons in its outgoing part of the trajectory (i.e. when it has gone beyond the point of closest approach). In simple terms: an electron is slower than the projectile; when it is left behind, it cannot any longer reach and cross the potential barrier.

## IV Corrections to the model

This picture suggests a straightforward remedy: a term must be inserted in eq. (8) to account for the diminished capture efficiency. This is accomplished formally through rewriting $`N_\mathrm{\Omega }w(t,\stackrel{~}{v})N_\mathrm{\Omega }`$, with $`w1`$. We have put into evidence that $`w`$ can in principle be function of time and of the impact velocity. The simplest correction is made by assuming a perfect efficiency for $`\stackrel{~}{v}<1`$, $`w(t,\stackrel{~}{v}<1)=1`$, while, for $`\stackrel{~}{v}>1`$, no electrons can be collected after that the distance of minimum approach has been reached: $`w^+w(t>0,\stackrel{~}{v}>1)=0`$. This can appear too strong an assumption, since those electrons which are by the same side of the projectile with respect to the nucleus, and which are close to their turning point may still be captured. In fig. (2) we can compare the original data with those for $`w^+=0`$ (dashed line). The sharp variation of $`\sigma `$ at $`\stackrel{~}{v}=1`$ is obviously a consequence of the crude approximations done choosing $`w`$ which has a step–like behaviour with $`v`$.

To get further insight, we plot in fig. 3 the quantity $`bP_l(b)`$ versus $`b`$ for the collision $`\mathrm{Na}^++\mathrm{Na}(28\mathrm{d})`$. The impact velocity is $`\stackrel{~}{v}=1`$. The symbols are the CTMC results of ref. . Solid line is the model result for $`w^+=1`$; dotted line, the result for $`w^+=0`$; dashed line, an intermediate situation, with $`w^+=1/2`$. Striking features are, for all curves, the nearly perfect accordance of the value $`b3000`$ at which $`P_l=0`$ (it is $`b_m`$ according to our definition). The behaviour at small $`b`$’s ($`P_l1/2`$) is well reproduced for $`w^+=1`$ while it is slightly underestimated by the two other curves. On the other hands, only by setting $`w^+=0`$ it is possible to avoid the gross overestimate of $`P_l`$ near its maximum.

It is thus evident that the agreement is somewhat improved in the region $`\stackrel{~}{v}1`$ by letting $`w^+=0`$. However, the high–velocity behaviour is still missed by the model, which predicts a power–law behaviour $`\sigma v^1`$, while the actual exponent is higher. Within our picture, this suggests that also the capture efficiency $`w^{}=w(t<0)`$ must be a decreasing function of $`\stackrel{~}{v}`$. An accurate modelization of the processes which affect this term is difficult, and we were not able to provide it. However, some semi–qualitative arguments can be given. Let us review again the process of capture as described in section II and shown in fig. (1): if $`\stackrel{~}{v}>1`$, an electron at time $`t`$ can be in the loss cone and still not to be lost, since within a time span $`\mathrm{\Delta }t\rho _m/v`$ the position of the loss cone has shifted of such an amount that only those electrons which were closer to the saddle point than a distance $`v_e\mathrm{\Delta }t`$ could be caught. The fraction of these electrons is $`\mathrm{\Delta }t(2/T_{em})\rho _m(2/vT_{em})`$. This correction gives an additional $`1/v`$ dependence, thus now $`\sigma 1/v^2`$.
As an exercise, we try to fit experimental data using $`w`$ as a free parameter instead that a function to be determined by first principles. We choose one of the simplest functional forms:

$$w=\frac{1+|\beta |^m}{1+|\stackrel{~}{v}\beta |^m},$$
(32)
with $`\beta ,m`$ free parameters to be adjusted. This form gives the two correct limits: $`w1,\stackrel{~}{v}0`$, and $`w0,\stackrel{~}{v}\mathrm{}`$. The parameter $`\beta `$ is not really needed; it has been added to reach a better fit. Its meaning is that of a treshold velocity, at which the capture efficiency begins to diminish. In fig. (2) we plot the fit obtained with $`\beta =0.2,m=4`$ (dotted line): this is not meant to be the best fit, just a choice of parameters which gives a very good agreement with data. We see that the suggested corrections are still not enough to give the right power–law, if one needs to go to some extent beyond the region $`\stackrel{~}{v}=1`$.

## V Other comparisons

### A Iodine - Cesium collisions

We apply now our model to the process of electron capture

$$\mathrm{I}^{q+}+\mathrm{Cs}\mathrm{I}^{(q1)+}+\mathrm{Cs}^+$$
(33)
with $`q=6÷30`$. This scattering process has been studied experimentally in . It is particularly interesting to study in this context since it has revealed untractable by a number of other OBM’s, including that of (I) (for a discussion and results, see ). The impact energy is chosen equal to $`1.5\times Z_p`$ keV: since it corresponds to $`\stackrel{~}{v}<<1`$, we can safely assume $`w=1`$. The Cesium atom is in its ground state with the optical electron in a $`s`$ state.
In fig. 4 we plot the experimental points together with our estimates. In this case the fit is excellent. It is important to notice that this agreement is entirely consequence of our choice of limiting integration to $`R`$ given by eq. (24): to understand this point, observe that because of the very high charge of the projectile, the exponential term in eq. (25) is small ($`F`$, by direct inspection, is increasing with $`Z_p`$) and thus $`P_l1`$. The details of the model which are in $`F`$ are therefore of no relevance. The only surviving parameter, and that which determines $`\sigma `$, is $`R_m^{}`$. It can be checked by directly comparing our fig. 4 with fig. 1 of ref. , where results from model (I) are shown, which differ from ours just in replacing eq. (24) with eq. (6). There, the disagreement is severe.

### B Ion - Na($`n=3`$) collisions

As a final test case we present the results for collisions H–Na(3s,3p). They are part of a set of experiments as well as numerical simulations involving also other singly–charged ions: He, Ne, and Ar (see  and the references therein and in particular ; ref.  presents numerical calculations for the same system). In fig. 5 we plot the results of our model together with those of ref. . Again, we find that only by neglecting $`w^+`$ some accordance is found. The low–energy wing of the curve is strongly underestimated for Na(3s), while the agreement is somewhat better for Na(3p). Again, the slope of $`\sigma `$ for relative velocities higher than 1 could not be reproduced.
We do not show results for other ions: they can be found in fig. 3 of ref. . What is important to note is that differencies of a factor two (and even larger for 3s states) appear between light (H<sup>+</sup>, He<sup>+</sup>) and heavy (Ne<sup>+</sup>, Ar<sup>+</sup>) ions which our model is unable to predict. We can reasonably conclude therefore: (i) that the present model is not satisfactory for $`v/v_e<<1`$ (it was already pointed out in sec. IV) and for $`v/v_e>1`$ ; (ii) the structure of the projectile must be incorporated into the model otherwise different ions with the same charge should cause the same effect, at odds with experiments. As emphasized in  the energy defect $`\mathrm{\Delta }E`$ of the process is a crucial parameter: captures to states with $`\mathrm{\Delta }E0`$ are strongly preferred. Obviously, the value of $`\mathrm{\Delta }E`$ depends on the energy levels structure of the recombining ion.

## VI Summary and conclusions

We have developed in this paper a classical OBM for single charge exchange between ions and atoms. The accuracy of the model has been tested against three cases, with results going from moderate–to–good (sec. III and IV), excellent (sec. V.A), and poor–to–moderate (sec. V.B). As a rule of thumb, the model can be stated to be very well suited for collisions involving highly charged ions at low velocities.
The model is based upon a previous work , and adds to it a number of features, which we go to recall and discuss: (i) the finite excursion from the nucleus permitted to the electrons; (ii) the redefinition of the fraction of lost electrons $`dt/T_{em}dt(2/T_{em})`$; (iii) a more accurate treatment of the small impact parameter region for symmetrical collisions; (iv) the explicit-altough still somewhat approximate-treatment of the capture from $`l>0`$ states; (v) a correction to the capture probability due finite impact velocity. Let us discuss briefly each of these points:
Point (i) and (ii) contribute a major correction: in particular, (i) is essential to recover that excellent agreement found in section V.A, while (ii) accounts for the correct $`bP_l`$ behaviour at small $`b`$’s (see fig. 2).
Point (iii) is unimportant for actual computations, but corrects an inconsistency of the model.
Point (iv) has been studied in less detail, in part for the lack of experimental data on which doing comparisons.
Point (v): a good theoretical estimate of $`w`$ should be of the outmost importance for developing a really accurate model of collision at medium-to-high impact velocity. In this paper we have just attempted a step towards this direction which, however, has allowed to recover definitely better results.
Finally we recall from sec. V.B that the treatment of the projectile–or better the process of the electron-projectile binding–is an aspect which probably awaits for main improvements. We just observe that it is a shortcoming of all classical methods, that they cannot easily deal with quantized energy levels.

## Acknowledgments

It is a pleasure to thank the staff at National Institute for Fusion Science (Nagoya), and in particular Prof. H. Tawara and Dr. K. Hosaka for providing the data of ref. .