format: text

text:
# Wannier-function description of the electronic polarization and infrared absorption of high-pressure hydrogen

## I INTRODUCTION

### A Wannier functions

The electronic structure of periodic solids is usually described, in the independent-electron approximation, in terms of the extended Bloch eigenfunctions. An alternative representation is provided by the Wannier functions (WFs), which are localized, with a typical spread of the order of the atomic dimensions; they can be obtained via a unitary transformation of the Bloch states belonging to an isolated band or to a composite group of bands (i.e., bands that may be connected among themselves by degeneracies, but are separated from all others by energy gaps). For some purposes the latter description is advantageous: for instance, the WFs constructed from the states in the valence bands provide an intuitive, “chemical-like” localized picture of bonding and dielectric properties of insulators.

The major drawback of the Wannier representation is the strong nonuniqueness of the WFs: their average location, shape, and spread all depend on the arbitrary choice of gauge. In practice, this indeterminacy can be resolved by working with the set of WFs which is most localized according to some sensible criterion. A certain degree of arbitrariness still remains regarding which measure of localization to use, and in fact several alternatives have been proposed in the literature. We follow the approach of Ref. , which amounts to minimizing the sum of the quadratic spreads of the WFs (see Sec. II). We will use density-functional theory in the local density approximation (LDA), complemented by a tight-binding analysis, to investigate in terms of well-localized WFs the electronic structure and dielectric properties of compressed molecular hydrogen.

### B Compressed molecular hydrogen

Solid hydrogen under pressure has attracted considerable attention over the years, and the main goal has been to try to metallize it; this is expected to occur at high enough pressures, either by band gap closure in a molecular phase, or by molecular dissociation, whichever occurs first. However, up to the highest pressures reached so far ($`340`$ GPa), hydrogen appears to remain both molecular and insulating. Nevertheless, a rich phase diagram has emerged, with three distinct phases unambiguously identified using Raman and infrared (IR) spectroscopy.

The precise crystal structure of the high-pressure phases (phases II and III) has not been determined experimentally, and conclusive theoretical predictions have proven quite difficult, due to the quantum effects associated with the protons. The purpose of the present work is not to propose new candidate structures, but rather to make some very general points, illustrated on a couple of particularly simple prototype structures (Fig. 1), which were chosen mainly for clarity. It is hoped that, even if none of them turns out to be the correct structure of phase III (which is likely to be the case), they manage to capture some of its relevant features. For instance, ab-initio calculations at megabar pressures and low temperatures tend to favor structures in which the centers of the molecules form an hcp or, more generally, a triangular lattice (possibly with a small distortion). It is believed that in phase III the molecules are orientationally ordered, with their axes tilted away from the $`c`$-axis, as such canted structures tend to be more stable. Moreover, the resulting lowering of symmetry gives rise to IR-active vibron modes; indeed, one of the signatures of phase III, above 150 GPa, is a strong IR absorption peak in the vibron frequency range, which contrasts with the much weaker absorption found in phases I and II.

### C Organization

The paper is organized as follows. In Sec. II we briefly review the method used for constructing well-localized WFs; in Sec. III we investigate the permanent dipole moments of the “Wannier molecules” as a function of crystal structure and pressure; the results are presented in terms of a “static effective charge” vector associated with each molecule. In Sec. IV we look at the vibron-induced fluctuations of those dipoles, which can be quantified in terms of a “vibron-induced effective charge” vector on each molecule. From the effective charges the strength of the vibron IR absorption is calculated and compared with experiment, and the relative importance of the “static” and “dynamical” charge contributions is ascertained. Sec. V deals with the spatial distribution of the WFs and the effect of molecular overlap on the dielectric properties of the compressed solid. In Sec. VI we investigate the nonuniqueness associated with the definition of well-localized WFs in the dense solid. The tight-binding analysis is presented in Sec. VII. In Sec. VIII we give a discussion of our results.

Atomic units are used for all quantities except pressures, which are in gigapascal (GPa), and energies, in electron-volt (eV). Densities are expressed in terms of $`r_s`$, defined as $`(4\pi r_s^3/3)a_0^3=V/N`$, where $`N`$ is the number of protons in the volume $`V`$, and $`a_0`$ is the Bohr radius. Since the LDA tends to underestimate the pressures, in order to convert from $`r_s`$ to pressure we use an experimental equation of state extrapolated to high pressures. The LDA calculations were performed using a plane-wave cutoff of $`90`$ Ry and the bare Coulomb potential of the protons. The self-consistent calculations with 4 atoms per cell used a $`(11,11,11)`$ Monkhorst-Pack mesh for the Brillouin zone sampling. After self-consistency was achieved, the “maxloc” WFs were determined starting from the Bloch states again calculated on a $`(11,11,11)`$ mesh. In all the calculations we have used the following parameters for the $`Cmc2_1`$ structure described in Fig. 1: $`r_{\mathrm{bond}}=1.445`$ a.u., $`c/a=1.576`$, and the tilt angle $`\theta =54.0^{}`$. For $`C2/m`$ the parameters are $`r_{\mathrm{bond}}=1.456`$ a.u., $`c/a=1.588`$, and $`\theta =69.5^{}`$. In both cases the structures were obtained by minimizing the enthalpy at a fixed LDA pressure of $`100`$ GPa, with a resulting density of $`r_s=1.52`$ (which experimentally corresponds to about $`115`$ GPa, according to Ref. ). The same parameters were used at all other densities.

## II Maximally-localized Wannier functions

A set of WFs $`\{w_n(𝐫𝐑)\}`$, each labeled by a different Bravais lattice vector $`𝐑`$, can be constructed from the Bloch eigenstates $`\{\psi _{n𝐤}\}`$ in band $`n`$ using the unitary transformation

$$w_n(𝐫𝐑)=\frac{v}{8\pi ^3}_{\mathrm{BZ}}e^{i𝐤𝐑}\psi _{n𝐤}𝑑𝐤,$$
(1)
where $`v`$ is the volume of the unit cell of the crystal and the integral is over the Brillouin zone. Except for the constraint $`\psi _{n,𝐤+𝐆}=\psi _{n𝐤}`$ for all reciprocal lattice vectors $`𝐆`$, the overall phases of the Bloch functions $`\psi _{n𝐤}=e^{i𝐤𝐫}u_{n𝐤}`$ are at our disposal. However, a different choice of phases (or “gauge”),

$$u_{n𝐤}e^{i\phi _n(𝐤)}u_{n𝐤},$$
(2)
does not translate into a simple change of the overall phases of the WFs; their shape and spatial extent will in general be affected, while the location of their centers of charge remains invariant modulo a lattice vector $`𝐑`$. If the band is isolated, Eq. 2 is the only allowed type of gauge transformation for changing the WF $`w_n(𝐫)`$ associated with that band. In the case of a composite group of bands, the allowed transformations are of the more general form

$$u_{n𝐤}\underset{m}{}U_{mn}^{(𝐤)}u_{m𝐤},$$
(3)
where $`U_{mn}^{(𝐤)}`$ is a unitary matrix that mixes the bands at every wave vector $`𝐤`$. Under this transformation the individual Wannier centers can shift, but their sum over the group of bands is preserved modulo a lattice vector . Once a measure of localization has been chosen and the group of bands specified, the search for the corresponding set of “maximally-localized” WFs becomes a problem of functional minimization in the space of the matrices $`U_{mn}^{(𝐤)}`$. The strategy of Ref.  is to minimize the sum of the quadratic spreads of the Wannier probability distributions $`\{|w_n(𝐫)|^2\}`$, given by

$$\mathrm{\Omega }=\underset{n}{}\left(r^2_n𝐫_n^2\right),$$
(4)
where the sum is over the chosen group of bands (in our application they will be the valence bands), and $`𝐫_n=𝐫|w_n(𝐫)|^2𝑑𝐫`$, etc. Interestingly, the resulting “maximally-localized” (or “maxloc”) WFs turn out to be real, apart from an arbitrary overall phase.

In numerical calculations the Bloch states $`\psi _{n𝐤}`$ are computed on a regular mesh of $`k`$-points in the Brillouin zone; the integral in Eq. 1 is then replaced by a sum over the points in the mesh. In Ref.  an expression was derived for the gradient of the spread functional $`\mathrm{\Omega }`$ with respect to an infinitesimal rotation $`\delta U_{mn}^{(𝐤)}`$ of the set of Bloch orbitals, in terms of the Bloch functions in such a mesh. The only information needed for calculating the gradient are the overlaps $`u_{m𝐤}|u_{n,𝐤+𝐛}`$, where $`𝐛`$ are vectors connecting each mesh point to its near neighbors. Once the gradient is computed, the minimization can then proceed via a steepest-descent or conjugate-gradient algorithm.

Since the Bloch eigenstates at different $`k`$-points are initially computed by independent numerical matrix diagonalizations, their phases are unrelated. As a consequence, the WFs obtained directly from them using the discretized version of Eq. 1 will be poorly localized, or not localized at all. In practice the following strategy is used for preparing a better set of Bloch states as the starting point for the minimization: one chooses a set of localized “trial functions” in the unit cell, which constitute a rough initial guess at the WFs; for solid hydrogen we use Gaussians on the centers of the molecules. Then, a unitary rotation among the initial Bloch orbitals is made in order to maximize their projections onto these trial functions (the detailed procedure is described in Eqs. 62-64 of Ref. ). For a reasonable choice of the width of the Gaussians (we have used a r.m.s. width of 1 Å), the resulting rotated WFs are already extremely close to the “maxloc” ones, as discussed in Secs. VI and VII A.

## III Equilibrium Wannier dipoles

The neutral entity composed of two paired nuclei and the occupied “maxloc” Wannier orbital centered around them forms a “Wannier molecule” in the bulk of the solid. In the low-density limit the “maxloc” WFs become nonoverlapping and coincide with the ground state bonding orbitals of isolated $`\mathrm{H}_2`$ molecules; however, at the high pressures we are interested in, there is an appreciable overlap between neighboring WFs. In what follows we will sometimes loosely refer to the “maxloc” Wannier molecules in the dense solid simply as “molecules”. One should keep in mind, however, that had we chosen a different measure of localization, the resulting “maximally-localized” WFs in the dense system would in general differ somewhat from the ones we obtain. This nonuniqueness is intrinsic to WFs, and can be viewed as a manifestation of the ambiguity that always arises when trying to define “molecules” in a dense medium, either in terms of WFs or by other means. These issues will be discussed in Secs. VI and VIII B.

Fig. 2 shows a contour plot of a Wannier orbital for $`Cmc2_1`$ at $`r_s=1.52`$. The central positive contour with a large amplitude represents the molecular bond. The lowering of symmetry due to the crystalline environment is clear from the shape of the outer “corona” formed by the negative lobes, which have an antibonding character. These so-called “orthogonality tails” appear when the molecules overlap, due to the orthogonality requirement between different WFs; they are concentrated around the twelve nearest molecules, which allows for an efficient orthogonalization between neighboring WFs. We will argue in Sec. V that these overlapping orthogonality tails strongly influence the dielectric properties, and in particular the vibron IR activity.

An important effect of the anisotropic crystalline environment is that the molecules becomes polarized under the self-consistent internal electric fields inside the solid. In commonly used treatments of the dielectric properties of molecular crystals, the so-called Clausius-Mossotti approximation is assumed: the system is modeled as a sum of nonoverlapping molecular charge distributions which become polarized in the local field produced by the surrounding molecules; the bulk polarization is then the sum of the individual, nonoverlapping, molecular dipole moments, which can be straightforwardly calculated from the bulk charge density $`\rho (𝐫)`$. Such a description becomes inappropriate whenever there is significant molecular overlap: the electron density becomes different from zero everywhere, and as a result the net dipole moment becomes dependent on the particular choice made for the unit cell. In the case of molecular crystals, this is expected to occur, for instance, when the system is strongly compressed; under such circumstances a more careful treatment of the macroscopic polarization is required. According to the Berry phase theory of bulk polarization, whenever such overlap effects are significant, the macroscopic polarization $`𝐏_{\mathrm{mac}}`$ of an insulating crystal cannot be extracted from the bulk $`\rho (𝐫)`$, and is instead given, in the independent-electron approximation, by a Berry phase of the occupied Bloch states. This is a gauge-invariant quantity, and it is identical to another invariant, the total sum of the dipoles $`𝐝(n)`$ of the “Wannier molecules”:

$$𝐏_{\mathrm{mac}}=\frac{1}{v}\underset{n=1}{\overset{M}{}}𝐝(n),$$
(5)
$$𝐝(n)=2e𝐫\left|w_n(𝐫)\right|^2𝑑𝐫,$$
(6)
where $`M`$ is the number of valence bands and $`e`$ is the magnitude of the electron charge. In Eq. 6 the origin is chosen midway between the two paired protons, to cancel their contributions. The factor of 2 comes from spin-degeneracy (each occupied WF carries two electrons). We stress that Eq. 5 does not rely on the Clausius-Mossotti approximation at all, and it remains exact even when the “Wannier molecules” overlap strongly, as long as the system remains insulating. The decomposition of $`𝐏_{\mathrm{mac}}`$ into individual Wannier dipoles is a powerful analysis tool, allowing us to go beyond the Berry phase approach used in Ref. , which only gives the net polarization of the unit cell.

The permanent Wannier dipoles can be used to assign to each molecular charge distribution a “static effective charge” vector:

$$𝐪^\mathrm{s}(n)=\frac{𝐝(n)}{r_{\mathrm{bond}}(n)},$$
(7)
where $`r_{\mathrm{bond}}(n)`$ is the equilibrium bond length. This quantity, which vanishes in the low density limit of isolated molecules, measures the spontaneous symmetry-breaking charge transfer which occurs in the compressed solid whenever the two atoms in a molecule occupy crystallographically inequivalent sites. We emphasize that this definition is somewhat arbitrary, for the reasons discussed at the beginning of this section (but see Sec. VI), and therefore $`𝐪^\mathrm{s}(n)`$ does not relate directly to any measurable quantity. However, it is a sensible definition, which reduces to the natural one in the extreme ionic limit where the electron distribution is strongly concentrated around the ions. In the next section we will use it to decompose the vibron effective charge into “static” and “dynamical” contributions, with the aim of understanding the origin of the strong IR absorption in phase III.

The location of the centers of the “maxloc” WFs reflects the symmetry properties of the crystal; this is apparent from the form of the vectors $`𝐪^\mathrm{s}(n)`$, shown in the first set of columns in Tables I and II for the two structures studied. The first set of columns in Table III lists their explicit values for $`r_s=1.52`$, both with the molecules on-site and after allowing them to move away from the ideal hcp sites. For comparison we also report the values calculated using an electric quadrupolar (EQ) model. In the $`C2/m`$ structure there is a center of inversion between the two molecules in the primitive cell, leading to a cancellation of their permanent dipoles; in the lower symmetry $`Cmc2_1`$ the $`y`$-components of the individual dipoles still average to zero over the primitive cell, but the $`z`$-components add up, yielding a small spontaneous polarization along the $`c`$ axis (see Fig. 1). In general $`|q_y^\mathrm{s}|>>|q_z^\mathrm{s}|`$, so that the dipole moments make an angle with the molecular axes, i.e., $`q_{}^\mathrm{s}`$ in Table III is nonzero, although it is smaller than $`q_{}^\mathrm{s}`$. This agrees with the EQ model, where for the hcp-centered (on-site) structures the quadrupolar field at the center of the molecules is along $`y`$, so that the dipole moment along $`z`$ is solely due to the small anisotropy in the polarizability. The EQ model also predicts larger dipoles in the $`C2/m`$ than in the $`Cmc2_1`$ structure (although the effect is not nearly as pronounced as in the LDA WFs), as well as an increase in their magnitude as the molecules move off-site. Although some of the qualitative features of the LDA Wannier dipoles are captured by the EQ model, its predictions are not reliable: for instance, it does not reproduce the change in sign of $`q_z^\mathrm{s}(n)`$ for $`Cmc2_1`$ at low pressures; other discrepancies can be seen in Table III, most notably in the vibron effective charges, and are discussed in Sec. IV B.

From Fig. 2 we can already see that $`|𝐪^\mathrm{s}(n)|`$ will be small, since the Wannier distribution is fairly symmetric with respect to the center of the paired protons. The dependence of the static charges versus $`r_s`$ is plotted in Figs. 3 and 4; as expected they vanish in the low density (large $`r_s`$) limit, and even at the highest pressures ($`210`$ GPa) they are only a few percent of the electron charge, indicating that, at least in the structures under consideration, the ionicity of the molecules remains quite small, contrary to some proposals. Notice also that at high pressures $`𝐪^\mathrm{s}(n)`$ becomes quite sensitive to the crystal structure (see Table III). This is not surprising, since it is totally induced by the crystal field. At $`165`$ GPa ($`r_s=1.45`$) the permanent dipole moment of an $`\mathrm{H}_2`$ Wannier molecule in the hcp-centered $`Cmc2_1`$ structure becomes 0.075 a.u., i.e. more than $`1/10`$ of the dipole of an isolated water molecule ($`0.74`$ a.u.).

## IV  vibron infrared activity

### A Wannier-function description

The investigation of the number of IR-active lattice modes and their oscillator strengths in candidate structures is a useful guide in the search for the structures of the high-pressure phases. Here we will focus on the vibron IR absorption observed in phase III. Its relatively strong intensity is somewhat puzzling, since the stretching mode of the isolated $`\mathrm{H}_2`$ molecule is Raman-active but IR-forbidden, and this has stimulated a large number of studies.

IR absorption is caused by the coupling of light to the change in $`𝐏_{\mathrm{mac}}`$ induced by the lattice modes. In the basic theory of IR absorption in molecular crystals, the Clausius-Mossotti approximation of nonoverlapping molecules is assumed. The modern theory of polarization treats rigorously the situation where that approximation breaks down, as discussed in the previous section. There, we decomposed the spontaneous macroscopic polarization into contributions from individual Wannier molecules; here we will do the same for the vibron-induced fluctuations in $`𝐏_{\mathrm{mac}}`$. For every vibron mode $`\nu `$ we will assign to each Wannier molecule a vibron-induced “effective charge” vector (compare with Eq. 7):

$$𝐪^\nu (n)=\frac{𝐝(n)}{u_\nu },$$
(8)
where $`u_\nu `$ is the normal coordinate. The above expression is very similar to the definition of the Born effective charge associated with the stretching mode of a diatomic molecule. The vector $`𝐪^\nu (n)`$ measures the vibron-induced symmetry-breaking charge transfer; with the help of Eq. 7, it can be decomposed into two parts (see Eq. 20 of Ref. ):

$$𝐪^\nu (n)=𝐪^\mathrm{s}(n)+u_\nu \frac{𝐪^\mathrm{s}(n)}{u_\nu }.$$
(9)
Since $`𝐪^\mathrm{s}(n)`$ is the static charge, we will call the second term on the r.h.s. the dynamical charge.

Experimentally two vibrons have been detected in phase III: the lower frequency mode appears in the Raman spectrum, and the higher-frequency one in the IR spectrum. Both the $`C2/m`$ and the $`Cmc2_1`$ structures have two vibron modes: one in which the two molecules in the primitive cell vibrate in-phase $`(\nu =\mathrm{i})`$, and a higher frequency mode in which they vibrate out-of-phase $`(\nu =\mathrm{o})`$. The effective charges $`𝐪^\nu (n)`$ were calculated using Eq. 8 by changing the molecular bond length by small amounts $`\delta u_\nu `$ in the range $`[0.0015,0.0035]`$ a.u., after checking that such displacements yield essentially linear changes in the Wannier dipoles. Table III shows their values for $`r_s=1.52`$; as in the case of the static charges, we have in general that $`q_{}^\nu >q_{}^\nu `$. Figs. 3 and 4 plot the static and dynamical charges versus $`r_s`$. Since these originate from the interactions between molecules, their magnitudes vanish in the low-pressure (large $`r_s`$) limit, and increase as pressure goes up. The most striking feature is that the dynamical terms increase with pressure much more rapidly than the static ones; at the highest pressures studied they are already 3.3 to 6.7 times larger, depending on the structure and on the vibron mode. Since this appears to be a rather general feature, it is also likely to occur in the yet-undetermined structure of phase III; the observed strong vibron IR activity is probably caused by this increase of the dynamical charges. Their dominant role had been previously inferred from the strong anisotropy of the atomic Born effective charge tensors.

Figs. 3 and 4 and Table III also show that the vibron-induced fluctuations in the individual molecular dipoles ($`𝐪^\nu (n)`$), although clearly mode-dependent, are comparable for the two vibrons. (Interestingly, this is not so for the EQ model–see Section IV B.) The important difference occurs only after adding the contributions from the two molecules, and can be seen in Tables I and II. In the in-phase mode the large $`y`$-components cancel between the two molecules in the primitive cell, resulting in a weak IR activity (which actually vanishes in C2/m, since the small $`z`$-components also cancel). By contrast, in the out-of-phase mode the large $`y`$-components add up, resulting in a large net $`𝐏_{\mathrm{mac}}/u_\nu `$, and thus in a strong IR activity.

We emphasize again that all the quantities in Eq. 9 are gauge-dependent, like the Wannier functions themselves (but see Sec. VI). The gauge-invariant, measurable quantity is the net “vibron effective charge” vector, obtained by averaging $`𝐪^\nu (n)`$ over all molecules in a primitive cell:

$$\overline{𝐪}^\nu =\frac{1}{M}\underset{n=1}{\overset{M}{}}𝐪^\nu (n)=\frac{1}{n_{\mathrm{mol}}}\frac{𝐏_{\mathrm{mac}}}{u_\nu },$$
(10)
where $`n_{\mathrm{mol}}`$ is the number of molecules per unit volume; the vibron oscillator strength is proportional to $`n_{\mathrm{mol}}|\overline{𝐪}^\nu |^2`$. The calculated values of $`|\overline{𝐪}^\nu |`$ versus $`r_s`$ are plotted in Fig. 5 together with the experimental results. For the out-of-phase mode in the $`Cmc2_1`$ structure the LDA calculation yields values very close to the experimentally measured IR absorption in phase III, but on the other hand the IR activity of the in-phase mode, although weaker, would still be observed, which is not the case. As for $`C2/m`$, the IR absorption is too strong compared to experiment. Hence it seems that neither structure is likely to be the correct one for phase III (this is also supported by the large number of observed libron modes in phase III, which is incompatible with structures with such small primitive cells). Nevertheless, the above results for these structures allow us to make an important general point. They show that large permanent molecular dipoles are not required in order for strong vibron IR absorption to occur, contrary to what has been sometimes stated in the literature. In fact, in both structures the magnitude of the permanent dipoles (static charge) is far too small to account by itself for the measured absorption. However, once the dynamical charge transfer is accounted for, the resulting IR activity becomes even larger than the one measured in phase III.

Finally, Table III and Fig. 5 show that the displacement away from the hcp sites significantly increases the vibron-induced charges, as well as the static charges. This is yet another example of the strong sensitivity of the effective charges to the crystal structure, which may help explain the large difference in the intensity of IR absorption between phase III and the lower pressure phases.

### B Comparison with the EQ model

Table III also contains the values of the vibron-induced charge vectors on each molecule, as calculated from the EQ model. Like the static charges, they are significantly smaller than those obtained from the WFs. Another important difference is that in the EQ model the vibron-induced charges on the individual molecules are more than one order of magnitude smaller for the out-of-phase than for the in-phase mode, whereas their WF counterparts are comparable for the two vibrons. The reason is the following: in the EQ model the vibron-induced change in a molecular dipole can be written to first order as

$$\delta d_{}\delta \alpha _{}E_{}+\alpha _{}\delta E_{},$$
(11)
where $`\alpha `$ is the molecular polarizability and $`E`$ is the quadrupolar electric field on the molecular site (for definiteness we look at the dipole along the molecular axis; the same analysis applies to the perpendicular component). The first term on the r.h.s. is equal for the two vibrons, so that the difference between their effective charges arises from the second term. Choosing the isolated-molecule parameters from Refs.  and , it turns out that the two terms have a very similar magnitude. But whereas in the in-phase mode they have the same sign, in the out-of-phase mode they have opposite signs, so that their contributions largely cancel, resulting in a much smaller molecular effective charge. As a consequence, in the $`Cmc2_1`$ structure the in-phase oscillator strength comes out larger than the out-of-phase, which is the opposite of the LDA result. These discrepancies between the LDA WFs and the EQ model are likely to be related at least in part to the rather delocalized nature of the induced dipoles, which will be discussed in the next section, whereas the EQ model assumes point-like (infinitely localized) molecules.

## V spatial extent of the Wannier molecules

In the previous sections we focused our attention on information that can be extracted from the location of the centers of charge of the WFs. Here we will examine in detail their spatial distribution at high pressure, in particular their spatial extent. This will allow us to investigate the effects associated with the overlap between neighboring WFs; such effects are expected to be significant at megabar pressures, as suggested by the large bandwidths, in excess of 20 eV (see Fig. 12).

### A Spread of the Wannier charge and dipole distributions

The spread of the Wannier charge and dipole distributions are presented in Fig. 6. For $`Cmc2_1`$ at $`r_s=1.52`$ ($`r_s=2.0`$), a radius of 1.72 a.u. (2.25 a.u.), half the shortest intermolecular distance $`d_{\mathrm{H}_2\mathrm{H}_2}`$, encloses about 72% (85%) of the charge and only 17% (28%) of the $`y`$-component of the dipole. This suggests that already at $`r_s=2.0`$ ($`13`$ GPa) the overlap between nearby WFs is far from negligible. It is also clear that the dipole is significantly more spread out than the charge, with very large contributions arising from the orthogonality tails in the overlap region, where the Wannier charge density is very small. The longer range of the dipole is to be expected, due to the factor $`𝐫`$ in Eq. 6, but the large magnitude of the effect is somewhat surprising. Fig. 7 shows even more clearly that for both pressures a rather small fraction of the total charge, located in the overlapping tails of the WF, is responsible for most of the dipole. Also striking is the fact that, at $`r_s=1.52`$, up to a radius of 7 a.u. the accumulated $`d_z(1)`$ remains positive, whereas the converged value is negative; this suggests that the agreement in sign with the EQ model (see Table III) may be fortuitous, since in that model the dipole is caused by the electric field at the center of the molecule.

By analogy with the radially integrated Wannier charge and dipole distributions (Fig. 6), one can plot the derivative of these quantities with respect to the normal coordinate of a vibron mode (Fig. 8). At high pressures the contributions from the overlapping tails to the change in the dipole moment are very significant, even more so than for the equilibrium dipole. In other words, in the dense solid the dynamical charge transfer processes responsible for the IR activity are very delocalized.

The effect of a vibron on the Wannier charge distribution is depicted in the upper panel of Fig. 9: charge is depleted from the inner part and accumulates in the outer part of the molecule. Note that the charge transfer occurs mainly along the molecular axis, and is essentially symmetrical with respect to the molecular center, as one would expect from stretching an isolated molecule, which does not break the symmetry between the two atoms. This kind of charge transfer alone would lead to a zero net change in the (vanishing) molecular dipole, and hence to no IR absorption. The contribution of a Wannier molecule to the IR activity of the crystal comes from the comparatively small odd part (with respect to its center) of the charge transfer. Since this is barely visible in the upper panel of Fig. 9, in the lower panel we have removed the large even part. It is interesting to note that near the paired protons most of the odd part is oriented roughly perpendicularly to the molecular axis, making a small angle with the $`c`$ axis, such that it gives a small positive contribution to $`q_y^\mathrm{o}(1)`$ and a larger negative contribution to $`q_z^\mathrm{o}(1)`$. This observation is supported by Fig. 8, which shows that for small radius the accumulated $`q_z^\mathrm{o}(1)`$ is negative and larger than the accumulated $`q_y^\mathrm{o}(1)`$, which is positive. For large radius $`q_y^\mathrm{o}(1)`$ changes sign and ends up overtaking $`q_z^\mathrm{o}(1)`$, and the net molecular vibron charge vector has a larger projection along the molecular axis than perpendicularly to it (see Table III).

The results of this section should be relevant for models that attempt to account for the dielectric properties of compressed hydrogen. For instance, it seems unlikely that models based on point-like objects, such as the electric quadrupole (EQ) model, contain all the important ingredients that lead to the strong IR absorption in the highly compressed phase III. The reasons are twofold: at such high densities (i) the electrostatic interactions are expected to differ substantially from the ideal quadrupolar one (and in fact the validity of a multipole expansion becomes questionable when the molecular charges overlap significantly), and (ii) the “classical” treatment of polarization, based on the bulk $`\rho (𝐫)`$, becomes inadequate. We note that although the EQ model can account for both static and dynamical charges (see Sec. IV B), both effects are then due to local fields and polarizabilities (see also Fig. 2 of Ref.  and associated discussion regarding local versus nonlocal mechanisms).

The central conclusion of the preceding analysis is that the contributions to the induced molecular dipoles (and their fluctuations) arising from the overlapping tails of the Wannier orbitals, which extend well beyond the nearest neighbor molecules, are crucial. It is instructive to contrast this state of affairs with what happens in liquid water: there, the contribution from the orthogonality tails is negligible. This difference may stem from the fact that an isolated water molecule is polar, so that the effect of the liquid environment is only to modify a previously existent dipole moment, whereas in solid hydrogen the molecular dipole is totally induced. Induced dipoles tend to be rather extended because the outer regions of the molecules are most easily polarizable; thus their contribution to the dipole can be large, even though $`\rho (𝐫)`$ is small, because of the $`𝐫`$ factor in Eq. 6. In conclusion, the WF analysis strongly suggests that the Clausius-Mossotti picture of nonoverlapping dipoles breaks down rather dramatically for solid hydrogen at megabar pressures.

### B Measuring the molecular overlap

The overlap between the charge distributions of neighboring WFs can be quantified as

$$O_{mn}=\frac{\left|w_m(𝐫)\right|^2\left|w_n(𝐫)\right|^2𝑑𝐫}{\left(\left|w_m(𝐫)\right|^4𝑑𝐫\right)^{1/2}\left(\left|w_n(𝐫)\right|^4𝑑𝐫\right)^{1/2}}.$$
(12)
For $`Cmc2_1`$ the largest value of $`O_{mn}`$ is 0.005 at $`r_s=2.0`$ and 0.021 at $`r_s=1.52`$; the latter value is still quite small, roughly twice the value for WFs located on nearby water molecules in liquid water. Thus, by inspection of $`O_{mn}`$ alone one would not suspect that the overlapping tails are so much more important for the dipole moments in compressed solid hydrogen than in liquid water. The reason is that $`O_{mn}`$ measures the overlap between charge distributions, whereas in this system the dipoles are much more spread out.

Another indication that overlap effects are important comes from the well-known fact that a very large number of $`k`$-points is required to converge the total-energy calculations in compressed hydrogen. In fact, if the molecules were strictly nonoverlapping a single $`k`$-point would suffice for computing all physical properties. In Fig. 10 is shown the static charge (i.e., the dipole) for WFs obtained using different meshes of $`k`$-points. It is clear that a dense mesh is required for converging this quantity. This results from the fact that, when using a discrete mesh, the WFs are actually periodic in real space, with a periodicity which is inversely proportional to the spacing between neighboring points. Hence the need for a fine sampling of the Brillouin zone is just a manifestation of the large contributions to the Wannier dipole arising from the tails far away from the “home” unit cell. In solid hydrogen a dense mesh of $`k`$-points is expected to be even more important for the dielectric properties than for the total energy, since $`𝐏_{\mathrm{mac}}`$ is particularly sensitive to the Wannier tails.

### C Quadratic spread and localization length

Another way of quantifying the spatial extent of the WFs is in terms of their quadratic spread. According to Eq. 4, the r.m.s. width of the Wannier probability distribution, averaged over all three Cartesian directions and over the occupied WFs, is $`\overline{\lambda }=(\mathrm{\Omega }/3M)^{1/2}`$. This quantity is plotted versus $`r_s`$ in Fig. 11. Notice that it increases with increasing pressure (decreasing $`r_s`$), i.e., the Wannier molecules become more extended upon compression, which is the opposite of what happens in the usual models of molecular solids. This is an overlap effect, caused by the orthogonality requirement. It is due to the enhancement of the outer corona shown in Fig. 2, which is not included in the definition of “molecules” in those models. It can also be viewed as a result of the gap reduction with pressure.

Also plotted in Fig. 11 is the electronic localization length $`\overline{\xi }=(\mathrm{\Omega }_\mathrm{I}/3M)^{1/2}`$, where $`\mathrm{\Omega }_\mathrm{I}<\mathrm{\Omega }`$ is the gauge-invariant part of the spread of the WFs. $`\overline{\xi }^2`$ measures the mean-square quantum fluctuation of the macroscopic polarization, normalized in such a way as to be finite for insulators, and diverging when the band gap closes. Notice that in the low-density limit $`\overline{\lambda }\overline{\xi }`$; that happens because there is only one occupied WF per molecule, and can be understood by comparing Eqs. 14 and 15 of Ref. . As expected, at high pressures $`\overline{\lambda }`$ and $`\overline{\xi }`$ are larger for the $`C2/m`$ structure, which has the smaller band gap: at $`r_s=1.52`$ $`\overline{\lambda }`$ has increased by 13% (20%) in $`Cmc2_1`$ ($`C2/m`$) with respect to the low-density (isolated molecule) value.

## VI Uniqueness of well-localized Wannier molecules

As discussed earlier, the WFs are strongly nonunique; in particular, it is only the sum of all the Wannier dipoles over a primitive cell which is physically meaningful (Eq. 5), whereas the individual dipoles (and hence $`𝐪^\mathrm{s}(n)`$ and $`𝐪^\nu (n)`$) are gauge-dependent. Nevertheless, in this work we have been looking at the individual dipoles of the “maxloc” WFs in an attempt to extract from them useful physical information. The underlying assumption is that in practice well-localized WFs are fairly unique. For solid hydrogen this is obviously true in the low density limit, where they reduce to the bonding orbitals of the individual molecules. Here we will discuss to what extent that assumption holds for the compressed solid as well.

A systematic way of assessing the degree of uniqueness of well-localized WFs would be to implement different localization criteria and then compare the resulting “maximally-localized” WFs. We have not attempted such a detailed study; instead we have performed a simpler test, which is a first step in that direction. As mentioned in Sec. II, in the method we are using an initial guess is made for the localized WFs, with the help of “trial functions”, which in our case are bond-centered Gaussians (we will call the resulting orbitals “projected WFs”). Their localization is then enhanced by minimizing the quadratic spread $`\mathrm{\Omega }`$ (Eq. 4), yielding the “maxloc” WFs. Since the projected WFs are totally oblivious of the localization criterion that one later uses to further localize them, it seems reasonable to assume that the difference between the projected and the “maxloc” WFs is an upper bound to the differences that would occur between WFs obtained using any two “sensible” localization criteria.

Let us consider the $`Cmc2_1`$ structure at $`r_s=1.52`$, for which the average r.m.s. width of the “maxloc” WFs is $`\overline{\lambda }=1.11`$ a.u. If we choose the r.m.s. width of the initial bond-centered Gaussians to be $`1.89`$ a.u. (1 Å), the resulting projected WFs are essentially indistinguishable from the “maxloc” ones: for instance, the curves corresponding to those in Fig. 6 are virtually identical, and the individual Wannier dipoles remain the same to at least six significant digits! This is compelling evidence for a high degree of uniqueness of well-localized WFs in this system, at least for the high-symmetry configurations that we studied. If we double the width of the initial Gaussian, some differences start to appear. They remain barely visible in the accumulated radial integral of the probability (upper panel of Fig. 6) but are noticeable, although still relatively small, in the radially integrated dipole (lower panel of Fig. 6). For instance, the large $`y`$-component of the dipole changes by around 2%.

## VII Tight-binding analysis

In this section, we investigate whether the essential physics of the WFs in high-pressure H<sub>2</sub> phases can be captured by a simpler tight-binding (TB) approach. As is well known, the TB approximation provides a simple, computationally inexpensive method for computing electronic structure effects, and has the additional advantage that its output is easily interpreted in terms of a local, real-space picture. Thus, TB is a natural approach to explore here, where we want to study the dielectric properties of H<sub>2</sub> phases from just this kind of local point of view. We confirm below that a previously-proposed $`sp^3`$ TB model  provides a good description of the occupied bands in these systems, and show how the operations of constructing the WFs and computing their contributions to dielectric properties (such as electric polarization) can be carried out in the TB framework. Finally, using the fact that the TB representation automatically provides an atom-by-atom and orbital-by-orbital decomposition, we obtain useful insights into the nature of the WFs and their contributions to the dielectric properties.

### A Tight-binding formalism

In the TB method, the Bloch functions $`\psi _{n𝐤}`$ are expanded in a basis of atomic-like orbitals $`\varphi _{il}`$ as

$$\psi _{n𝐤}(𝐫)=\underset{il}{}C_{n𝐤}(il)e^{i𝐤𝐫_{il}}\varphi _{il}(𝐫).$$
(13)
Here $`l`$ labels the unit cell located at $`𝐑_l`$, $`i`$ labels an orbital on the atom at $`𝐫_{il}=𝐑_l+\tau _i`$ (where $`\tau _i`$ specifies the relative position of the atom within the unit cell), and the vector of coefficients $`C_{n𝐤}(il)`$ forms the TB representation of the Bloch function. Our goal is to carry out a unitary transformation to a set of $`M`$ localized WFs

$$w_\alpha (𝐫)=\frac{1}{N}\underset{n𝐤}{}U_{\alpha n}^{(𝐤)}\psi _{n𝐤}(𝐫)$$
(14)
associated with a set of $`M`$ occupied bands, where the $`U_{\alpha n}^{(𝐤)}`$ are $`𝐤`$-dependent $`M\times M`$ unitary matrices that will be fixed by the requirement of maximal localization . Introducing the TB representation of the WF

$$w_\alpha (𝐫)=\underset{il}{}W_\alpha (il)\varphi _{il}(𝐫),$$
(15)
it follows that

$$W_\alpha (il)=\frac{1}{N}\underset{n𝐤}{}U_{\alpha n}^{(𝐤)}e^{i𝐤𝐫_{il}}C_{n𝐤}(il).$$
(16)
The essential ingredients needed for the construction of the maximally-localized WFs , or for the computation of the Berry-phase polarization , are inner products $`u_{n𝐤}|u_{n^{}𝐤^{}}`$ between the cell-periodic part of the Bloch functions

$$u_{n𝐤}(𝐫)=e^{i𝐤𝐫}\psi _{n𝐤}(𝐫)$$
(17)
at nearby $`k`$-points in the Brillouin zone. In principle, the calculation of the $`u_{n𝐤}|u_{n^{}𝐤^{}}`$ requires a detailed knowledge of the basis orbitals $`\varphi _{il}`$, which has been done in Ref. . However, in the spirit of minimal empirical TB, we make the approximation

$$u_{n𝐤}|u_{n^{}𝐤^{}}=\underset{il}{}C_{n𝐤}^{}(il)C_{n^{}𝐤^{}}^{}(il).$$
(18)
\[When completing the circuit across a Brillouin zone boundary, the relation $`C_{n,𝐤+𝐆}(il)=e^{i𝐆𝐫_{il}}C_{n𝐤}(il)`$ should be used to translate by a reciprocal lattice vector $`𝐆`$.\]

Eq. (18) can be derived via a Taylor expansion of the exponential factor $`\mathrm{exp}[i𝐤(𝐫𝐫_{il})]`$, with the following assumptions: (i) that the TB basis orbitals are orthonormal, $`\varphi _{il}|\varphi _{i^{}l^{}}=\delta _{ii^{}}\delta _{ll^{}}`$; (ii) that the position operator is diagonal in the TB basis, $`\varphi _{il}|𝐫|\varphi _{i^{}l^{}}=𝐫_{il}\delta _{ii^{}}\delta _{ll^{}}`$; and (iii) that matrix elements of higher powers of the position operator are likewise trivial,

$$\varphi _{il}|x^py^qz^r|\varphi _{i^{}l^{}}=x_{il}^py_{il}^qz_{il}^r\delta _{ii^{}}\delta _{ll^{}}.$$
(19)
Conditions (i) and (ii) are actually special cases of (iii) and all all are quite artificial in that they cannot be satisfied for actual basis functions. For example, while an $`sp^3`$ hybrid on a given atom should have its charge center displaced from the geometric center of the atom, condition (ii) does not allow this effect to be captured. Similarly, the spread $`\varphi _{il}|r^2|\varphi _{il}\varphi _{il}|𝐫|\varphi _{il}^2`$ of an individual basis orbital is taken to vanish, according to condition (iii). Nevertheless, Eq. (18) is the logical extension of the empirical TB philosophy, in which one tries to avoid introducing any additional parameters beyond those needed to parameterize the Hamiltonian itself. Despite its simplicity, the tests presented below demonstrate that this approach captures much of the interesting complexity of the WFs, at least for the systems under study here. A similar approximation was previously shown to allow for reasonably accurate TB calculations of dynamical effective charges in semiconductors .

In practice, we work entirely within the TB representation. First the $`C_{n𝐤}(il)`$ are determined on a regular mesh of $`k`$-points by solving the standard secular equation involving the Hamiltonian matrix $`H_{ii^{}}^{(𝐤)}`$. Then the electric polarization can be computed by inserting Eq. (18) into the formalism of Ref. . Similarly, an “optimal” set of unitary matrices $`U_{\alpha n}^{(𝐤)}`$ can be obtained by inserting Eq. (18) into the formalism of Ref. , and from these, the WFs $`W_\alpha (il)`$ obtained via Eq. (16). The resulting WFs are optimal in the sense of being maximally localized in real space, i.e., of minimizing Eq. (1).

In the “maxloc” method , one usually begins by choosing a set of localized “trial functions,” and making a preliminary unitary rotation among the Bloch orbitals in order to maximize their projections onto these trial functions, as discussed in the LDA context at the end of Sec. II. In the TB context, we have found the following natural way of constructing the trial Wannier functions. Since we have two molecules per cell, we want to carry out the $`2\times 2`$ rotation that makes one state have most of its projection on the first molecule, and the other have most of its projection on the second molecule. To do this, we consider the difference $`\mathrm{\Delta }\text{P}=\text{P}_1\text{P}_2`$ of projection operators $`\text{P}_1`$ and $`\text{P}_2`$ onto the first and second molecule, respectively. (In the TB basis, $`\mathrm{\Delta }𝐏`$ is just a diagonal matrix with $`\pm 1`$ diagonal entries.) Then, at each $`𝐤`$, we diagonalize $`\mathrm{\Delta }𝐏`$ in the space of the two Bloch states, and set the phase of each eigenvector by requiring that its inner product with an even linear combination of $`s`$ orbitals on the two atoms comprising the molecule should be real and positive. We find that the unitary transformation $`\stackrel{~}{U}_{\alpha n}^{(𝐤)}`$ obtained in this way turns out to be an excellent approximation to the “maxloc” $`U_{\alpha n}^{(𝐤)}`$ which minimizes Eq. (1). In fact, subsequent minimization typically only leads to changes of WF coefficients of order one part in $`10^5`$, and so that in practice it is not even necessary to carry out the maxloc minimization procedure. This is consistent with our similar experience in the LDA context as discussed at the end of Sec. VI.

### B Details of the tight-binding model

The TB parameterization we used is the one proposed by Chacham et al. . These authors showed that the main characteristics of the electronic structure of high-density solid hydrogen at megabar pressures could be reproduced by using a minimal orthogonal TB basis comprised of $`s`$, $`p_x`$, $`p_y`$, and $`p_z`$ orbitals on each hydrogen atom. The intermolecular matrix elements are taken as

$`V_{ss}`$ $`=`$ $`V_{ss}(d_0)e^{\alpha (1d/d_0)}`$ (20)
$`V_{pp\sigma }`$ $`=`$ $`V_{pp\sigma }(d_0)e^{\beta (1d/d_0)}`$ (21)
$`V_{sp}`$ $`=`$ $`V_{sp}(d_0)e^{[(\alpha +\beta )/2](1d/d_0)}`$ (22)
$`V_{pp\pi }`$ $`=`$ $`0`$ (23)
with dimensionless constants $`\alpha =5.76`$ and $`\beta =2.52`$, where $`d`$ is the interatomic distance and $`d_0=3.79`$Å is the equilibrium hcp lattice constant at zero pressure. Given in Table IV are the $`V_{ss}(d_0)`$, $`V_{pp\sigma }(d_0)`$, $`V_{sp}(d_0)`$, and the intramolecular matrix elements $`V_{ss}`$, $`V_{sp}`$ and $`V_{pp\sigma }`$ which are independent of $`d`$. The intra-atomic parameter $`ϵ_sϵ_p=20`$eV.

In addition to this original tight-binding scheme, we have also tested an extended scheme that includes a correction designed to incorporate the effects of the quadrupolar electrostatic fields arising from neighboring molecules. In this extended “TB+Q” scheme, the molecules are first modeled as point quadrupoles centered at the molecular sites (mid-bond positions). The quadrupole moment tensor for each molecule is taken from the free-molecule calculations of Ref.  by assuming a linear dependence upon the bond length in the range of 1.4-1.6 a.u. The total quadrupolar electric field is then evaluated at each molecular site, and the electrostatic potential shift on each atom in the molecule is calculated by assuming a linear extrapolation to the atomic position. Finally, the diagonal elements (self-energies) of the TB Hamiltonian matrix are modified by adding these energy shifts, and the solution of the secular equation then proceeds as usual.

### C Band structure

We have applied the TB model to the same $`Cmc2_1`$ and $`C2/m`$ candidate H<sub>2</sub> structures studied with LDA methods in earlier sections, and confirmed that this TB model does a good job of reproducing the critical features of the electronic band structure. For convenience, we present only results on the $`Cmc2_1`$ geometry; the corresponding results for $`C2/m`$ are qualitatively similar. As already indicated in Sec. I C, our $`Cmc2_1`$ structure has $`r_s=1.52`$ a.u., $`r_{\mathrm{bond}}=1.445`$ a.u, $`c/a=1.576`$, and the tilt angle $`\theta =54.0^{}`$. The two-molecule (four-atom) unit cell is illustrated in Fig. 1. Use of the $`sp^3`$ TB basis leads to a $`16\times 16`$ TB Hamiltonian matrix.

Figure 12 shows the good agreement between TB and LDA band structures for this geometry. The agreement in the occupied valence-band region (lowest two bands) is excellent, and the resemblance in the conduction-band region is also reasonable. The TB model predicts a gap closure at a density of 0.3962 mol/cm<sup>3</sup>, consistent with the results from other studies . The band structure is hardly affected at all if the TB+Q theory is used in place of simple TB.

The WFs are constructed for the two occupied bands by using the definition in Eq. 16. A $`10\times 10\times 10`$ $`k`$-point mesh is used in calculation. Since there are only two WFs in the unit cell, and these are related to each other by a symmetry, it suffices to analyze just one of them. In Table V we analyze the spatial distribution of the WF by decomposing into contributions coming from the “home molecule” and the first two nearest-neighbor shells of molecules in real space. The home molecule is labeled as “Neighbor 0”, the next six neighboring molecules form the first shell at a radius of $`3.4342`$ a.u. ($`0.9767`$ in units of lattice constant), and the second shell comprises the next six neighboring molecules at a radius of $`3.5162`$ a.u. (one lattice constant). Table V gives a clear picture of the spatial structure of the WFs. We can see that the WFs have about 85% of their probability on the home molecule, 91% (cumulatively) inside the first shell, and 97% up to the second shell. In the TB framework, the contribution to the WFs can be very easily decomposed further into the $`sp_3`$ tight-binding basis orbitals, as shown in Fig. 13. While the $`s`$-orbital contribution is $``$5 times larger than that of the $`p`$ orbitals, we find that the $`s`$ contribution is almost entirely localized to the home molecule. On the other hand, although the $`p`$ orbitals give a smaller total contribution, they play a much bigger role in the tail region that determines the spatial distribution of the WFs.

### D Tight-binding Wannier functions

The above quantities have been recomputed using the TB+Q theory in place of the simple TB theory, but the differences are not significant.

According to the modern understanding , the electronic contribution to the polarization can be equivalently expressed either in terms of a Berry phase (BP) computed from the Bloch functions, or in terms of the displacements of the Wannier centers (i.e., from the molecular dipole moments). These two approaches are compared for each of the three different computational schemes (TB, TB+Q, and LDA) in Table VI. It can be seen that the bulk polarizations $`P_z`$ computed from the WF and BP approaches are in good agreement with each other ($``$5%) for all three schemes, indicating good internal consistency. (The small discrepancies can be traced mainly to incomplete $`k`$-point convergence.)

The advantage of the WF approach is that the decomposition into molecular dipole moments can give some insight into the microscopic origins of the dielectric properties. As in the LDA calculations, we find larger $`d_y`$ components and smaller $`d_z`$ components, with a pattern of signs determined by symmetry requirements. The interpretation of the vibron infrared activities in terms of this picture has already been discussed in Sec. IV.

Unfortunately, the level of agreement between the TB and LDA results for the molecular dipole moments is somewhat disappointing. We find a TB $`d_y`$ value that has the right sign, and the correct order of magnitude, relative to the LDA value, but the actual values differ by a factor of about two. The TB $`d_z`$ value even has the wrong sign, but this is related to the fact that the LDA $`d_z`$ value happens to come out very small ($``$10 times smaller than for $`d_y`$; see also Fig. 6). Thus, it is not surprising that the relative TB error in $`d_z`$ is large, even though the absolute TB error is actually smaller for $`d_z`$ than for $`d_y`$. Because the simple TB theory does not include any charge self-consistency, it was hoped that the extension to the TB+Q theory might improve the results by incorporating a leading (quadrupolar) Coulomb contribution. Table VI shows that the changes from TB to TB+Q are in the right direction, and there is some improvement in the $`d_z`$ (and therefore $`P_z`$) values, but the $`d_y`$ discrepancy is hardly affected.

Thus, while the TB theory gives a picture that is qualitatively correct, it is clear that there is room for improvement. Possible avenues for future investigation may be to consider nonorthonormal TB basis sets, to parameterize and include off-diagonal matrix elements of the position operators between TB basis functions, to include $`s^{}`$ or $`d`$ orbitals in the TB basis, or to include other Coulomb contributions (e.g., a self-consistent inclusion of the field arising from the induced dipoles).

Finally, we calculated the shell-by-shell spatial decomposition of the molecular dipole moments in the TB scheme. The results appear in Fig. 14. (By symmetry, all $`d_x`$ contributions are identically zero.) Comparing Fig. 14 with Fig. 13, one sees that although the contribution up to the first neighbor shell is 91% for the density, it is only 62% and 60% for $`d_y`$ and $`d_z`$, respectively. Up to the second shell, the contribution is 97% for the density and 95% for $`d_y`$, but only 67% for $`d_z`$. (Of course, $`d_z`$ is the only component that survives in the summation giving the bulk polarization.) Thus, the dipoles are found to be very delocalized, spanning over quite a few neighboring molecules, in agreement with the LDA results.

To summarize this section, we have demonstrated that an empirical TB framework allows for a very useful qualitative (and often semiquantitative) analysis of the WFs in systems such as the H<sub>2</sub> phases under study here. It is typical of the empirical TB approach that one cannot insist on quantitative accuracy at the level of first-principles schemes. However, the TB approximation has proven enormously useful over the years because of its simplicity, transparency, and ease of application. These features often allow for insightful modeling of simple systems, or for efficient calculations of large and complex systems where ab-initio schemes would not be practical. For example, one could easily use the present scheme for a computationally efficient analysis of the local dielectric structure of more complex H<sub>2</sub> crystal structures  or of supercell realizations of disordered H<sub>2</sub> systems. We expect that the coupling of Wannier and TB approaches will prove to be a useful strategy in a wide variety of other materials systems as well.

## VIII discussion

### A Wannier functions and the Clausius-Mossotti approximation

As discussed by Harrison, WFs provide a rigorous formulation of the “extended bond orbitals” which appear in tight-binding models. It should be noted that orbitals more localized than WFs can be constructed by removing the orthogonality constraint, which gives rise to the long-range tails; the resulting orbitals correspond to the “bond orbitals” in Harrison’s picture. However, the straightforward connection to the polarization is then lost, since $`𝐏_{\mathrm{mac}}`$ is no longer simply the sum of the dipole moments of those orbitals, and additional cross-terms connecting orbitals in different cells appear. Therefore, as far as dielectric properties are concerned, the “maxloc” WFs seem to provided the most localized description of the electronic structure. In particular, the spatial extent of their dipole moments (Figs. 6 and 14) provides a natural length scale to compare with the intermolecular distance in order to assess the validity of the Clausius-Mossotti approximation of nonoverlapping dipoles.

### B Identifying “molecules” in the dense solid

There is a vast literature dealing with useful ways of identifying individual “atoms” inside a molecule, or individual “molecules” in a dense medium (see Refs.  and references cited therein). In this paper we have advocated using “maxloc” WFs as a useful computational definition of “$`\mathrm{H}_2`$ molecules” in the insulating molecular solid (similarly, in Refs.  WFs were used for defining “water molecules” in liquid water). In spite of having some counterintuitive features (becoming larger under pressure) and some conceptual limitations (being defined only in the independent-electron framework; not being totally unique, since they depend on the measure of localization), they have the following important conceptual advantages. (i) When the Clausius-Mossotti approximation breaks down - which is when the ambiguity in identifying individual “molecules” appears - the sum of the dipoles of the Wannier molecules still gives the bulk polarization exactly (Eq. 5). By contrast, any definition of “molecules” based upon a direct partition of the bulk electronic charge density necessarily yields an incorrect result, since away from the Clausius-Mossotti limit the information about the bulk $`𝐏_{\mathrm{mac}}`$ is not in $`\rho (𝐫)`$. For example, in Ref.  it was found that different schemes for partitioning the charge density yield very different molecular dipoles in liquid water. We would expect the situation with such approaches to be even more severe in the case of solid hydrogen, since the dipoles are smaller and originate mostly in the Wannier tails. (ii) At high densities the WFs interpenetrate one another, so that the charge density at a given point is a sum of contributions from different molecules. Thus, effects related to molecular overlap are naturally discussed in the Wannier representation.

### C Intramolecular versus intermolecular charge transfer

There has been some debate about whether the strong vibron IR activity in compressed solid hydrogen is due mainly to intramolecular or intermolecular charge transfer. As a result of the ambiguity in defining “molecules” in the solid, the question is to some extent ill-posed. Of course, there is no charge transfer between “Wannier molecules”, since their charge is fixed. However, a heuristic argument can be attempted: the fraction of the vibron-induced molecular effective charge originating in the central part of the WF can be viewed as the “intramolecular” contribution (polarization of the molecular bond), whereas the contributions from the orthogonality tails in the overlap regions are of “intermolecular” origin. As the pressure increases, the “intramolecular” part of the WF becomes smaller and contains less charge (and presumably becomes less polarizable), whereas the opposite happens to the outer corona. The results of Secs. V and VII suggest that as far as polarization–and hence IR activity–are concerned, at megabar pressures the “intermolecular” contribution associated with the outer corona is dominant.

### D Summary

Using a method for computing well-localized Wannier functions, we have presented a “chemical-like” localized picture of the electronic structure of solid molecular hydrogen, and used it to investigate the dielectric properties of the compressed system. This approach is particularly well-suited for studying the effects of molecular overlap, which become increasingly more important at high pressures. We found the somewhat surprising result that already at moderate pressures the orthogonality tails of the WFs in the overlap regions give rise to most of the induced dipole moments on the “Wannier molecules”; this clearly indicates a breakdown of the Clausius-Mossotti approximation. Under those circumstances the electric polarization cannot be extracted from the electronic charge density in the unit cell, and the Berry-phase/Wannier-function theory must be used instead. The present approach clarifies the origin of the strong vibron IR activity in phase III and identifies the dominant mechanism: even though the permanent dipoles of the molecules in our prototype structures are too small to account for the vibron oscillator strength, the vibron-induced dipole fluctuations are of the right order of magnitude in $`Cmc2_1`$, and actually too large in $`C2/m`$. In other words, in the strongly compressed solid the dynamical contribution to the vibron effective charge dominates the static one. This conclusion seems to be supported by the fact that, even though several libron modes have been identified in phase III, no strong libron IR activity has been reported. If the molecules had any significant spontaneous polarization, it should manifest itself in the librational IR absorption. Thus, we see that well-localized Wannier functions provide a useful definition of “$`\mathrm{H}_2`$ molecules” in the dense solid, which can be used to gain important insight into the microscopic mechanisms of its dielectric response.

## Acknowledgments

I.S. and R.M.M. acknowledge financial support from DOE Grant No. DEFG02-96-ER45439, X.Z. and D.V. from NSF Grant DMR-9981193. I.S. acknowledges financial support from FCT (Portugal).


text:
# Lattice QCD Calculations of Hadron Structure: Constituent Quarks and Chiral Symmetry

## 1 INTRODUCTION

Previous experimental reports  of the strange-quark contribution to the nucleon magnetic moment, $`G_M^s(0)`$, have involved the poorly understood nucleon anapole moment. However, combining the original proton data with new data from parity violation experiments on the deuteron now allow an independent isolation of $`G_M^s(0)`$. The preliminary analysis  indicates $`G_M^s(0)`$ is smaller than the earlier result  of

$$G_M^s(0.1\mathrm{GeV}^2)=+0.61\pm 0.17\pm 0.21\mu _N,$$
(1)
but positive values for $`G_M^s(0)`$ are still favoured.

In this report, we elucidate the impact of this experimental result on our understanding of the concept of a constituent quark. In particular, the constituent quark approach places properties such as the magnetic moment contribution of a quark to a baryon intrinsic to the constituent quark itself and therefore not dependent upon the environment in which the quark resides.

The extraction of $`G_M^s(0)`$ is based on the approximation of charge symmetry in the nucleon moments; i.e. equal $`u`$ and $`d`$ current quark masses and negligible electromagnetic corrections. This approximation is expected to be valid at the 1% level  and is certainly sufficient for the analysis of $`G_M^s(0)`$.

An examination of the symmetries manifest in the QCD path integral for current matrix elements reveals various relationships among the quark sector contributions . Current matrix elements of hadrons, such as the magnetic moment of the proton, are extracted from the three-point function, a time-ordered product of three operators. Generally, an operator exciting the hadron of interest from the QCD vacuum is followed by the current of interest, which in turn is followed by an operator annihilating the hadron back to the QCD vacuum. In calculating the three point function, one encounters two topologically distinct ways of performing the electromagnetic current insertion. Figure 1 displays skeleton diagrams for these two possible insertions (with Euclidean time increasing to the right). These diagrams may be dressed with an arbitrary number of gluons and quark loops. The left-hand diagram illustrates the connected insertion of the current to one of the “valence”111Note that the term “valence” used here differs from that commonly used in the framework of deep-inelastic scattering. Here “valence” simply describes the quark whose quark flow line runs continuously from $`0x_2`$. These lines can flow backwards as well as forwards in time and therefore have a sea contribution associated with them . quarks of the baryon. In the right-hand diagram the external field produces a $`q\overline{q}`$ pair which in turn interacts with the valence quarks of the baryon via gluons. It is important to realize that within the lattice QCD calculation of the loop diagram on the right in Fig. 1 there is no antisymmetrization (Pauli blocking) of the quark in the loop with the valence quarks. For this reason, in general only the sum of the two processes in Fig. 1 is physical.

## 2 CHARGE SYMMETRY EQUALITIES

Under the assumption of charge symmetry, the three-point correlation functions for octet baryons leads to the following equalities for electromagnetic current matrix elements :

$`p=`$ $`e_uu_p+e_dd_p+O_N,n`$ $`=e_du_p+e_ud_p+O_N,`$
$`\mathrm{\Sigma }^+=`$ $`e_uu_{\mathrm{\Sigma }^+}+e_ss_\mathrm{\Sigma }+O_\mathrm{\Sigma },\mathrm{\Sigma }^{}`$ $`=e_du_{\mathrm{\Sigma }^+}+e_ss_\mathrm{\Sigma }+O_\mathrm{\Sigma },`$
$`\mathrm{\Xi }^0=`$ $`e_ss_\mathrm{\Xi }+e_uu_{\mathrm{\Xi }^0}+O_\mathrm{\Xi },\mathrm{\Xi }^{}`$ $`=e_ss_\mathrm{\Xi }+e_du_{\mathrm{\Xi }^0}+O_\mathrm{\Xi }.`$ (2)
Here, $`O`$ denotes the contributions from the quark-loop sector – shown on the right-hand side of Fig. 1. The baryon label represents the magnetic moment. Subscripts allow for environment sensitivity implicit in the three-point function . For example, the three-point function for $`\mathrm{\Sigma }^+`$ is the same as for the proton, but with $`d`$ replaced by $`s`$. Hence, the $`u`$-quark propagators in the $`\mathrm{\Sigma }^+`$ are multiplied by an $`s`$-quark propagator, whereas in the proton the $`u`$-quark propagators are multiplied by a $`d`$-quark propagator. The different mass of the neighboring quark gives rise to an environment sensitivity in the $`u`$-quark contributions to observables, which means that the naive expectations of the constituent quark model $`u_p/u_{\mathrm{\Sigma }^+}=u_n/u_{\mathrm{\Xi }^0}=1`$ may not be satisfied . This observation should be contrasted with the common assumption that the quark magnetic moment is an intrinsic-quark property which is independent of the quark’s environment.

Focusing now on the nucleon, we note that for magnetic properties, $`O_N`$ contains sea-quark-loop contributions from primarily $`u`$, $`d`$, and $`s`$ quarks. In the SU(3)-flavor limit ($`m_u=m_d=m_s`$) the charges add to zero and hence the sum vanishes. However, the heavier strange quark mass allows for a result which is non-zero. By definition

$`O_N`$ $`=`$ $`{\displaystyle \frac{2}{3}}{}_{}{}^{\mathrm{}}G_{M}^{u}{\displaystyle \frac{1}{3}}{}_{}{}^{\mathrm{}}G_{M}^{d}{\displaystyle \frac{1}{3}}{}_{}{}^{\mathrm{}}G_{M}^{s},`$ (3)
$`=`$ $`{\displaystyle \frac{{}_{}{}^{\mathrm{}}G_{M}^{s}}{3}}\left({\displaystyle \frac{1{}_{}{}^{\mathrm{}}R_{d}^{s}}{{}_{}{}^{\mathrm{}}R_{d}^{s}}}\right),\text{where}{}_{}{}^{\mathrm{}}R_{d}^{s}{\displaystyle \frac{{}_{}{}^{\mathrm{}}G_{M}^{s}}{{}_{}{}^{\mathrm{}}G_{M}^{d}}},`$ (4)
and the leading superscript, $`\mathrm{}`$, reminds the reader that the contributions are loop contributions. Note that, in deriving Eq.(4), we have set $`{}_{}{}^{\mathrm{}}G_{M}^{u}={}_{}{}^{\mathrm{}}G_{M}^{d}`$, corresponding to $`m_u=m_d`$ . In the constituent quark model $`{}_{}{}^{\mathrm{}}R_{d}^{s}=m_d/m_s0.65`$.

With no more than a little accounting, the strange-quark loop contributions to the nucleon magnetic moment, $`G_M^s`$ may be isolated from (2) and (4) in the following two phenomenologically useful forms,

$$G_M^s=\left(\frac{{}_{}{}^{\mathrm{}}R_{d}^{s}}{1{}_{}{}^{\mathrm{}}R_{d}^{s}}\right)\left[2p+n\frac{u_p}{u_{\mathrm{\Sigma }^+}}\left(\mathrm{\Sigma }^+\mathrm{\Sigma }^{}\right)\right],$$
(5)
and

$$G_M^s=\left(\frac{{}_{}{}^{\mathrm{}}R_{d}^{s}}{1{}_{}{}^{\mathrm{}}R_{d}^{s}}\right)\left[p+2n\frac{u_n}{u_{\mathrm{\Xi }^0}}\left(\mathrm{\Xi }^0\mathrm{\Xi }^{}\right)\right].$$
(6)
As we explained above, under the assumption that quark magnetic moments are not environment dependent, these ratios (i.e. $`u_p/u_{\mathrm{\Sigma }^+}`$ and $`u_n/u_{\mathrm{\Xi }^0}`$) are taken to be one in many quark models. Incorporating the experimentally measured baryon moments leads to:

$$G_M^s=\left(\frac{{}_{}{}^{\mathrm{}}R_{d}^{s}}{1{}_{}{}^{\mathrm{}}R_{d}^{s}}\right)\left[3.673\frac{u_p}{u_{\mathrm{\Sigma }^+}}\left(3.618\right)\right],$$
(7)
and

$$G_M^s=\left(\frac{{}_{}{}^{\mathrm{}}R_{d}^{s}}{1{}_{}{}^{\mathrm{}}R_{d}^{s}}\right)\left[1.033\frac{u_n}{u_{\mathrm{\Xi }^0}}\left(0.599\right)\right],$$
(8)
where all moments are expressed in nuclear magnetons $`(\mu _N)`$. (Note that the measured magnetic moments are all known sufficiently accurately that the experimental errors play no role in our subsequent analysis.) We stress that these expressions for $`G_M^s`$ are exact consequences of QCD, under the assumption of charge symmetry. Equation (8) provides a particularly favorable case for the determination of $`G_M^s`$ with minimal dependence on the valence-quark ratio.

If one considers the quark model suggestions of $`u_n/u_{\mathrm{\Xi }^0}=1`$ and $`{}_{}{}^{\mathrm{}}R_{d}^{s}=0.65`$ in (8), one finds $`G_M^s=0.81\mu _N`$, a significant departure from the experimental indication of positive values.

Equating (7) and (8) provides a linear relationship between $`u_p/u_{\mathrm{\Sigma }^+}`$ and $`u_n/u_{\mathrm{\Xi }^0}`$ which must be satisfied within QCD. Figure 2 displays this relationship by the dashed and solid line, the latter corresponding to values for which $`G_M^s(0)>0`$ when $`{}_{}{}^{\mathrm{}}R_{d}^{s}`$ is in the anticipated range $`0<{}_{}{}^{\mathrm{}}R_{d}^{s}<1`$. Since the line does not pass through the point $`(1,1)`$ corresponding to the simple quark model assumption of universality, the experimentally measured baryon moments are signaling that there must be an environment effect exceeding 12% in both ratios or approaching 20% or more in at least one of the ratios. Moreover, a positive value for $`G_M^s(0)`$ requires an environment sensitivity exceeding 70% in the $`u_n/u_{\mathrm{\Xi }^0}`$ ratio. Hence the experimental indication that $`G_M^s(0)>0`$ challenges the intrinsic magnetic moment concept which is fundamental to the constituent quark model.

## 3 CHIRAL CORRECTIONS

One of the major challenges at present in connecting lattice calculations of hadronic properties with the physical world is that computational limitations restrict the accessible quark masses to values an order of magnitude larger than the physical values. At such large masses one is far from the region where chiral perturbation theory is applicable. Yet one knows that for current quark masses near zero there is important non-analytic structure (as a function of the quark mass) which must be treated correctly if we are to compare with physical hadron properties. Our present analysis of the strangeness magnetic form factor has been made possible by a recent breakthrough in the treatment of these chiral corrections for the nucleon magnetic moments . In particular, a study of the dependence of the nucleon magnetic moments on the input current quark mass, within a chiral quark model which was fitted to existing lattice data, suggested a model independent method for extrapolating baryon magnetic moments which satisfied the chiral constraints imposed by QCD. We briefly summarize the main results of that analysis:

* a series expansion of $`\mu _{p(n)}`$ in powers of $`m_\pi `$ is not a valid approximation for $`m_\pi `$ larger than the physical mass,
* on the other hand, the behavior of the model, after adjustments to fit the lattice data at large $`m_\pi `$ is well determined by the simple Padé approximant:

$$\mu _{p(n)}=\frac{\mu _0}{1\frac{\chi }{\mu _0}m_\pi +cm_\pi ^2},$$
(9)
as illustrated in Fig. 3.
* Eq.(9) not only builds in the usual magnetic moment of a Dirac particle at moderately large $`m_\pi ^2`$ but has the correct leading non-analytic (LNA) behavior of chiral perturbation theory

$$\mu =\mu _0+\chi m_\pi ,$$
(10)
with $`\chi `$ a model independent constant.
* fixing $`\chi `$ at the value given by chiral perturbation theory and adjusting $`\mu _0`$ and $`c`$ to fit the lattice data yielded values of $`\mu _p`$ and $`\mu _n`$ of $`2.85\pm 0.22\mu _N`$ and $`1.96\pm 0.16\mu _N`$, respectively, at the physical pion mass. These are in remarkably good agreement with the experimental values – certainly much closer than the usual linear extrapolations in $`m_q`$.

Clearly it is vital to extend the lattice calculations of the baryon magnetic moments to lower values of $`m_\pi `$ than the 600 MeV used in the study just outlined. It is also vital to include dynamical quarks. Nevertheless, the apparent success of the extrapolation procedure gives us strong encouragement to investigate the same approach for resolving the strange quark contribution to the proton magnetic moment.

Our focus here is to extrapolate existing lattice QCD estimates of valance quark contributions to baryon moments $`(u_p,u_n,u_{\mathrm{\Sigma }^+},u_{\mathrm{\Xi }^0})`$ to the physical mass regime. The leading non-analytic pieces of these various contributions to the baryon magnetic moments are proportional to

$$\beta \frac{m_N}{8\pi f_\pi ^2}m_\pi \chi m_\pi ,$$
(11)
where the pion decay constant $`f_\pi =93`$ MeV. Isolating a particular quark flavour contribution only requires setting the electric charge of all other quark flavours to zero. Separation of the contributions of the LNA terms into valance and sea-quark sectors has been resolved in Ref. . The coefficients $`\beta `$ and $`\chi `$ for various quark sectors and baryons are summarized in Table 1. In removing the $`u`$ sea-quark contribution from the $`u`$-quark sector, one eliminates the contributions of the $`\overline{u}`$ sea-quark to the LNA term. Since the $`z`$-component of angular momentum in the pion-baryon system is positive, the $`\overline{u}`$ contribution to the moment is negative. Hence its removal suppresses the non-analyticity in the neutron moment and enhances the non-analyticity in $`\mathrm{\Xi }^0,p,`$ and $`\mathrm{\Sigma }^+`$ where the pion cloud acts to enhance the $`u`$-quark sector contribution to these baryon moments.

For the calculation of $`G_M^s(0)`$, the LNA chiral behavior of $`u_n`$ and $`u_{\mathrm{\Xi }^0}`$ is crucial, as the corresponding coefficients are of opposite sign. As a result, the ratio of these $`u`$-quark contributions, extrapolated with the correct chiral behavior, will be completely different from the linearly extrapolated ratio. That the signs of these chiral coefficients are opposite may be traced back to the charge of the predominant pion cloud associated with these baryons. The chiral behavior in the neutron case is dominated by the transition $`np\pi ^{}n`$, where there is a contribution from a $`\overline{u}`$ loop. On the other hand, for the $`\mathrm{\Xi }^0`$ the LNA $`u`$-quark contribution is dominated by the process $`\mathrm{\Xi }^0\mathrm{\Xi }^{}\pi ^+\mathrm{\Xi }^0`$. In this case the virtual $`\pi ^+`$ involves a valence $`u`$-quark – and hence the sign change.

In order to extrapolate the lattice data for valence quarks to the physical pion mass we use the same approach which worked so well for the neutron and proton moments in Ref. . That is, we fit the lattice data for the valence quark of flavor $`q`$ in baryon $`B`$ with the form:

$$\mu (m_\pi )=\frac{\mu ^0}{1\frac{\chi }{\mu ^0}m_\pi +cm_\pi ^2}.$$
(12)
Here $`\chi `$ are the modified LNA chiral coefficients for the valence quarks given in Table 1 (last column) and $`\mu ^0`$ and $`c`$ are fitting parameters.

Figures 4 and 5 illustrate the extrapolations of the lattice data. The results of these fits shown at $`m_\pi =140`$ MeV are the results of the extrapolation procedure, including fitting errors. Naive linear extrapolations are also shown to help emphasize the dominant role of chiral symmetry in the extrapolation process.

Clearly the effect of the chiral corrections on the extrapolated values of the key magnetic moment ratio $`u_n/u_{\mathrm{\Xi }^0}`$ at $`m_\pi =140`$ MeV is dramatic. Instead of the value $`0.72\pm 0.46`$ obtained by linear extrapolation one now finds $`1.51\pm 0.37`$. It is vital that this ratio is much more consistent with the range of values found necessary to yield a positive value of $`G_M^s`$, as illustrated in Fig. 2 – values which seemed quite unreasonable in the constituent quark model. By comparison, the ratio $`u_p/u_{\mathrm{\Sigma }^+}`$ at $`m_\pi =140`$ MeV is quite stable, changing from $`1.14\pm 0.08`$ in the case of linear extrapolation to $`1.14\pm 0.11`$ when the correct chiral extrapolation is used. The lattice point illustrated in Fig. 2 is consistent with the constraint obtained from charge symmetry. This suggests that any charge symmetry breaking in the experimentally measured moments is small.

The effect of the chiral extrapolation on the ratio $`u_n/u_{\mathrm{\Xi }^0}`$ is the most dramatic, changing it from $`0.72\pm 0.46`$ to $`1.51\pm 0.37`$. Thus Eq. (8), which seemed certain to guarantee $`G_M^s<0`$, no longer does. Indeed, using the recent estimate  of $`{}_{}{}^{\mathrm{}}R_{d}^{s}=0.55`$, we obtain $`G_M^s=0.16\pm 0.18\mu _N`$. While this is still negative it does permit a small positive value within the error.


text:
# Comptonization in X-ray bright neutron star globular cluster systems

## REFERENCES

Bałucińska-Church M., et al., 1999, A&A, 349, 495

Barret D., et al., 1996, ApJ, 473, 963

Barret D., et al., 2000, ApJ, 533, 329

Boella G., et al., 1997, A&AS, 122, 299

Church M., et al. 1999, AA&, 338, 556

Guainazzi M., et al., 1998, A&A, 339, 802

King A., 1995, in ”X-ray binaries”, Lewin G., van Paradijs J., van den Heuvel E.P.J. eds., 439

in’ t Zand J., et al., 1999a, A&A, 345, 100

in ‘t Zand J., et al., 1999b, A&A, 347, 891

Parmar A.N., et al. 1999, A&A, 351, 225

Parmar A.N., et al., 2000, A&A, 356, 175


text:
# KPD 1930+2752 – a candidate Type Ia supernova progenitor.

## 1 Introduction

Type Ia supernovae (SNe Ia) are one of the most important tools for observational cosmology because there appears to be a relatively small spread in their peak optical brightness around $`M_\mathrm{V}=19.6`$ and they can be seen out to cosmological distances (z$``$1) so they can be used to measure cosmological parameters, e.g., the cosmological constant, $`\mathrm{\Lambda }`$ (Riess et al. 1998, Perlmutter et al. 1999). However, the peak optical brightnesses of SNe Ia are not uniform, they are correlated with the shape of the lightcurve and vary by about one magnitude. Meaningful measurements of cosmological parameters require this variation to be calibrated, e.g., the non-zero value of $`\mathrm{\Lambda }`$ measured by Perlmutter et al. is required by supernovae at z$``$0.5 being about 0.3 magnitudes too bright compared to a non-accelerating ($`\mathrm{\Lambda }`$=0) Universe. The corrections to peak brightnesses have to be empirical because it is still not yet clear what causes SNe Ia.

Type Ia supernovae near maximum light show no hydrogen or helium lines but do show strong silicon lines. The absence in the spectrum of the two most common elements in the Universe dramatically reduces the number of potential progenitors, as does their appearance in old stellar populations, e.g., elliptical galaxies. All the most likely models for progenitors feature an accreting white dwarf (Leibundgut 2000; Branch et al. 1995) which ignites carbon in its core either because it has reached the Chandrasekhar mass (1.4$`M_{}`$) or because ignition of accumulated helium causes compression of the core and a so-called “edge-lit detonation”. This explains the fast rise times for SNe Ia, the lack of hydrogen and helium and the fairly uniform peak brightness. To initiate the explosion, the white dwarf must accrete material from a companion star. Two models for the companion star which have gained popularity in recent times are super-soft sources and double degenerates.

Super-soft sources are white dwarfs which accrete hydrogen from a non-degenerate star at a rate just sufficient to support steady nuclear burning on the surface of the white dwarf (Kahabka & van den Heuvel 1997). This leads to an accumulation of material on the white dwarf but it is not clear whether or not the accretion rate stays within the required range sufficiently long for a SNe Ia explosion to result or that there are a sufficient number of these binaries to explain the observed rate of SNe Ia. Neither is it clear that these binaries have a sufficiently long lifetime to be observed in elliptical galaxies (Yungelson & Livio 2000; Yungelson et al. 1995).

The double degenerate model posits two white dwarfs with an orbital period of a few hours which merge due to the loss of gravitational wave radiation. Two drawbacks with this model have been a poor understanding of how such a detonation might be initiated and the lack of observed progenitors. There are now many double degenerates known (Marsh, Dhillon & Duck 1995; Moran, Marsh & Bragaglia, 1997; Moran, Maxted & Marsh 2000) but none have both a sufficiently short orbital period and a total mass in excess of 1.4$`M_{}`$. However, at least one good candidate for the progenitor of an edge-lit detonation is known (WD 1704+481.2; Maxted et al. 2000).

KPD 1930+2752 was identified as a subdwarf-B (sdB) star in the Kitt Peak – Downes survey of UV excess objects near the Galactic plane (Downes 1986). Photometry in the Cousins BVRI system by Allard et al. (1994) and Strömgren photometry by Wesemael et al. (1992) revealed nothing exceptional about this star other than that it has a low reddening and that there is no evidence for companion to this star. Of the 100 subdwarfs in the study of Allard et al., 31 show evidence for a companion which, when some estimate of the selection effects is made, suggests that more than half of these stars are binaries.

The binary fraction of sdB stars is a matter of some interest because the properties of these stars suggest they have lost a substantial fraction of their mass, perhaps due to interactions with a companion star (Heber 1986; Saffer et al. 1994; Iben & Livio 1993). In addition to the composite spectrum binaries identified by Allard et al., several binary sdB stars have been identified from the Doppler shift of the spectral lines due to the orbital motion (Saffer, Livio & Yungelson 1998). This technique has the advantage of being sensitive to the companion star whatever its type and can detect sdB stars with white dwarf companions which would be missed by almost any photometric technique. Direct evidence for white dwarf companions to sdB stars is seen in the eclipsing sdB – white dwarf binary KPD 0422+5421 (Orosz & Wade 2000).

Other clues to the properties of sdB stars comes from the EC 14026 stars – sdB stars showing p-mode pulsations with periods of a few minutes (Koen et al. 1998; Fontaine et al. 1998; O’Donoghue et al. 1999). High speed photometry by Billères et al. (2000) identified KPD 1930+2752 as the 14th EC 14026 star known. The photometry also showed variability with a period of $`2^\mathrm{h}17\mathrm{m}`$ with an amplitude of 1.4 percent. This is much longer than the pulsation periods of EC 14026 stars. The lightcurve folded on this period shows an almost sinusoidal variation with two minima per cycle, but with one minimum being slightly deeper than the other. This was interpreted as being the ellipsoidal variation due to the rotation of a star distorted by the presence of a companion star.

In this paper we present spectroscopy of the the H$`\alpha `$ and HeI 6678Å spectral lines which confirms that KPD 1930+2752 is a binary star with an orbital period of $`2^\mathrm{h}17^\mathrm{m}`$. We also show that the total mass of the binary exceeds the Chandrasekhar mass and conclude that KPD 1930+2752 is the first good candidate Type Ia supernova progenitor which will explode due to the accretion of helium and other elements heavier than hydrogen onto a white dwarf on an astrophysically interesting timescale.

## 2 Observations and reductions

Observations were obtain with the Isaac Newton Telescope at the Observatario Roque de los Muchachos on the Island of La Palma. We used the Intermediate Dispersion Spectrograph with the 500mm camera, a 1200 lines/mm grating and a TEK charged-coupled device detector to obtain 25 spectra of KPD 1930+2752 covering 400Å around the H$`\alpha `$ line with a dispersion of 0.39Å per pixel. The resolution measured from the full-width at maximum of the arc lines is 0.9Å. All the spectra were obtained in a single run of observations of just over 2 hours on the morning of 17 April 2000. The slit width used was 0.97 arcsec, which was well matched to seeing estimated from the spatial profile of the spectra of around 1.3 arcsec. Observations of a CuNe arc were obtained before and after the run of observations and every 25 minutes in-between. The exposure time used for all the spectra was 300s.

Extraction of the spectra from the images was performed automatically using optimal extraction to maximize the signal-to-noise of the resulting spectra (Horne 1986). The arcs associated with each stellar spectrum were extracted using the same weighting determined for the stellar image to avoid possible systematic errors due to tilted spectra. The wavelength scale was determined from a fourth-order polynomial fit to measured arc-line positions. The standard deviation of the fit to the 8 spectra lines was typically 0.09Å. The wavelength scale for an individual spectrum was determined by interpolation to the time of mid-exposure from the fits to arcs taken before and after the spectrum to account for the small amount of drift in the wavelength scale ($`<0.1`$Å) due to flexure of the instrument. Statistical errors on every data point calculated from photon statistics are rigorously propagated through every stage of the data reduction.

## 3 The spectroscopic orbit

Visual inspection of the H$`\alpha `$ and HeI 6678Å spectral lines shows clearly the Doppler motion expected from a binary star with an orbital period of $`2^\mathrm{h}17^\mathrm{m}`$ with a semi-amplitude of $``$350 km s-1, confirming the interpretation of Billères et al. To obtain a more accurate value of for the semi-amplitude of the orbit, we created a model of the spectrum composed of three Gaussian profiles for the H$`\alpha `$ line and a single Gaussian profile for the HeI 6678Å line. We used a least-squares fit to the first spectrum to obtain an initial estimate of the shape of the model spectrum. We then fixed the shape of the model spectrum and varied only the position of the lines in a least-squares fit to each of the spectra to obtain an initial estimate of the radial velocity from each spectrum. These radial velocities are shown in Fig. 1. The spectra were normalised prior to fitting using a quadratic fit to the continuum either side of the H$`\alpha `$ and HeI 6678Å spectral lines. Only data in the range 6500–6675Å were included in the fit.

We then used a simultaneous least-squares fit to all the spectra to determine the best profile shape and to determine the spectroscopic orbit simultaneously. The position of the model profile varies from spectrum-to-spectrum according to the radial velocity predicted by $`\gamma +K\mathrm{sin}((TT_0)/P)`$, where $`\gamma `$ is the systemic velocity, $`K`$ is the projected orbital speed, $`T`$ is the time of mid-exposure of the spectrum and $`P`$ is the orbital period. We fixed the value of $`P`$ at the value given by Billères et al. (8217.8s = 0.095111d) since this value is much more accurately determined by their photometry than can be done from our spectroscopy. The time $`T_0`$ corresponds to the point in the orbit when the sdB star is closest to the observer. The smearing of the spectra due to the motion during the exposure is incorporated into the model. We used independent values of $`\gamma `$ for the H$`\alpha `$ and HeI 6678Å spectral lines to allow for any difference due to pressure shifts. There are many free parameters in this fit so we first varied the parameters in the spectroscopic orbit while holding the model profile fixed and then vice versa so all the parameters were near their optimum values. We then varied all the parameters to find the solution given in Table 1. The fit to the data is very good judged both by the value of $`\chi ^2`$ and visual inspection of the residuals.

## 4 Physical properties of the binary

The projected orbital speed, $`K=349.3\pm 2.7`$ km s-1, and orbital period, $`2^\mathrm{h}17^\mathrm{m}`$, immediately imply a minimum mass for the unseen companion star of 0.42$`\pm `$0.01$`M_{}`$. More realistically, the measured effective temperature and surface gravity of KPD 1930+2752 (Teff=33 000K, $`\mathrm{log}g=5.61`$) is typical for sdB stars and places it squarely in the region of the Teff – $`\mathrm{log}g`$ plane occupied by models of core helium-burning stars with masses of 0.5$`M_{}`$ and very thin hydrogen envelopes ($`<0.02M_{}`$, Saffer et al. 1994), in which case the mass of the companion star is at least 0.97$`\pm `$0.01$`M_{}`$. The lightcurve of KPD 1930+2752 observed by Billères et al. after removal of the signal due to pulsations is shown in Fig. 2. The quasi-sinusoidal signal with unequal minima characteristic of a star distorted by a close companion is apparent. We can produce a model lightcurve using the physical parameters derived above and assuming that the orbital inclination is 90. The radius of the sdB star implied by the surface gravity and canonical mass is 0.18$`\pm `$0.01$`R_{}`$ and the separation of the stars is 0.98$`R_{}`$. The orbit is far too small to contain a normal star of 0.97$`\pm `$0.01$`M_{}`$ or more, so we assume that the unseen companion is a white dwarf star. Other parameters of the model are the gravity darkening exponent of the sdB star, for which we assume the standard value appropriate for radiative stars, and the limb-darkening. The lightcurve was obtained with a blue-sensitive detector so we use a linear limb-darkening coefficient of 0.29 which is the mean of the values for U and B filters for a $`\mathrm{log}g=5`$, Teff=33000 model atmosphere given by Diaz-Cordoves, Claret & Gimenez (1995). The precise value of the limb-darkening or gravity darkening exponent used has very little effect on the lightcurve. An additional effect included in the model is the Doppler boosting due to the high orbital velocity. This effect increases the total flux by a factor $`(1v(t)/\mathrm{c})^3`$ where $`v`$ is the radial velocity at time $`t`$ and c is the speed of light. This effect is counteracted by the Doppler shift, which reduces the effect by a factor $`(1v(t)/\mathrm{c})^2`$ for observations on the Rayleigh-Jeans tail of a black-body spectrum, which is a good approximation to the spectrum of KPD 1930+2752 in the optical region. The overall effect is to make the maxima of the lightcurve asymmetric. This is shown in Fig. 2 by plotting a model lightcurves both with and without Doppler boosting. The agreement between either model and the observed lightcurve is excellent. Note that that there has been no attempt made to optimize the parameters of the model, the lightcurve is prediction based purely on the measured surface gravity, effective temperature and orbital velocity together with the assumption that the mass of the sdB star is 0.5$`M_{}`$ and that the inclination is 90. One effect excluded from our model is the transit of the companion star across the face of the sdB star. If the companion is a white dwarf star, its radius will be approximately 0.01$`R_{}`$ so the eclipse depth will be about $`(0.01/0.18)^2=0.3`$ percent. There is a hint of just such a feature at the correct phase (0.5) in the lightcurve, but improved photometry will be required to confirm this feature.

## 5 Discussion

KPD 1930+2752 is very similar to the star KPD0422+5421, which is an sdB–white dwarf binary with an orbital period of 2.16h and an projected orbital speed of 237$`\pm `$18 km s-1 (Orosz & Wade 2000). The ellipsoidal variation in KPD0422+5421 has almost exactly the same amplitude as that seen KPD 1930+2752 and the transit of the white dwarf across the face of the sdB star is seen in high quality lightcurves.

Of course, the larger orbital velocity of KPD 1930+2752 implies a larger total mass for this binary than KPD0422+5421. The total mass of KPD 1930+2752 assuming that the sdB star has a mass of 0.5$`M_{}`$ is at least $`1.47\pm 0.01`$$`M_{}`$. The significance of this result is that the total mass exceeds the Chandrasekhar limit for white dwarfs (1.40$`M_{}`$, Hamada & Salpeter 1961). KPD 1930+2752 will merge within about 200 million years due to a combination of orbital shrinkage through gravitational wave radiation and the evolutionary expansion of the sdB star. Thus, KPD 1930+2752 is the first star to be discovered which is a good candidate for the progenitor of a Type Ia supernova due to accretion of helium and other elements heavier than hydrogen onto the white dwarf which then exceeds the Chandrasekhar mass.

## 6 Conclusion

We have confirmed the conclusion of Billères et al. that KPD 1930+2752 is a binary star in which the sdB star shows ellipsoidal variability and the orbital period is $`2^\mathrm{h}17^\mathrm{m}`$. The amplitude of the orbital motion (349.3$`\pm `$2.7 km s-1) combined with the canonical mass for sdB stars (0.5$`M_{}`$) implies a total mass for the binary of 1.47$`\pm `$0.01$`M_{}`$. The unseen companion star is almost certainly a white dwarf star. The binary will merge within $``$200 million years due gravitational wave radiation. The accretion of helium and other elements heavier than hydrogen onto the white dwarf which then exceeds the Chandrasekhar mass (1.4$`M_{}`$) is a viable model for the cause of Type Ia supernovae. KPD 1930+2752 is the first star to be discovered which is a good candidate for the progenitor of a Type Ia supernova of this type which may explode on an astrophysically interesting timescale.

## Acknowledgements

PFLM was supported by a PPARC post-doctoral grant. The Isaac Newton Telescope is operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias.