# Polarization and Strong Infra-Red Activity in Compressed Solid Hydrogen

## Abstract

Under a pressure of $`150`$ GPa solid molecular hydrogen undergoes a phase transition accompanied by a dramatic rise in infra-red absorption in the vibron frequency range. We use the Berry’s phase approach to calculate the electric polarization in several candidate structures finding large, anisotropic dynamic charges and strongly IR-active vibron modes. The polarization is shown to be greatly affected by the overlap between the molecules in the crystal, so that the commonly used Clausius-Mossotti description in terms of polarizable, non-overlapping molecular charge densities is inadequate already at low pressures and even more so for the compressed solid.

The stretching mode of an isolated hydrogen molecule does not absorb infra-red (IR) radiation, since the electric dipole moment vanishes by symmetry and remains zero upon stretching. In the molecular solid the symmetry is lower, and for some structures first-order IR activity of vibrational modes is allowed. Indeed, some IR absorption occurs in the broken-symmetry phase (BSP or phase II) . More surprising is the dramatic increase in the intensity of IR absorption which occurs when crossing the boundary between the BSP phase and phase III (H-A) at around 150 GPa. At 85 K the vibron oscillator strength goes from $`f0.0002`$ at 140 GPa to $`f0.0134`$ at 167 GPa.

A considerable effort has been made to try to explain this remarkable behaviour. Hemley and co-workers proposed that the increase in molecular overlap leads in phase III to the formation of charge-transfer states between neighboring molecules, with the coupling between a vibron and charge-transfer excitations giving rise to charge oscillations between the molecules. Baranowski invoked a strong ionization of the molecules in phase III. Mazin el al neglected the molecular overlap and treated the molecules as point objects which become polarized in the quadrupolar field of the other molecules (EQQ model ). Edwards and Ashcroft suggested that a dielectric instability associated with a charge-density wave causes an enhancement of the molecular dipole moments in phase III. Some theories focus on purely static charge, others on purely dynamic (i.e., displacement-induced) and the others  take the two contributions to be comparable. There is also no agreement on the relative importance of intramolecular  and intermolecular charge transfer . To settle these issues, a careful treatment of the bulk electric polarization, which is the central quantity at play in the physics of IR absorption, is needed.

For that purpose we have used the “modern theory of polarization” , which shows how to compute the polarization of a periodic insulating system as a Berry’s phase derived from the electronic wave functions. The information about the polarization is therefore in the phases of the wave functions, not in the periodic density alone, whose dipole moment is in general ill-defined, depending on the choice of unit cell. The density-based picture is valid only in the Clausius-Mossotti limit of localized, polarizable units with non-overlapping charge densities and breaks down whenever significant charge delocalization prevents an unambiguous assignment of the electron density to particular atoms, molecules or ions in the solid. Such overlap effects are known to play an important role in the absorption of light by lattice vibrations in ionic systems , which typically have a much larger overlap than molecular crystals. But at 150 GPa solid hydrogen has been compressed by almost a factor of ten, hence the motivation for applying the modern theory of polarization to this problem.

Using the Berry’s phase approach we have calculated the Born effective charge tensors $`Z^{}`$. In order to compare with experiments done on polycrystalline samples, we used the expression for the oscillator strength $`f(j)`$ of the contribution of a TO phonon mode j to the dielectric function averaged over all directions,

$$f(j)=\frac{4\pi }{\mathrm{\Omega }\omega (j)^2}\left(\frac{1}{3}\underset{\alpha =1}{\overset{3}{}}\xi _\alpha (j)^2\right),$$

where $`\xi _\alpha (j)\frac{1}{\sqrt{m_\mathrm{p}}}_\mathrm{k}_\beta Z_{\alpha \beta }^{}\left(k\right)e_\beta (k;j)`$, $`k`$ indexes the atoms in the unit cell of volume $`\mathrm{\Omega }`$, $`e_\beta (k;j)`$ are the normalized eigenvectors associated with the mode, $`\omega (j)`$ is its angular frequency and $`m_\mathrm{p}`$ is the proton mass.

The crystal structure in the high-pressure phases has not yet been determined experimentally. For clarity we have investigated two of the simplest candidate structures for the compressed solid which have IR-active vibron modes and are also energetically favorable. The parameters $`a`$ (lattice constant), $`c/a`$, bond length and $`\theta `$ of $`Cmc2_1`$ and $`C2/m`$ (see Fig. 1) were optimized by minimizing the enthalpy at fixed pressure using a variable cell shape method, and they are collected in Table I for several densities. In some cases the centers of the molecules were allowed to move away from the hcp sites by symmetrical amounts in the yz plane, which does not lower the symmetry while allowing for the manifestation of the instability reported in Ref.. The calculations were done in the local density approximation (LDA) to the density functional theory, with clamped nuclei.

FIG. 1. The $`Cmc2_1`$ structure viewed along the c-axis (left) and in the $`yz`$ plane (right). The centers of the molecules lie on hcp sites, and the molecules in the two sublattices are tilted away from the c-axis by opposite angles $`\theta `$ and $`\theta `$. The $`C2/m`$ structure is identical except that the two molecules are tilted in the same direction by an angle $`\theta `$.

| Structure | P | $`r_\mathrm{s}`$ | c/a | $`r_{\mathrm{bond}}`$ | $`\theta `$ | $`f_{\mathrm{in}\mathrm{phase}}`$ | $`f_{\mathrm{out}\mathrm{phase}}`$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| $`C2/m`$ | 115 | 1.52 | 1.588 | 1.456 | 69.5 | 0 | 0.169 |
| $`C2/m^{}`$ | 115 | 1.52 | 1.583 | 1.460 | 70.6 | 0 | 0.589 |
| $`Cmc2_1`$ | 13 | 2.0 | 1.576 | 1.445 | 54.0 | $`10^6`$ | $`9\times 10^6`$ |
| $`Cmc2_1`$ | 115 | 1.52 | 1.576 | 1.445 | 54.0 | 0.00079 | 0.0114 |
| $`Cmc2_1^{}`$ | 115 | 1.52 | 1.574 | 1.451 | 56.6 | 0.014 | 0.049 |
| $`Cmc2_1`$ | 152 | 1.47 | 1.571 | 1.443 | 55.2 | 0.0026 | 0.026 |
| $`Cmc2_1`$ | 167 | 1.45 | 1.570 | 1.443 | 55.6 | 0.004 | 0.034 |
| $`Cmc2_1`$ | 180 | 1.43 | 1.569 | 1.443 | 55.9 | 0.0059 | 0.044 |
| Phase III | 167 |  |  |  |  | undetected | 0.0134 |

Centers of the molecules are off-site (see text). The off-site displacement in a.u. of the upper molecule in Fig. 1 is $`\delta y=0.1`$, $`\delta z=0.01`$ ($`C2/m`$), and $`\delta y=0.07`$ ($`Cmc2_1`$).

TABLE I. Optimized parameters (except for $`r_s=2.0`$) and vibron oscillator strengths for several pressures. Pressures are in GPa and are obtained using the equation of state of Ref., angles are in degrees, and the other quantities in atomic units (a.u.). The experimental results for phase III are taken from Ref.. “undetected” means that the Raman-active vibron is not observed in the IR spectra.

The existence of a center of inversion between the two molecules in the unit cell of $`C2/m`$ leads to the relation $`Z^{}(1)=Z^{}(2)=Z^{}(3)=Z^{}(4)`$ between the effective charge tensors of the four atoms. For the optimized structure at 115 GPa we obtained, for atom 2 before and after the relaxation away from the hcp sites, respectively:   

$`Z^{}\left(\begin{array}{ccc}<0.005& 0& 0\\ 0& 0.53& 0.18\\ 0& 0.08& 0.07\end{array}\right)`$, $`Z^{}\left(\begin{array}{ccc}0.05& 0& 0\\ 0& 0.96& 0.28\\ 0& 0.15& 0.16\end{array}\right)`$.   

Although $`Z_{xx}^{}`$ is non-zero by symmetry, we found it to be very small before the off-site relaxation. The pronounced anisotropy of $`Z^{}`$ implies large dynamic charges, since displacements of rigid ions (static charges) give rise to isotropic diagonal $`Z^{}`$ tensors. The smallness of $`Z_{xx}^{}`$ relative to the components in the yz submatrix suggests that most of the charge transfer occurs in the planes of the molecules, since a displacement $`\delta x`$ of an atom in the x direction changes the distances to other atoms in that plane to second order only. This analysis establishes that intramolecular and/or in-plane intermolecular charge transfer are the dominant effects, but it does not distinguish between the two . The small displacement away from the hcp sites has a large effect on $`Z^{}`$, revealing a strongly nonlinear dependence of the polarization on the positions of the molecules.

Such large effective charges give rise to very strong IR absorption by the vibron mode in which the two molecules in the primitive cell vibrate out-of-phase (the in-phase mode is Raman active but IR-inactive in this structure). Before the relaxation, the oscillator strength at 115 GPa is about 13 times larger than the experimentally measured value in phase III at 167 GPa, and it increases by a factor of 3.5 upon relaxation (see Table I). Due to closure of the LDA band gap, this structure was not studied at significantly higher pressures.

In $`Cmc2_1`$ the symmetry of the structure dictates the following form for the effective charge tensors:   

$`Z^{}(1)=\left(\begin{array}{ccc}a& 0& 0\\ 0& b& c\\ 0& d& e\end{array}\right)`$, $`Z^{}(2)=\left(\begin{array}{ccc}a& 0& 0\\ 0& b& c^{}\\ 0& d^{}& e\end{array}\right)`$,   
$`Z^{}(3)=\left(\begin{array}{ccc}a& 0& 0\\ 0& b& c\\ 0& d& e\end{array}\right)`$, $`Z^{}(4)=\left(\begin{array}{ccc}a& 0& 0\\ 0& b& c^{}\\ 0& d^{}& e\end{array}\right)`$.   

| $`r_\mathrm{s}`$ | a | b | c | d | e | c’ | d’ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 2.0 | -0.0178 | -0.0033 | 0.028 | 0.044 | 0.0005 | 0.042 | 0.040 |
| 1.52 | -0.022 | -0.14 | 0.069 | 0.193 | -0.0088 | 0.208 | 0.271 |
| $`1.52^{}`$ | -0.050 | -0.30 | -0.014 | 0.091 | -0.0087 | 0.243 | 0.348 |
| 1.47 | -0.020 | -0.20 | 0.069 | 0.216 | -0.015 | 0.260 | 0.347 |
| 1.45 | -0.018 | -0.229 | 0.068 | 0.223 | -0.017 | 0.283 | 0.380 |
| 1.43 | -0.017 | -0.254 | 0.065 | 0.228 | -0.020 | 0.307 | 0.413 |

Centers of the molecules are off-site (see text).

TABLE II. Effective charges in atomic units for the optimized $`Cmc2_1`$ at several densities.   

The calculated values for several densities are shown in Table II. The anisotropy is again very large, which implies large dynamic charges, but the overall values are somewhat smaller than for $`C2/m`$. Both vibrons are IR-active in this structure: the out-of-phase vibron has a large oscillator strength in reasonable agreement with the measured values in phase III (see Table I), whereas the in-phase vibron is about an order of magnitude weaker, which would still be enough to be detected experimentally, if $`Cmc2_1`$ was the correct structure for that phase. The intensity of the absorption grows rapidly with pressure, in agreement with the behavior in phase III. At 167 GPa we obtain for the Szigeti effective charge $`q^{}`$ of the out-of-phase vibron $`q^{}0.047e`$, comparable to the measured value of $`q^{}0.032e`$, and its pressure derivative is $`dq^{}/dP2.3\times 10^4`$ e/GPa, close to the experimental value of $`2.1\times 10^4`$ e/GPa. The effect of off-site relaxation of the molecules is again significant, with the oscillator strengths of the vibrons increasing 18-fold (in-phase) and 4-fold (out-of-phase) at 115 GPa.

FIG. 2. (a) electron density $`n(𝐫)`$ in a.u. for $`Cmc2_1`$ along a line joining atoms 1(3) and 2(4) of Fig. 1. (b) asymmetry along that line, defined as $`a(\xi )=\left[n(\xi )n(\xi )\right]/2n_0`$, where $`n_0`$ is the maximum electronic density in the hydrogen atom. For $`r_s=3.4`$ $`a(\xi )`$ was multiplied by $`\left(3.4/2.0\right)^4`$ for comparison with $`r_s=2.0`$. The molecular centers are on hcp sites, and therefore the intramolecular $`a(\xi )`$ is very small, and yet the vibron oscillator strength for $`r_s=1.52`$ is quite large.

More insight into the physics of polarization can be obtained by looking at the electron density. In Fig. $`2a`$ we show the density along the axis of a molecule in the $`Cmc2_1`$ structure for several pressures. The asymmetry induced by the neighboring molecules can be clearly seen in the higher pressure curve, and is more visible in Fig. $`2b`$, where we have plotted the asymmetry function $`a(\xi )`$. If the molecules can be considered to be separate polarizable objects, as in the Clausius-Mossotti model, then the asymmetry of the molecular density will determine the polarization. For example, in the widely used EQQ model  the molecular dipoles are induced by the electric field at the centers of the molecules due to the quadrupole moments $`Q`$ of the other molecules. In this model the intramolecular asymmetry scales as $`1/r_s^4`$ (assuming that $`Q`$ and the molecular polarizability do not change much with pressure). Indeed this is the case at low and moderate pressures ( $`r_s2.0`$), as can be seen in Fig. $`2b`$ from the agreement in the intramolecular region between the curve for $`r_s=2.0`$ and the curve for $`r_s=3.4`$ scaled by $`\left(3.4/2.0\right)^4`$. However, the scaling clearly does not hold up to the megabar range, since the intramolecular $`a(\xi )`$ would have to be around 3 times larger for $`r_s=1.52`$ than for $`r_s=2.0`$, whereas it is found to be slightly smaller and decreasing with pressure in that range. Therefore in the megabar range the EQQ interactions do not account for the local field acting on the molecules, contrary to the model used in Ref..

In Ref. $`a(\xi )`$ in the intramolecular region was interpreted as an estimate of the dipole moment of the molecule in the solid, which leads to the prediction that the vibron oscillator strength should scale roughly as the square of its magnitude. Instead, we have found that although the intramolecular asymmetry is comparable for $`r_s=1.52`$ and $`r_s=2.0`$, the oscillator strengths of the vibrons are $`800`$ (in-phase) and $`1300`$ (out-of-phase) times larger for $`r_s=1.52`$ (see Table I). The origin of this surprising result can already be seen in Fig. 2, where it is clear that for both these densities $`a(\xi )`$ is large in the intermolecular region, where the density is small, which is related to the fact that the outer regions of the molecules are most easily polarizable. Thus, when considering polarization effects, molecular overlaps may be important, even if they are small, causing the dipole moments of the molecules to be ill-defined and compromising the Clausius-Mossotti picture. 

At low enough densities the Claussius-Mossotti model must be correct; however, to our knowledge there has never been a thorough investigation of how low the density must be for this limit to apply. This is now possible using the modern theory of polarization. Here we present results for the $`Cmc2_1`$ structure as a function of volume. (We emphasize that this is not meant to represent real hydrogen at low pressure, where the molecules are not oriented but instead behave as quantum rotors, but is sufficient to establish the desired points regarding the sensitivity of the polarization to the molecular overlap.) In this structure there is a net spontaneous polarization $`\mathrm{P}_0`$ along the c-axis, and in Fig. 3 we plot $`\mathrm{P}_0`$ as a function of $`r_s`$, together with the value given by the EQQ model and the dipole moment for a particular choice of unit cell. The cell we used goes from $`z=0`$ to $`z=c`$ in Fig. 1, which is the obvious choice in the low pressure limit since the electronic density will vanish at the boundaries and the dipoles will become well-defined in that limit. For $`r_s1.52`$ atoms 1 and 3 in Fig. 1 have an excess of electron charge with respect to atoms 2 and 4, as inferred from $`a(\xi )`$. In a Clausius-Mossotti framework this would lead to a net polarization (sum of the molecular dipole moments) along minus z, and that is indeed the result of the EQQ model (and also of the dipole moment for $`r_s>2.5`$), but between $`r_s1.7`$ ($`50`$ GPa) and $`r_s3.0`$ (around atmospheric pressure) the Berry’s phase polarization has the opposite sign. Most surprising is the fact that the corrections to the models are so large even at low pressures (large $`r_s`$). Since the three quantities in Fig. 3 start to converge to each other only for $`r_s3.0`$, we conclude that the Clausius-Mossotti limit is reached only at negative pressures. A possible explanation for this observation is that at zero pressure there must still be some overlap of the wavefunctions in order to balance the attractive interactions, and only at negative pressure does this overlap become truly negligible.

FIG. 3. Spontaneous polarization per unit cell as a function of $`r_s`$ for the hcp-centered $`Cmc2_1`$. The structural parameters are the optimized parameters for $`r_s=1.52`$ (see Table I). The Berry’s phase result is compared with the value given by the EQQ model, using the quadrupole moment and polarizability of the H<sub>2</sub> molecule quoted in Ref., and with the “dipole moment” (see text) of the charge density in a unit cell.

In summary, we have presented the first ab-initio calculations of the intensity of IR absorption in solid hydrogen, showing that the intense vibron IR activity found in some structures is caused mainly by large dynamic charges. The strong dependence of the effective charges on the symmetry and structural parameters offers a plausible explanation for the difference in IR absorption between the two high-pressure phases, and the study of other structures with the present approach may help identify the crystal structure of those phases. The physics of polarization of solid hydrogen appears to be more subtle and less amenable to simple models than has been previously assumed. This seems to be caused by the very delocalized nature of the induced dipoles, which extend out to the exponentially decaying tails of the molecular charge density, rendering the bulk polarization very sensitive to even small overlaps between the molecules. These findings may be relevant for other molecular solids, since the Clausius-Mossotti picture is assumed in most treatments of the IR absorption in such systems.

It is a pleasure to thank J. L. Martins for providing the LDA program and for many illuminating discussions. We also benefited from conversations with N. W. Ashcroft, B. Edwards, A. F. Goncharov, R. Resta, S. Scandolo, I. F. Silvera and R. Zallen. The calculations were done in part at the NCSA. This work was supported by the DOE Grant No. DEFG02-96-ER45439 and the NSF grant No. DMR9422496. I.S. acknowledges financial support from JNICT (Portugal) and R.M.M. thanks the Aspen Center for Physics for its hospitality during part of this work.