text:
# Discovery of a 450 Hz QPO from the Microquasar GRO J1655-40 with RXTE

## 1 Introduction

The Galactic microquasar GRO J1655-40 is one of a handful of black hole systems in which high frequency quasiperiodic X-ray brightness oscillations have been recently observed with RXTE. Remillard et al. (1999) reported the presence of a weak, broad $`300`$ Hz QPO at times when the X-ray spectrum was hardest (dominated by a power law component) and the luminosity above $`0.2L_{Edd}`$. Recently, high frequency QPOs have also been observed in three black hole transients; 4U 1630-47 at 184 Hz, XTE J1858+226 at 170 Hz, and XTE J1550-564 between 100 - 283 Hz (see Remillard & Morgan 1998; Markwardt et al. 1999; Remillard et al. 1999a; Homan et al 2000). A stable 67 Hz QPO has also been observed from the Galactic microquasar GRS 1915+105 (Morgan, Remillard & Greiner 1997). The high frequencies of these QPO combined with their X-ray origin argue forcefully that they are produced in the innermost region of the accretion flow close to the black hole event horizon. If their production mechanism can be understood they stand to provide a wealth of information on black hole mass and spin as well as the structure of strongly curved spacetime. Their behavior may also provide tests of General Relativity (GR) in the strong field limit. A number of different models have been proposed to explain these oscillations, but all of them require that strong-field GR effects be taken into account (see Milsom & Taam 1997; Nowak et al. 1997; Wagoner 1998; Stella, Vietri & Morsink 1999; Merloni et al. 1999). Although pairs of simultaneously present high frequency QPOs have been observed ubiquitously from neutron star systems (see the review by van der Klis 2000), until now no black hole system has revealed a pair of simultaneous high frequency ($`\nu 50`$ Hz) QPO.

GRO J1655-40 is unique amongst these QPO sources in that it is the only one whose mass is known with reasonable precision. Orosz & Bailyn (1997) observed the source in quiescence, modelled the observed ellipsoidal variations using both quiescent and outburst data and obtained a mass $`M=7.02\pm 0.22M_{}`$. More recently, Shabaz et al. (1999) used radial velocity observations in quiescence to constrain the black hole mass to the range $`5.57.9M_{}`$. Since these results were obtained using only quiescent data they should be free from systematic effects due to X-ray heating of the secondary which can bias the mass determination (Phillips, Shahbaz & Podsiadlowski 1999). These mass constraints make GRO J1655-40 arguably the best object in which to test different models for the high frequency X-ray oscillations and thus to learn how to infer fundamental properties of black holes from X-ray variability measurements.

In this Letter we report the discovery in archival RXTE data of a strong and narrow $`450`$ Hz QPO in the hard X-ray flux from GRO J1655-40. We show that this QPO is present at the same time as a previously discovered 300 Hz QPO (Remillard et al. 1999), however, the 450 Hz QPO is not always seen when the 300 Hz QPO is detected and vice versa. This result is the first detection of a pair of high frequency QPO in a black hole system, and is also the highest frequency modulation yet seen from a black hole. We discuss the implications of our findings for models of X-ray variability in neutron star and black holes systems as well as the implications for black hole spin in GRO J1655-40.

## 2 Data Analysis

GRO J1655-40 is one of the best studied Galactic ‘microquasars’, a class of black hole binaries which produce superluminal radio jets and correlated X-ray activity (Mirabel & Rodriquez; Hjellming & Rupen 1994). The source was discovered by BATSE in July, 1994 (Harmon et al. 1995). A new outburst of GRO J1655-40 was discovered in April, 1996 by the All-Sky Monitor (ASM) aboard RXTE (Remillard et al. 1996). The outburst lasted more than 450 days and an extensive set of RXTE pointed observations were conducted. Remillard et al. (1999) reported an extensive timing analysis of 52 RXTE observations. They reported four types of QPO behavior spanning the range from 0.1 - 300 Hz. They found that the 300 Hz QPO as well as a 0.1 Hz QPO appeared only when the X-ray spectrum was dominated by a hard power law component.

We began investigating the X-ray timing properties of GRO J1655-40 as part of an effort to compare systematically the variability properties of black holes and neutron stars in the RXTE archive. The data we discuss here were obtained from August through November, 1996 as part of proposal P10255 (PI: Remillard) and are now public. The data modes employed for most of these observations included a high time resolution (sampling rate of 65,536 Hz) event mode recording events above PCA channel 35 ($`13`$ keV) and a single bit mode covering the lower energies in a single band from 2 - 12 keV. As part of our analysis we computed power spectra for each observation and data mode separately. We used 32 s intervals and 8 Hz frequency resolution in the range from 1 - 2048 Hz. We then averaged the individual spectra for each observation. We found almost immediately that some of the average power spectra computed using only the high energy event mode data contained a significant feature near 450 Hz. Figure 1 shows the power spectrum for observation 10255-01-10-00 (UT date, 09/09/96, Observation 13 in Table 2 of Remillard et al. 1999) in which we first saw the 450 Hz feature.

To estimate the significance of this feature we first rescaled the power spectrum so that the local mean near 450 Hz was 2 (the value expected for a Poisson process), and then computed the probability of obtaining a power $`P=P_{max}\times 512\times 256`$ from a $`\chi ^2`$ distribution with $`2\times 512\times 256`$ degrees of freedom. Here $`P_{max}`$ is the highest power in the QPO feature. We used this $`\chi ^2`$ distribution because we averaged in frequency by a factor of 512 and averaged 256 individual power spectra. This gives a chance probability of $`2\times 10^8`$ for the highest bin within the QPO profile, better than a 5$`\sigma `$ deviation. We note that this is also a conservative estimate since the QPO feature has not been averaged into a single bin at this frequency resolution.

We next modelled the power spectrum using a Gaussian profile and a power law component in the 100 - 1200 Hz range. The QPO is well fit by a Gaussian centered at $`\nu _0=449.3\pm 5`$ Hz and a width of $`\delta \nu =19.6\pm 4`$ Hz. This gives a coherence value $`Q=\nu /\delta \nu =23`$. The rms amplitude in the 13 - 27 keV band is $`4.8\pm 0.6\%`$.

Having found a significant QPO in one of the observations we then searched the remaining set of observations for a similar feature. We detected the 450 Hz QPO in four additional observations (10255-01-06-01 on 8/16/96, 10255-01-07-00 on 8/22/96, 10255-01-09-00 on 9/4/96, and 10255-01-17-00 on 10/27/96). Power spectra from the three observations with the strongest detections are shown in figure 2 (left panel). We modelled the QPO in each observation in which it was detected and found that its centroid did not change significantly over those observations. Thus, to within the precision of the measurements the frequency appears to be stable. Finally, we computed an average power spectrum for all observations which contained the 450 Hz QPO (see figure 2).

## 3 Discussion and Summary

The detection of a 450 Hz QPO in GRO J1655-40 has many important implications and also raises many interesting questions. One of the most important questions bears on the relationship, if any, between the 300 Hz QPO discovered by Remillard et al. (1999) and the 450 Hz QPO reported here? We detected the 450 Hz QPO in three observations in which Remillard et al. (1999) detected the 300 Hz QPO; 10255-01-06-01 on 8/16/96, 10255-01-07-00 on 8/22/16 and 10255-01-17-00 on 10/27/96. Note that for the latter observation, Remillard et al. report a possible occurrence of the 300 Hz QPO. To investigate whether both QPOs can appear simultaneously we computed an average power spectrum by combining these three observations. We computed two power spectra, one using only the 2 - 12 keV single bit data, and the other using just the 13 - 27 keV event data. Our resulting power spectra are shown in figure 3. We clearly detect both the 300 Hz QPO in the 2 - 12 keV band and the 450 Hz QPO in the hard band. We also investigated the average power spectra of each observation individually. Both the 300 Hz and 450 Hz QPOs are detected in the 8/16/96 and 8/22/96 data. The 450 Hz QPO is detected in the 10/27/96, but the 300 Hz QPO is only marginally detected, consistent with the findings of Remillard et al. (1999). During the 8/22/96 observation we were able to detect the 450 Hz QPO in several subintervals and we confirm that during this observation the 450 Hz frequency did not drift significantly. Although in this observation we can not track the frequency of the 300 Hz QPO with the same temporal resolution as the 450 Hz QPO, it’s stability on timescales of several hours is evident in other observations (for example 10255-01-04-00 on 8/1/96). This evidence, combined with the observed sharpness of both QPOs in the average power spectrum (see figure 3), provides strong evidence that the pair of peaks could not be produced by a single feature which drifted substantially in frequency.

Remillard et al. (1999) report 300 Hz QPO in three additional observations in which we did not detect the 450 Hz QPO. We also investigated the average power spectra of these observations. The results confirm the 300 Hz detections but also show no evidence for the 450 Hz QPO. This indicates that the two QPOs are not always detected together. Interestingly, the pairs of kilohertz QPO seen in neutron star systems are also not always detected at the same time (see van der Klis 2000).

The 300 Hz QPO has a typical amplitude in the 2 - 12 keV band of $`0.8\%`$. Such an amplitude would not be detectable in the 13 - 27 keV data. If the amplitude increased to about $`1.7\%`$ in the hard band then we would have been able to just detect the signal, so we can place a limit on the amplitude of the 300 Hz signal in the 13 - 27 keV band of about $`1.7\%`$. Similarly, if the 450 Hz QPO amplitude stayed the same down into the soft band then we would have detected it easily, so its amplitude must be a strong function of energy, and in fact, it cannot be much larger than about $`0.4\%`$ in the 2 - 12 keV band or it would have been detected. Not only are two high frequency signals present but they appear to have very different energy dependencies as well, perhaps indicating different physical mechanisms or production sites for each oscillation.

The detection of a 450 Hz QPO in GRO J1655-40 has interesting implications for black hole spin in this source. A black hole whose mass lies between $`5.57.9M_{}`$, the 95 % confidence limits from Shahbaz et al. (1999), must have a non-zero angular momentum in order for the orbital frequency at the ISCO to be greater than or equal to 450 Hz (see Bardeen, Press & Teukolsky 1972 for a discussion of the ISCO radius). In figure 4a (left) we have plotted the ISCO radii as a function of dimensionless angular momentum $`j=cJ/GM^2`$ using the mass limits for GRO J1655-40 from Shahbaz et al. (1999) and Orosz & Bailyn (1997). For each mass limit we also show the corresponding radii at which the orbital frequency would be 450 Hz. Even at the lower mass limit of Shahbaz et al. (1999) one requires $`j0.15`$ in order for the orbital frequency at the ISCO to be 450 Hz. Indeed, if the 450 Hz QPO is set by the orbital frequency at the ISCO, then $`j`$ could range from about 0.15 to 0.5. Only with $`M4.8M_{}`$ would a non-rotating hole be able to produce a 450 Hz orbital frequency at the ISCO.

Although the 450 Hz QPO may not actually be produced by Kepler motion of matter at some radius, the Kepler frequency is generally the highest characteristic variability frequency at a given radius. For example, the lowest order diskoseismic modes discussed by Perez, Silbergleit & Wagoner (1999), some of which have been suggested as the mechanism for high frequency black hole QPOs (Nowak et al. 1997), all have characteristic frequencies below that of the Kepler frequency at the ISCO. Therefore, if one of these modes is responsible for the 450 Hz oscillation the conclusion that the black hole has appreciable spin is unchanged, in fact, in such a case it would have to be spinning even faster than suggested by the orbital frequency at the ISCO. Moreover, Cui, Zhang, & Chen (1998) suggested that the 300 Hz QPO in GRO J1655-40 could be associated with the Lense-Thirring precession frequency close to the ISCO radius. For GRO J1655-40 this could only be achieved with a nearly maximally rotating black hole. With the discovery of a 450 Hz QPO it is now even less obvious which frequency should be associated with Lense-Thirring precession (see Mendez, Belloni & van der Klis 1998). If 300 Hz is indeed the precession frequency, then the 450 Hz QPO cannot be due to Kepler motion at the same radius as this would have a much higher frequency. However, if either of the QPO can be associated with Lense-Thirring precession the black hole must be nearly maximally rotating. Although controversial, spectral analysis by Zhang, Cui & Chen (1997) suggests that the inner accretion disk of GRO J1655-40 is hot and compact, perhaps requiring rapid black hole spin (see however, Merloni, Fabian & Ross 2000; Sobczak et al. 1999).

There has been much recent work directed at investigating the relationship between high frequency QPOs and characteristic noise frequencies in neutron stars and black holes. For example, Psaltis et al. (1999) showed that a correlation could be found between a pair of QPO frequencies across a broad range of source types and luminosities, including both neutron star and black hole systems. In the neutron star systems the pair of frequencies are the lower frequency kilohertz QPO and the 20 - 50 Hz QPO most commonly referred to as a horizontal branch oscillation (HBO). For the black hole systems the situation is a bit less certain because of the lack of high frequency features with properties unambiguously similar to the neutron star kHz QPO. In the so called relativistic precession models the observed QPO frequencies have been identified with the Keplerian, the periastron precession and nodal precession frequencies at some characteristic radius in the accretion disk (see Stella, Vietri & Morsink 1999; Psaltis & Norman 2000; Markovic 2000; Markovic & Lamb 2000). Psaltis et al. (1999) tentatively identified the 300 Hz QPO from GRO J1655-40 with the lower kHz QPO in neutron stars. Our detection of a higher frequency QPO in GRO J1655-40 (perhaps the analog to the upper kHz QPO in neutron stars?) would at first glance seem to strengthen this association. To test this we investigated whether the pair of high frequency QPO and the lower frequency ($`18`$ Hz) QPO could be consistently associated with the Keplerian, periastron precession and nodal precession frequencies of a $`7M_{}`$ black hole. We show in figure 4b a plot of the radial epicyclic frequency vs the Kepler frequency for the mass range appropriate for GRO J1655-40 and several different values of the dimensionless angular momentum $`j`$ (curved traces). In the relativistic models the frequency difference between the pair of high frequency QPOs is identified with the radial epicyclic frequency. We also plot the nodal precession frequency for the same set of angular momenta and mass range. The QPO data from GRO J1655-40 are shown with diamond symbols. A black hole with $`0.4<j<0.6`$ could account for the high frequency QPOs, but then the predicted nodal precession frequency is uncomfortably high to be associated with the $`18`$ Hz QPO. Whether or not hydrodynamic corrections can mitigate this apparent discrepancy remains to be seen (see Psaltis 2000).

Models for the Khz QPO in neutron star systems generally fall into two broad classes; beat frequency models (Miller, Lamb, & Psaltis 1998; Strohmayer et al. 1996), which require the spinning neutron star surface to produce a pair of frequencies; and relativistic disk models (Stella, Vietri, & Morsink 1999; Psaltis & Norman 2000; Nowak et al. 1997; Perez et al. 2000), which generate the frequencies in the disk alone and are therefore thought to be more generally applicable to both neutron star and black hole systems. Our discovery of a 2nd high frequency QPO in GRO J1655-40 proves that it is not necessary to have a hard surface (as in a neutron star) to produce a pair of high frequency peaks. This would seem to provide some support to the relativistic disk models, since this model can naturally produce pairs of peaks in either class of source, however, a straightforward application of the model seems to have some difficulty (see above). Caution regarding the identification of the QPOs in GRO J1655-40 with the better understood kHz QPO in neutron stars is warranted because of the different properties of the QPOs in GRO J1655-40 and the neutron star kHz QPOs.

###### Acknowledgements.

We thank Craig Markwardt, Cole Miller and Jean Swank for many helpful discussions and comments on the manuscript. We thank the referee, Ron Remillard, for his comments which helped us improve the paper.

## References

Bardeen, J. M., Press, W. H. & Teukolsky, S. A. 1972, ApJ, 178, 347

Cui, W., Zhang, S. N. & Chen, W. 1998, ApJ, 492, L53

Harmon, B. A. et al. 1995, Nature, 374, 703

Hjellming, R. M. & Rupen, M. P. 1995, Nature, 375, 464

Homan, J. et al. 2000, ApJ in press, (astro-ph/0001163)

Leahy, D. A. et al. 1983, ApJ, 266, 160

Markovic, D. 2000, MNRAS, submitted, (astro-ph/0009450)

Markovic, D. & Lamb, F. K. 2000, MNRAS, submitted, (astro-ph/0009169)

Markwardt, C. B., Swank, J. H. & Taam, R. E. 1999, ApJ, 513, L37

Mendez, M., Belloni, T., & van der Klis, M. 1998, ApJ, 499, L187

Merloni, A., Fabian, A. C. & Ross, R. R. 2000, MNRAS, 313, 193

Merloni, A., Vietri, M., Stella, L. & Bini, D. 1999, MNRAS, 304, 155

Miller, M. C., Lamb, F. K. & Psaltis, D. 1998, ApJ, 508, 791

Milsom, J. A., & Taam, R. E. 1997, MNRAS, 286, 359

Mirabel, I. F. & Rodriguez, L. F. 1994, Nature, 371, 46

Morgan, E. H., Remillard, R. A. & Greiner, J. 1997, ApJ, 482, 993

Nowak, M. A., Wagoner, R. V. Begelman, M. C. & Lehr, D. E. 1997, ApJ, 477, L91

Orosz, J. A. & Bailyn, C. D. 1997, ApJ, 477, 876

Perez, C. A., Silbergleit, A. S., Wagoner, R. V. & Lehr, D. E. 1997, ApJ, 476, 589

Phillips, S. N., Shahbaz, T. & Podsiadlowski, P. 1999, MNRAS, 304, 839

Psaltis, D. 2000 ApJ, submitted, (astro-ph/0010316)

Psaltis, D., Belloni, T. & van der Klis, M. 1999, ApJ, 520, 262

Psaltis, D. & Norman, C. 2000, ApJ, in press (astro-ph/0001391)

Remillard, R. A., et al. 1996, IAU Circ. 6393

Remillard, R. A., Morgan, E. H., McClintock, J. E., Bailyn, C. D. & Orosz, J. A. 1999, ApJ, 522, 397

Remillard, R. A. & Morgan, E. H. 1998, in The Active X-ray Sky, ed. L. Scarsi, H. Bradt, P. Giommi, & F. Fiore (Amsterdam: Elsevier), 316

Remillard, R. A., Morgan, E. H., McClintock, J. E., Bailyn, C. D. & Orisz, J. A. 1998, in Proc. 18th Texas Symp. on Relativistic Astrophysics, ed. A. Olinto, J. Frieman, & D Schramm (Singapore: World Scientific), 750

Shahbaz, T., van der Hooft, F., Casares, J., Charles, P. A. & van Paradijs, J. 1999, MNRAS, 306, 89

Sobczak, G. J., McClintock, J. E., Remillard, R. A., Bailyn, C. D. & Orosz, J. A. 1999, ApJ, 520, 776

Stella, L., Vietri, M. & Morsink, S. M. 1999, ApJ, 524, L63

Strohmayer, T. E., Zhang, W., Swank, J. H., Smale, A. P., Titarchuk, L., Day, C. & Lee, U. 1996, ApJ, 469, L9

van der Klis, M. 2000, ARAA, in press (astro-ph/0001167)

Zhang, S. N., Cui, W. & Chen, W. 1997, ApJ, 482, L155

## 4 Figure Captions

Figure 1: Average power spectrum for observation 10255-01-10-00 in the 13 - 27 keV energy band. The spectrum was computed by averaging a total of 256 individual spectra each computed from a 32 second interval. The frequency resolution is 16 Hz and the spectrum has been normalized following Leahy et al. (1983). The QPO at 450 Hz is clearly visible.

Figure 2: Average power spectra from three observations with the most significant detections of the 450 Hz QPO (top three traces). The spectra were computed in the same way as for figure 1. The centroid of the 450 Hz peak did not change significantly over these observations. The spectra have been shifted vertically for clarity. The average power spectrum at 2 Hz frequency resolution of all observations in which the 450 Hz QPO was detected is shown in the bottom trace. The lower frequency $`18`$ Hz QPO can also be seen.

Figure 3: Power spectra computed from the three observations in which we detected the 450 Hz QPO and in which Remillard et al. reported either detections or evidence of the 300 Hz QPO. The upper spectrum was computed from the 2 - 12 keV single bit mode data only while the lower spectrum was computed using only the 13 - 27 keV event mode data. Both the 300 and 450 Hz QPO are clearly detected simultaneously.

Figure 4: The radius of the inner most stable circular orbit (ISCO) as a function of dimensionless black hole angular momentum, $`j`$, for the mass limits for GRO J1655-40 of Shahbaz et al. (1999) and Orosz & Bailyn (1997) (left panel). The curves denoting the ISCO radii are labelled with the corresponding mass (5.5, 6.8, 7.2, or 7.9 $`M_{}`$). Also shown for each mass are the corresponding radii at which the orbital frequency is equal to 450 Hz. Curves corresponding to a particular mass are plotted with the same line style (ie. solid, dashed, etc.). For example, for a 5.5 $`M_{}`$ black hole one requires $`j0.15`$ in order for the Kepler frequency at the ISCO to be 450 Hz. The Radial epicyclic frequency (curved lines) and the nodal precession frequency (diagonal lines) vs Keplerian frequency for the Shahbaz et al. mass limits for GRO J1655-40 and for several different values of the dimensionless angular momentum $`j`$ (right panel). the solid, dotted and dashed curves correspond to $`j=0.2,\mathrm{\hspace{0.33em}0.4}`$, and 0.6, respectively. The curves for $`j=0.2`$ are labelled with their corresponding masses in order to indicate the sense of the mass dependence. The difference frequency between the 300 and 450 Hz QPO and the low frequency ($`18`$ Hz) QPO are plotted with diamonds. Note that the error bars are smaller than the symbol size.


format: text

text:
# 1 Model

## 1 Model

Consider $`N+1`$ spins positioned at the sites of crystalline ferromagnet lattice. Let us suppose that one of the site is located at the origin of the coordinates. The probability of finding the interaction field $`𝐇`$ to lie in the interval $`(𝐇,𝐇+d𝐇)`$ at the origin is

$$\delta \left[𝐇\underset{k=1}{\overset{N}{}}𝐡_k(𝐦_k,𝐫_k)\right]d𝐇,$$

where $`𝐡_k(𝐦_k,𝐫_k)`$ is interaction field induced by the atom with the magnetic moment $`𝐦_k`$ situated at the site $`𝐫_k`$, $`\delta `$ is Dirac’s delta function. Then the random field distribution function will be

$$W(𝐇)d𝐇=\delta \left[𝐇\underset{k=1}{\overset{N}{}}𝐡_k(𝐦_k,𝐫_k)\right]\underset{k=1}{\overset{N}{}}\tau _k\left(𝐦_k\right)d𝐦_kd𝐇,$$
(1)
where $`\tau _k\left(𝐦_k\right)`$ is the distribution over atom magnetic moments $`𝐦_k`$ function. In the approach of the Ising model, when the magnetic moments are aligned with the selected direction given by unite vector $`𝐬`$ and their magnitudes are equal to $`m_0`$, $`\tau _k`$ takes the following form

$$\tau _k\left(𝐦_k\right)=\alpha _k\delta \left(m_0𝐬𝐦_k\right)+\beta _k\delta \left(m_0𝐬+𝐦_k\right),$$

where $`\alpha _k`$ and $`\beta _k`$ are the probabilities of finding magnetic moment to be oriented along and opposite to chosen direction. Let us consider a case where the interaction field $`𝐡_k(𝐦_k,𝐫_k)`$ is a sum of exchange and dipole-dipole interaction:

$$𝐡_k(𝐦_k,𝐫_k)=J(r_k)\frac{𝐦_k}{m_k}+\frac{3(𝐫_k𝐦_k)𝐫_kr_k^2𝐦_k}{r_k^5},$$
(2)
where $`J(r_k)=J`$ being exchange constant if $`k`$th atom is nearest neighbor and $`J(r_k)=0`$ otherwise.

Let us replace $`\alpha _k`$ and $`\beta _k`$ by their average values $`\overline{\alpha }`$ and $`\overline{\beta }`$ which are given by the equation

$$\overline{\alpha }=\frac{\mathrm{exp}\left(m_0𝐬𝐇/k_BT\right)}{2\mathrm{cosh}\left(m_0𝐬𝐇/k_BT\right)}W(𝐇)𝑑𝐇,\overline{\alpha }+\overline{\beta }=1,$$
(3)
where $`k_B`$ is Boltzmann constant, $`T`$ is temperature.

The characteristic function $`A(𝝆)=W(𝐇)\mathrm{exp}(i𝝆𝐇)𝑑𝐇`$, being the Fourier transform of the function $`W(𝐇)`$, takes the form

$$A(𝝆)=\underset{k}{}\left[\overline{\alpha }\mathrm{exp}\left(i𝝆𝐡_k(m_0𝐬,𝐫_k)\right)+\overline{\beta }\mathrm{exp}\left(i𝝆𝐡_k(m_0𝐬,𝐫_k)\right)\right].$$
(4)
After the inverse transformation the distribution function is

$$W(𝐇)=\underset{j_1=0}{\overset{1}{}}C_1^{j_1}\underset{j_2=0}{\overset{1}{}}C_1^{j_2}\mathrm{}\underset{j_N=0}{\overset{1}{}}C_1^{j_N}\overline{\alpha }^{_{k=1}^Nj_k}\overline{\beta }^{N_{k=1}^Nj_k}\delta \left[𝐇\underset{k=1}{\overset{N}{}}(2j_k1)𝐡_k(m_0𝐬,𝐫_k)\right],$$
(5)
where $`C_1^{j_k}=1`$ are binomial coefficients.

One can obtain the free energy per atom using relation $`F=k_BT\mathrm{ln}\mathrm{exp}\left(\frac{m_0𝐬𝐇}{k_BT}\right)W(𝐇)𝑑𝐇`$. After integrating $`F`$ will be

$$F=k_BT\underset{k}{}\mathrm{ln}\left[\overline{\alpha }\mathrm{exp}\left(\frac{m_0𝐬𝐡_k(m_0𝐬,𝐫_k)}{k_BT}\right)+\overline{\beta }\mathrm{exp}\left(\frac{m_0𝐬𝐡_k(m_0𝐬,𝐫_k)}{k_BT}\right)\right]$$
(6)
and substituting (2) into Eq. (6) $`F`$ takes the form

$$F=k_BT\underset{k}{}\mathrm{ln}\left\{\overline{\alpha }\mathrm{exp}\left[\frac{1}{k_BT}\left(J(r_k)+m_0^2\frac{3(𝐫_k𝐬)^2r_k^2}{r_k^5}\right)\right]+\overline{\beta }\mathrm{exp}\left[\frac{1}{k_BT}\left(J(r_k)+m_0^2\frac{3(𝐫_k𝐬)^2r_k^2}{r_k^5}\right)\right]\right\}.$$
(7)
Note that the terms describing contribution of the external field $`𝐇_{\mathrm{ext}}`$ and single-ion anisotropy $`m_0𝐬𝐇_{\mathrm{ext}}+D(s_x^4+s_y^4+s_z^4)`$, $`D`$ being single-ion cubic anisotropy constant, may be added to the Eq. (7). One can calculate $`\overline{\alpha }`$ and $`\overline{\beta }`$ using equation derived from (3)  taking account of the exchange coupling only:

$$\overline{m}=2\overline{\alpha }1=\underset{j=0}{\overset{Z}{}}C_Z^j\overline{\alpha }^j(1\overline{\alpha })^{Zj}\mathrm{tanh}\left[\frac{J}{k_BT}(2jZ)\right],$$
(8)
where $`\overline{m}=\overline{\alpha }\overline{\beta }`$, $`Z`$ is number of nearest neighbors.

## 2 Results

Let us treat dipole interaction as small parameter with respect to exchange coupling ($`m_0^2J`$), also consider the case of high temperatures ($`m_0^2k_BT`$). Expand $`F`$ into a series in powers of $`m_0^2`$ and direction cosines $`s_x,s_y,s_z`$. Taking the identity for cubic crystals

$$s_x^4+s_y^4+s_z^4=12\left(s_x^2s_y^2+s_y^2s_z^2+s_z^2s_x^2\right)$$

into account and retaining first nonzero anisotropic terms of this expansion we will obtain the anisotropic part of free energy

$$F_{\mathrm{an}}=36\overline{\alpha }\overline{\beta }\frac{m_0^4}{k_BT}\frac{1}{(a/2)^6}\left[3L_1L_2+L_3\frac{1}{\overline{\alpha }\mathrm{exp}\left(J/k_BT\right)+\overline{\beta }\mathrm{exp}\left(J/k_BT\right)}\right]\left(s_x^2s_y^2+s_y^2s_z^2+s_z^2s_x^2\right),$$
(9)
$`a`$ being lattice constant,

$$L_1=\underset{\begin{array}{c}j,k,l=\mathrm{}\\ j^2+k^2+l^24\end{array}}{\overset{\mathrm{}}{}}\frac{j^4}{(j^2+k^2+l^2)^5},$$
(10)
$$L_2=\underset{\begin{array}{c}j,k,l=\mathrm{}\\ j^2+k^2+l^24\end{array}}{\overset{\mathrm{}}{}}\frac{j^2k^2}{(j^2+k^2+l^2)^5},$$
(11)
can be calculated numerically. For bcc lattice it was calculated that $`L_10.00454`$, $`L_20.04343`$, $`L_3=\frac{16}{243}`$ and for fcc lattice $`L_10.03200`$, $`L_20.12912`$, $`L_3=\frac{1}{8}`$.

As one can see at $`T=0`$ ferromagnet is totally ordered and $`F_{\mathrm{an}}`$ derived from Eq. (7) is equal to zero since quadrupole moment of this system is zero. With increase in temperature $`|F_{\mathrm{an}}|`$ raises and at $`T>T_c`$ $`F_{\mathrm{an}}0`$ because of used approximation of Ising model (magnetic moments of all atoms are parallel to chosen direction).

The method proposed can be used for estimation of the magnetic anisotropy induced by dipole-dipole interactions in two-dimensional and bulk materials.

This work was supported by Russian Ministry of Education. The authors would like to thank V. I. Belokon’ for fruitful discussion.


text:
# Symmetry-Breaking and Chaos in Electron Transport in Semiconductor Superlattices

## I Introduction

Since the seminal work of Esaki and Tsu  predicting negative differential conductivity in the current-voltage characteristics of a semiconductor superlattice (SSL), SSLs have been the subject of intense theoretical and experimental study. SSLs consist of thin layers of alternating semiconductors – typically 10-100 layers ranging from 10 to 500 Å in thickness – whose different band structures create a periodic potential along the growth axis, essentially forming a one-dimensional, artificial crystal. The superlattice period can be controlled during the growth process and is generally much larger than the period of the constituent materials. Consequently the conduction and valence bands break into narrow allowed and forbidden regions termed minibands, with typical widths ranging from 1 to 100 meV. For an introduction to SSL design and properties, see reference .

Recent experiments on SSLs in the presence of high-frequency fields have revealed many novel, nonlinear effects, including photon-assisted tunneling , a reduction of the current during laser irradiation , absolute negative conductivity , resonant increases in the current , and, most recently, deterministic chaos in the sequential resonant tunneling regime . Our previous work  suggests that chaos should also occur in the miniband regime for a SSL driven by an AC field. In this brief note, we recall some of the principal results of that work and report our current studies of symmetry-breaking effects in electron transport. Using numerical simulations we establish the possibility of spontaneous DC current and voltage generation in response to a purely AC drive. Our results cover the whole range of frequency/amplitude of the applied field, extending a recent investigation of the high frequency, strong field limit .

## II Balance Equations

We consider electrons subject to an external field $`E_{ext}=E_0\mathrm{cos}(\mathrm{\Omega }t)`$ using the tight-binding dispersion relation, $`ϵ(k)=\frac{\mathrm{\Delta }}{2}[1\mathrm{cos}(ka)]`$, for SSL period $`a`$, minibandwidth $`\mathrm{\Delta }`$, electron energy $`ϵ`$ and crystal momentum $`k`$. We assume that the frequency of the external field is not large enough to support Zener tunneling to higher minibands, and that the electric field is uniform across the SSL. To model the electron dynamics, we generalize the semiclassical balance equations of Ignatov et al.  by including the self-consistent electric field resulting from electron motion, obtaining ,

$`\dot{V}`$ $`=`$ $`eE_{tot}(t)/m(\epsilon )\gamma _vV`$ (1)
$`\dot{\epsilon }`$ $`=`$ $`eE_{tot}(t)V\gamma _\epsilon (\epsilon \epsilon _0)`$ (2)
$`\dot{E}_{sc}`$ $`=`$ $`4\pi j/\kappa \alpha E_{sc}`$ (3)
where $`E_{tot}=E_{ext}+E_{sc}`$ is the total electric field, including the self-consistent contribution; $`j=eNV`$ is the current density; $`V(t)`$ and $`\epsilon (t)`$ are the average electron velocity and energy, respectively; $`\gamma _v`$, $`\gamma _\epsilon `$, and $`\alpha `$ are phenomenological relaxation rates for the average velocity, average energy and self-consistent field, respectively; $`\kappa `$ is the average dielectric constant for the SSL; $`\epsilon _0`$ is the temperature-dependent equilibrium energy; $`N`$ is the carrier concentration, assumed uniform across the SSL; and $`m(\epsilon )=\frac{m_0}{12\epsilon /\mathrm{\Delta }},m_0=\frac{2\mathrm{}}{\mathrm{\Delta }a^2}`$, is the effective mass of the electrons.

The balance equations remain invariant under the symmetry transformation:

$$tt+\frac{2n+1}{2}T,VV,E_{sc}E_{sc},$$
(4)
where $`T`$ is the period of the applied field and $`n=0,1,2,\mathrm{}`$; in other words, if $`\{V(t),\epsilon (t),E_{sc}(t)\}`$ satisfy (1)–(3), then so do $`\{V(t+\frac{2n+1}{2}T),\epsilon (t+\frac{2n+1}{2}T),E_{sc}(t+\frac{2n+1}{2}T)\}`$. Physically, this represents the lack of preferred spatial direction for a SSL without DC bias; under this transformation the electron velocity and the self-consistent field reverse direction, but the electron energy remains unaltered.

Symmetric solutions also remain invariant under transformation (4); hence $`V(t+\frac{2n+1}{2}T)=V(t)`$ and $`E_{sc}(t+\frac{2n+1}{2}T)=E_{sc}(t)`$. For such cases the time-averages of both the velocity and the self-consistent field are zero. However there may also exist solutions which do not satisfy this symmetry, for which $`V`$ and $`E_{sc}`$ have non-zero average values; these symmetry-broken solutions correspond to the spontaneous generation of a DC current and bias.

In the limiting case $`\gamma _v=\gamma _\epsilon =0`$, (1)–(3) can be transformed to the equation describing the Resistively Shunted Junction (RSJ) model of a Josephson junction (JJ)  This model possesses a thoroughly studied symmetry analagous to (4) . When biased by an AC current of frequency $`\mathrm{\Omega }`$, the JJ may develop a DC component to its voltage difference satisfying the phase-locking condition $`Un\mathrm{\Omega }`$, where $`n`$ is an integer, $`U`$ is the voltage across the JJ, and angular brackets indicate time-averages .

## III Numerical Studies

In our earlier work , we showed that the balance equations give rise to chaotic behavior for a wide range of damping and driving parameters. Readers should consult this reference for a full description of the numerics and scalings used in the simulations. Briefly, however, there are three important frequencies: the drive frequency, $`\mathrm{\Omega }`$; the Stark frequency of the driving field, $`\omega _s=eaE_0/\mathrm{}`$; and the frequency of the cooperative oscillations associated with the self-consistent field, $`\omega _E=\sqrt{2\pi e^2Na^2\mathrm{\Delta }/\kappa \mathrm{}^2}`$ which is usually in the terahertz range. As it depends only on SSL material parameters, we use $`\omega _E`$ as the scale of inverse time, and throughout the rest of this article we consider the dimensionless frequencies $`\mathrm{\Omega }/\omega _E`$, $`\omega _s/\omega _E`$, $`\gamma _v/\omega _E`$, $`\gamma _\epsilon /\omega _E`$ and $`\alpha /\omega _E`$ without explicitly showing the scaling.

As in many studies of the JJ , we examine the dynamics for fixed damping while varying the parameters of the applied field. Figures 1a and 1b indicate the region of chaos and symmetry-breaking, respectively. Note the similar general shape of these regions; further work shows that for greater damping the region of chaos diminishes while the region of symmetry-breaking retains this general shape .

Fig. 1c plots the time-averaged self-consistent field, $`eaE_{sc}/\mathrm{}\omega _E`$, as a function of external frequency, $`\omega _s`$, for the frequency $`\mathrm{\Omega }=0.9`$. Phase-locked symmetry-breaking where $`eaE_{sc}/\mathrm{}\omega _En\mathrm{\Omega }/\omega _E`$ for $`n=\pm 1,\pm 2`$ is clearly visible. A careful examination reveals that $`eaE_{sc}/\mathrm{}\omega _E`$ also depends weakly on $`\omega _s`$; therefore the phase-locking relation is only approximately satisfied.

From (3) it immediately follows that $`VE_{sc}`$ for non-chaotic solutions. Consequently symmetry-breaking should be experimentally observable as a DC current and bias across an unbiased SSL illuminated by a terahertz laser. For representative SSL parameters of $`a=100\AA `$, $`N=3\times 10^{15}cm^3`$, $`\mathrm{\Delta }=22meV`$, and $`ϵ=13`$ , $`\omega _E=3.3\times 10^{12}s^1`$. Taking $`eaE_{sc}/\mathrm{}\omega _E_{max}=1.5`$ yields a DC bias per SSL period of $`3mV`$. Furthermore, a lateral area of $`8\mu m^2`$ would result in a current of $`1\mu A`$. Both these magnitudes should be experimentally detectable.

## IV Conclusion

We have studied the complex dynamics of semiclassical balance equations that model the motion of electrons in the miniband of a semiconductor superlattice exposed to an external ac field. The equations include the self-consistent field created by electron motion. Numerical simulations indicate the possibility of broken-symmetry solutions, which occur in roughly the same laser frequency and electric field ranges as chaotic dynamics. The symmetry-breaking engenders a DC bias and current that could be experimentally detectable.

We would like to contrast this work with recent theoretical  and experimental  work predicting and observing chaos for SSLs in the incoherent tunneling regime when the electric field breaks into domains of different magnitudes . Our model considers the coherent motion of electrons within a miniband in a spatially homogenous electric field. A complete theoretical understanding of when the applied field breaks into domains has not yet been attained, so it is encouraging that chaos is predicted in quite different models applicable in different regimes of SSL design and applied field characteristics.

We are especially grateful to Gennady Berman for his advice and continuing collaboration on many aspects of chaos in mesoscopic systems. We thank Sergei Turovets for a discussion on the notation of symmetry-breaking in dynamical systems. K.N.A. thanks the Department of Physics at The University of Illinois at Urbana-Champaign for hospitality. This work was partially supported by Linkage Grant No. 93-1602 from the NATO Special Programme Panel on Nanotechnology. K.N.A. was supported by INTAS(94-2058) and KRSF(6F0030). E.H.C. thanks the U.S. Department of Education for support by a GAANN Fellowship (DE-P200A40532), and J.C.M. thanks the U.S.-NSF for support under its REU program (NSF Grant No. PHYS93-22320).


text:
# Structural phase transition at high temperatures in solid molecular hydrogen and deuterium

## I INTRODUCTION

The possibility for hydrogen to undergo a pressure-induced phase transition from a proton-paired insulator to a monatomic metal was first suggested in 1935.  Since then a large body of experimental and theoretical works have been done to determine its phase diagram in the pressure ($`P`$) and temperature ($`T`$) plane. Despite studies at $`P`$ reported to be as high as 342 GPa,  the tenacious covalent-bond feature of a pair of protons does not allow observation of a monatomic phase so far.

For $`P`$ up to ca. 200GPa, the phase diagram has been investigated by several research groups at the room temperature and below, based on the optical measurement performed in the diamond anvil cell (DAC) devices. By now it is well established that the solid hydrogen exhibits at least three different molecular phases (the phases I$``$III), although some details are still in dispute.  (1) Phase I: At $`P<110`$GPa, each center of the H2 molecule occupies the lattice site of an hcp structure, but quantum rotational effects overcome librational barriers, leading to an orientationally disordered phase. (2) Phase II or broken-symmetry phase (BSP): At $`P`$ between 110GPa and 150GPa, anisotropic intermolecular interactions freeze the molecular rotations into an orientationally ordered phase. (3) Phase III or H-A phase: At $`P`$ above 150 GPa, a third phase is observed, expected to be another kind of orientationally ordered phase.

In 1996, Weir and coworkers gave one of the exciting results. They observed high electrical conductivity in the shock compressed H2 and D2 liquids and interpreted it as a transition from a semiconducting to a metallic diatomic fluid at 140 GPa and 3000K.  This experiment clearly demonstrates the importance of temperature effects in the search of metallization. For $`T`$ larger than 5000K, the effects are studied theoretically in the molecular, the dissociated, and the plasma regime of dense H2.  It is, however, still largely unclear how $`T`$ influences the state of solid H2 at $`T`$ higher than the room temperature. We need to know the phase diagram at this temperature region to connect the liquid state with the low-temperature phases for a comprehensive understanding of dense H2. Thus we focus our attention on the range of $`P`$ up to 200GPa and 300K$`<T<1000`$K in this paper.

In order to study the condensed H2 phases theoretically, several methods have been adopted at various levels of approximations to the ab initio Hamiltonian representing the coupled system of $`N_a`$ protons and $`N_a`$ electrons. These include the calculation of electronic energies in the local density approximation (LDA) or its refinements to the density functional theory for a variety of crystal structures with molecular orientations fixed to some particular configuration,  and the implementation of ab initio molecular dynamics treating protons as classical particles. We should note, however, that the strong quantum nature of light protons requires a more careful quantum-mechanical description of their zero-point motion. In fact, using either the first-principle path-integral molecular dynamics  or the quantum Monte Carlo (QMC) simulations,  the calculation treating protons as dynamic quantum particles has already been performed. Although desirable, this approach demands the computer resources very much. Thus, even in a recent paper,  the calculation is done only with as small as $`N_a=64`$. This size of $`N_a`$ is too small in order to obtain a correct equilibrium structure free from any bias of an initially assumed one through simulations.

A ten times increase in $`N_a`$ is possible if we approach the problem by adopting the hydrogen molecule as a basic ingredient in simulations rather than the proton-electron mixture. In this approach, the system is reduced to a quantum-mechanical problem of $`N`$ (=2$`N_a`$) molecules interacting to each other through an effective intermolecular pair potential $`V_{\mathrm{pair}}`$. Thus the QMC calculation is implemented only for the nuclear degrees of freedom  and the electronic degrees of freedom are implicitly taken into account in the choice of $`V_{\mathrm{pair}}`$. In the present paper, we adopt this approach and employ the finite-temperature path-integral Monte Carlo (PIMC) method to investigate both lattice and orientational transitions in the molecular phase with zero-point motions incorporated rigorously. We shall use an empirically determined $`V_{\mathrm{pair}}`$ which is the sum of the Hemley-corrected Silvera-Goldman potential  and the Runge-scaled Shaefer potential.  This choice of $`V_{\mathrm{pair}}`$ is known to give an equation of state in very good agreement with experiment as well as the I/II phase boundaries for both H2 and D2 over the whole experimentally investigated range of pressures.

This paper is organized as follows: We shall explain our theoretical model in more detail in Sec. II, including a review of trials to determine $`V_{\mathrm{pair}}`$, a description of our choice of $`V_{\mathrm{pair}}`$, a brief summary of the PIMC method with a constant-pressure ensemble, and some computational details. In Sec. III our results are shown for both solid H2 and D2, indicating a temperature-induced solid-solid structural phase transition. The structure of the new high-temperature phase, the solid-liquid phase transition, and the phase diagram are also given here. Section IV contains a discussion on the role of zero-point motion of protons in the solid-solid structural phase transition and finally in Sec. V we summarize our results.

## II THEORETICAL MODEL

### A Molecule base

Within the temperature and pressure range that a molecule can be regarded as a basic ingredient, a quantum molecular solid can be described by the Hamiltonian $`H_{\mathrm{nuclear}}`$ as the sum of $`T_{\mathrm{nuclear}}`$ and $`V_{\mathrm{nuclear}}`$, where

$`T_{\mathrm{nuclear}}={\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \underset{i}{\overset{N}{}}}_{R_i}^2+{\displaystyle \frac{\mathrm{}^2}{2I}}{\displaystyle \underset{i}{\overset{N}{}}}𝐋_i^2,`$ (1)
and

$`V_{\mathrm{nuclear}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{\overset{N}{}}}V_{\mathrm{pair}}(𝐑_{ij},𝛀_i,𝛀_j).`$ (2)
Here $`𝐑_i`$ is the center-of-mass position vector of the $`i`$th molecule, $`𝛀_i`$ its orientation vector, $`𝐋_i`$ its angular momentum operator, and the intermolecular separation vector $`𝐑_{ij}𝐑_i𝐑_j`$. The molecular mass and moment of inertia are denoted by $`m`$ and $`I`$, respectively. For a linear molecule with the intramolecular bond length fixed, $`I`$ is a scalar constant. The values of $`m=3676`$ and $`7352`$ (atomic units) for H2 and D2, respectively, and rotational constants $`B\mathrm{}^2/2I=84.98`$ and 42.92 K for H2 and D2, respectively, are used.  Among the the nuclear degrees of freedom, we include both translation and rotation modes of molecules in the kinetic energy $`T_{\mathrm{nuclear}}`$, but we consider implicitly the intramolecular vibration mode along with the electronic degrees of freedom by judiciously choosing $`V_{\mathrm{pair}}`$ the pair-wise sum of which gives the potential energy $`V_{\mathrm{nuclear}}`$.

We have enough evidence to support the present approach in our research range ($`0200`$GPa and $`01000`$K). First, we can neglect the molecular dissociation, because its fraction is less than 0.01% from the same analysis as the one leading to about 5% at 140GPa and 3000K.  Secondly, ab initio quantum-chemistry calculations at $`N`$=2 lead us to the conclusion that the dependence of $`V_{\mathrm{pair}}`$ on the bond length $`r_0`$ is weak, even though $`r_0`$ shrinks by more than 10% from 1.4 bohr (0.74Å) as two H2 molecules come closer.  This assures the validity to employ an $`r_0`$-independent $`V_{\mathrm{pair}}`$. The notion of $`V_{\mathrm{pair}}`$ is relevant if $`r`$ the intermolecular distance is larger than about 2.3 bohr (1.2Å), because each molecule is estimated to possess a repulsive core of the radius of 0.6Å.  Finally the effect of zero-point intramolecular vibrations (vibrons) can be taken into account in the form of $`V_{\mathrm{pair}}`$ as we shall illustrate in Sec. IV by using a toy model.

### B H2-H2 interaction

Rather than determining $`V_{\mathrm{pair}}`$ from first principles, we shall resort to a phenomenological approach which is well reviewed in the literature.  In general, the multi-dimensional function $`V_{\mathrm{pair}}(𝐑_{12},𝛀_1,𝛀_2)`$ can be expressed as the sum of the isotropic pair potential $`V_0(R_{12})`$ and the anisotropic potential $`V_{ani}(𝐑_{12},𝛀_1,𝛀_2)`$. A variety of analytic forms have been proposed for $`V_0(R)`$ and their relevance has been tested in the past. Four of them are illustrated in Fig. 1. Because the Lennard-Jones (LJ) two-parameter form (the dotted-dashed curve in Fig. 1) has an unphysically strong repulsive core, more sophisticated forms have been proposed. The most commonly used form is either BUCK potential $`V_{\mathrm{BUCK}}`$  or SG potential $`V_{\mathrm{SG}}`$,  plotted by dashed and dotted curves, respectively, in Fig. 1. Although they have almost the same analytic form, $`V_{\mathrm{BUCK}}`$ is fitted to the gas-phase data, while $`V_{\mathrm{SG}}`$ to the solid-phase data. In order to include the three-body interactions effectively in the solid environment, a further refinement is introduced to the original form of $`V_{\mathrm{SG}}`$ by incorporating an additional repulsive term, leading to the final form of SG potential. Compared to $`V_{\mathrm{BUCK}}`$, the net effects contributed from the three-body interaction are a hardening of $`V_{\mathrm{SG}}`$ at small $`R`$ and a slight raising of the well depth. Their concrete forms can be found in Refs. , and .

This $`V_{\mathrm{SG}}(R)`$ works well for the solid H2 and D2 under ambient pressures. Under high pressures, however, an additional softening effect found by Hemley et al. is needed to correctly describe the enhanced many-body contributions in the short-range region caused by the dense solid environment. This effect can be incorporated by an ad hoc short-range correction $`V_{\mathrm{SR}}(R)`$ to $`V_{\mathrm{SG}}(R)`$,  leading to the definition of the Hemley-corrected SG potential $`V_{\mathrm{SGH}}(R)`$ (the solid curve in Fig. 1), where the actual form of $`V_{\mathrm{SR}}(R)`$ can be found in Refs.  and . This softening correction at small $`R`$ is very important in reproducing the equation of state (EOS) derived experimentally from the liquid D2 shock-wave data,  as assessed by Cui et al. in Ref.  in which the EOS for D2 at $`T=20`$K and for $`P`$ up to 40GPa is calculated for testing the intermolecular potential. From the test, the EOS is proved to be insensitive to the details of $`V_{ani}`$, but it is very sensitive to the choice of $`V_0(R)`$ which controls the translational motion of molecules. In fact, $`V_{\mathrm{SGH}}(R)`$ gives results in excellent agreement with the experimental data,  while $`V_{\mathrm{SG}}(R)`$ does not. The same conclusion is also drawn at $`T=300`$K from a similar test; for $`P`$ up to 200GPa, $`V_{\mathrm{SGH}}(R)`$ provides the results for the EOS in line with the experimentally obtained fitting formula.

Only a few works have been done for the assessment of $`V_{ani}`$. Using the ab initio quantum chemical method to evaluate the contributions such as the long-range electronic quadrupole-quadrupole interaction and the atom-diatom scattering term, Schaefer et al. have given a potential for $`V_{ani}`$.  Comparing the results obtained from the electronic-structure calculation in the LDA, Runge et al. have found that this Schaefer potential $`V_{\mathrm{Schaefer}}`$ is too repulsive in the dense solid. Thus, in order to soften it, they have introduced a scaling factor to $`V_{\mathrm{Schaefer}}`$ to provide $`V_{ani}`$ as $`\alpha _nV_{\mathrm{Schaefer}}`$.  This factor $`\alpha _n`$ is linear in the nearest neighbor spacing $`R_{nn}`$ to correctly describe the dense-solid environment as $`\alpha _n=0.61+0.31(R_{nn}/R_{nn}^00.5)`$ with $`R_{nn}^0=3.789`$Å, the H2 equilibrium zero-pressure nearest-neighbor spacing. Both the bare $`V_{\mathrm{Schaefer}}`$ and the scaled one for $`V_{ani}`$ have been tested in Ref.  to find that this difference is crucial in reproducing the I/II phase line for the solid D2. The best agreement with the experimental data for $`P`$ up to 150GPa is achieved by using the scaled anisotropic potential combined with $`V_{\mathrm{SGH}}(R)`$ for $`V_0(R)`$. The bare $`V_{\mathrm{Schaefer}}`$ predicts the transition to occure at much lower pressures. This difference was also observed in the fixed-lattice PIMC study.

For those reasons mentioned above, we shall choose the sum of $`V_{\mathrm{SGH}}(R)`$ and the scaled Schaefer potential as $`V_{\mathrm{pair}}`$ in this paper. We believe that this $`V_{\mathrm{pair}}`$ is most reliable for the study of solid molecular H2 and D2 in our research range, in particular, for $`P<150`$GPa.

### C Path Integral Monte Carlo method with constant pressure ensemble

With the Hamiltonian thus provided, we model the bulk solid by a simulation cell, which is periodically duplicated in all three spatial dimensions to minimize surface and finite-size effects. Initial size and geometry of the simulation cell are chosen to accommodate a particular number of molecules ($`N`$) and hcp structure. Our calculations are performed mostly on $`N`$ =288 and from an initial simulation cell determined by two basis vectors (ap and bp) forming a 600 angle and by the third one (cp) perpendicular to them with the length ratio of $`a_p:b_p:c_p=1:1:4\sqrt{6}/9`$. The packing pattern is ABAB… to form the hcp lattice structure. Extensive testing on $`N`$ = 96, 150 and 392 has also been made.

In studying the effect of heating solid H2 isobarically at high pressures, we perform the path-integral Monte Carlo calculations with a constant-pressure ($`NPT`$) ensemble, instead of a simpler constant-volume ($`NVT`$) one, to avoid a bias of the restrictive cell geometry with a predetermined crystal structure. Implementation of the $`NPT`$ ensemble is achieved by extra independent Metropolis moves of three basis vectors ap, bp and cp, which generates a Markov chain of states having a limiting distribution proportional to $`\mathrm{exp}(H/T+N\mathrm{ln}V)`$ at $`T`$ with the enthalpy $`H`$($`=PV+E_s`$), where $`E_s`$ is the energy expectation value $`H_{\mathrm{nuclear}}`$ of the configuration $`s`$ with $`s`$ representing a set of scaled coordinates. This enables us to monitor volume changes and therefore to observe a possible first-order phase transition directly.

In order to avoid the “minus-sign problem” in QMC studies inherent in fermions, we confine ourselves to studying p-H2 and o-D2. After 4000 Monte Carlo steps for equilibration, statistical averages are collected from every second step, to a total of about 4000 data points. Because we focus on higher temperatures than the room temperature, relatively smaller partition number in the imaginary-time axis $`(M)`$ is enough. At a fixed $`P`$, we start with simulations at the room temperature at which $`M`$ is increased up to 16 to ensure the convergence. We come to know that $`M=8(5)`$ is enough for $`P=30`$GPa (higher $`P`$). Once its large enough value is determined, $`M`$ is held fixed during heating at the given $`P`$. This allows us to use the equilibrium configuration at low temperature as a start-up one for higher temperature, resulting in continuous heating of our sample. Because other computation techniques are documented in Ref. , we omit further details of them here.

## III RESULTS

We observe the equilibrium structure in real space directly and also monitor the pair distribution function $`g(R)`$ that represents the conditional probability of finding other molecules at the distance $`R`$ from a specified molecule at the origin:

$$g(R)=\frac{1}{4\pi R^2\rho }\underset{i<j}{}\delta (R_{ij}R),$$
(3)
with $`\rho `$ the density of system. This function shows a set of well-defined peaks characteristic to the configuration of neighboring molecules around the one at the origin. The function $`O(R)`$ defined in Ref.  is used to measure the correlation in molecular orientation.

The equilibrium distribution of 288 ($`=6\times 6\times 8`$) molecules in real space projected on both $`xy`$ and $`xz`$ planes are shown in Fig. 2 for p-H2 at 100GPa and at two different temperatures. Although it fluctuates with the standard deviation indicated by the painted region in Fig. 2(iii), on the average, each molecule occupies the lattice site of hcp at 450K, while that of an orthorhombic Cmca structure at 650K.

For detecting this transition quantitatively, we plot $`g(R)`$ at $`P=100`$GPa with the increase of $`T`$ in Fig. 3. Below 560K the second peak in $`g(R)`$ at $`R2.72`$Å, a characteristic feature of hcp, is clearly seen. But the peak disappears at higher $`T`$, which is a sign of the structural transition from hcp to Cmca. Thus we identify the transition temperature $`T_{tr}`$ as 560K with a statistical error of 20K. (By “statistical” we mean that the difference of states cannot be seen clearly if that in $`T`$ is less than 20K due to statistical fluctuations.) This structural change is always seen in the system with $`N=`$ 96, 150, or 392. We find a large $`N`$ dependence of $`T_{tr}`$, leading to an error of 40K or larger incurred at extrapolation of $`T_{tr}`$ at $`N\mathrm{}`$. This error cannot be reduced further, because calculations with a much larger $`N`$ are not feasible at present.

With the further increase of $`T`$ from $`T_{tr}`$, the peaks in $`g(R)`$ become less sharp. Eventually at $`T`$ above 875K (which we identify the melting temperature $`T_m`$ with a statistical error of 25K, about the same size of an error due to extrapolation in $`N\mathrm{}`$) $`g(R)`$ exhibits characteristics of a liquid phase. This first-order melting transition can be much better identified by the discontinuous change with $`T`$ in either the enthalpy $`H`$ or volume $`V`$ (the inset of Fig. 3).

By evaluating both $`T_{tr}`$ and $`T_m`$ at other pressures, we obtain the phase diagram of dense p-H2 in the $`TP`$ plane as shown in Fig. 4(a). The new solid phase is labeled as phase I’. Since $`O(R)`$ vanishes, its structure is orientationally disordered. The I’/I and liquid/I’ phase boundaries are indicated by solid triangles and open circles with error bars, respectively. Large curvature of the phase line for $`P<50`$GPa reflects the rapid change of the volume in that region. A Simon equation of $`P=1.2\times 10^5T_m^{2.355}`$ (the dotted curves in Fig. 4) fits our results rather well. A lack of experiment on $`T_m`$ at high pressures makes a precise comparison between theory and experiment difficult, but we obtain a reasonably good agreement with the data  for $`P<4`$GPa, as can be seen in Fig. 4(b). The data for n-H2 up to 15GPa,  are always higher than our results, reflecting the fact that the presence of o-H2 makes $`T_m`$ higher.  In Fig. 5(a) we show our result for dense o-D2 exhibiting a similar phase diagram. The difference in $`T_{tr}`$ between these two systems becomes smaller with the increase of $`P`$, but p-H2 gives definitely a higher $`T_{tr}`$ than o-D2. As for $`T_m`$, we find no meaningful difference.

## IV QUANTUM HARDENING EFFECT

For understanding the mechanism to bring about the phase I’, we have performed another simulation with $`M=1`$ (the “classical reference” system in which the molecules are treated as classical particles) and found an analogous phase diagram as shown in Fig. 5(b). In this system, only the temperature-induced random motion of molecules can play a role in transforming a compact hcp structure into a less compact Cmca one and finally into a liquid. Thus we can identify the thermal fluctuations as the basic driving force of the present phase transition.

Examining the changes of various physical quantities at $`T`$ around $`T_{tr}`$, we found that the most important change occurs in the potential energy. The energy increases with $`T`$ in both the classical (Fig. 6(a)) and the H2 (Fig. 6(b)) systems due primarily to the fluctuation of molecules to shorten the nearest-neighbor distance $`R_01.9`$Å, because $`V_{\mathrm{pair}}`$ is very repulsive at $`R`$ near $`R_0`$ and it changes very rapidly with $`R`$. (See Fig. 1.) We should note that the increasing rate is larger in the compact hcp phase than that in the less compact Cmca phase. In this regard, the phase transition is induced by the potential-energy gain as indicated in Fig. 6. No such gain is seen for the translational kinetic energy as shown in the inset of Fig. 6(b).

Eventually the reason why $`T_{tr}`$ in the classical system is much lower than that in solid p-H2 must be ascribed to the quantum zero-point motion of protons, the only ingredient that is not included in the classical system. Such a motion is, in a sense, regarded as a kind of vacuum polarization the effect of which, in the context of high-energy physics, can be renormalized into fundamental physical quantities such as mass and interaction constant. We envisage that this must also be the case here. In order to substantiate this, we consider a toy model, namely, an interacting two-particle system in a harmonic potential in one dimension, described by the Hamiltonian as

$$H=\frac{p_1^2}{2m}+\frac{Kx_1^2}{2}+\frac{p_2^2}{2m}+\frac{Kx_2^2}{2}+V(x_1x_2),$$
(4)
from which we can derive the Hamiltonian for the relative motion as

$$H_r=\frac{p^2}{2\mu }+\mu \omega ^2x^2+V(x),$$
(5)
where $`\mu =m/2,\omega =\sqrt{K/2\mu }`$ and $`x=x_1x_2`$. By writing $`\psi (x)=\psi _0(x)\phi (x)`$ and $`E=\omega /2+\epsilon `$ with $`\psi _0(x)=\sqrt[4]{\mu \omega /\pi }\mathrm{exp}(\mu \omega x^2/2)`$, we can cast Eq. (5) into

$$\left(\frac{p^2}{2\mu }+\stackrel{~}{V}(x)\right)\phi (x)=\epsilon \phi (x),$$
(6)
where the effective interaction $`\stackrel{~}{V}(x)`$ is related to the bare one through

$`\stackrel{~}{V}(x)V(x)\mu ^1{\displaystyle \frac{\psi _0^{^{}}}{\psi _0}}{\displaystyle \frac{\phi ^{^{}}}{\phi }}=V(x)+\omega x{\displaystyle \frac{\phi ^{^{}}}{\phi }},`$ (7)
We see that the effect of the zero-point oscillation is included effectively in $`\stackrel{~}{V}(x)`$, implying that the physics can be captured in terms of $`\stackrel{~}{V}(x)`$.  In this sense, $`V_{\mathrm{pair}}`$ can include the effect of intramolecular vibrons.

In general the closed-packed structure is known to become more stable with a harder $`V_{\mathrm{pair}}`$.  In the above toy model, we readily know that, if $`V(x)`$ is repulsive, $`x\phi ^{^{}}/\phi `$ is positive, leading $`\stackrel{~}{V}(x)`$ to be more repulsive than $`V(x)`$. This indicates that the zero-point motion of protons can make $`V_{\mathrm{pair}}`$ effectively more repulsive than that in the classical system in which $`\stackrel{~}{V}(x)`$ is reduced to $`V(x)`$ owing to $`\omega =\sqrt{K/m}0`$ as $`m\mathrm{}`$.

This hardening effect due to the molecular fluctuations can be understood more visually through Fig. 7. At the pressure range considered in this paper, the repulsive part of the interaction plays a central role in determining the potential energy. The fluctuation shortens one distance between the neighboring molecules, but at the same time it lengthens another one. Thus we have to average these two opposite effects. Because the potential changes very rapidly with $`R`$ the intermolecular distance, $`\mathrm{\Delta }V_1`$ the increase of the potential due to the shortening is larger than $`\mathrm{\Delta }V_2`$ the decrease due to the lengthening. Therefore the average potential energy, $`[V(R\delta R)+V(R+\delta R)]/2`$, becomes larger than $`V(R)`$ the potential energy without the fluctuation. As shown in the inset (a) of Fig. 7, the larger the fluctuation, the more repulsive the potential energy. This explains the behavior of the potential energy in Fig. 6; the solid H2 has a larger potential energy than the classical reference system at the same pressure and temperature, because it provides larger fluctuations due to the combination between the quantum and thermal effects. This also manifests itself in the much larger potential-energy gain in Fig. 6(b) than that in Fig. 6(a). Thus we can conclude that the molecules in the quantum system feel an effectively more repulsive potential than those in the classical system. This “quantum hardening effect” of $`V_{\mathrm{pair}}`$ explains why hcp becomes more stable in p-H2.

With the increase of pressure, the deviation of the I/I’ phase lines between solid H2 and the classical system is getting larger. In order to understand this result, we plot $`[V(R\delta R)+V(R+\delta R)]/2`$ as a function of $`R`$ with a given $`\delta R`$ in the inset (b) of Fig. 7. The average potential energy increases with the decrease of $`R`$, a behavior of $`R`$ seen with the increase of pressure. Thus we conclude that the quantum hardening effect becomes stronger with increasing pressure and this explains the above result well.

In order to study this quantum hardening effect more quantitatively, we have performed another simulation in the classical system with the pair potential, defined as

$$V_{\mathrm{pair}}=V_{\mathrm{SG}}+\eta V_{SR}+V_{ani},$$
(8)
where the parameter $`\eta `$ is introduced to measure the softness of $`V_{\mathrm{pair}}`$; $`\eta =0`$ for the hardest SG potential and $`\eta =1`$ for the softest SGH one. In Fig. 5(b), the open squares with error bars show the I/I’ phase boundary obtained using $`V_{\mathrm{pair}}`$ with $`\eta =0.61`$. Notice that this more repulsive potential than the SGH one makes the I/I’ phase boundary of the reference classical system almost the same as the one for the quantum solid H2 with the SGH potential (the double-dotted-dashed curve) at $`P=100`$GPa. We need a softer (harder) $`V_{\mathrm{pair}}`$ to better fit the data at lower (higher) pressure, but this behavior is consistent with the above-mentioned pressure dependence of the quantum hardening effect.

## V SUMMARY

Utilizing an intermolecular interaction potential $`V_{\mathrm{pair}}`$ proved to be reliable at high pressures in the PIMC simulations, we have found the structural phase transition from orientationally disordered hcp to $`Cmca`$ induced by thermal fluctuations in both dense solid H2 and D2 before melting. The result shows that a new phase exists at high temperatures in the solid molecular systems. The potential-energy gain due to the enhancement of fluctuations with increasing temperature is the main source to bring about the structural phase transition. The strong quantum zero-point motion plays an important role in determining $`T_{tr}`$ the transition temperature through the quantum hardening effect on the potential energy. The melting temperature $`T_m`$ is also estimated for full understanding the behavior of the solids and for the evaluating the temperature range of the new phase.

Through additional simulation studies, we have also examined to see how our conclusion about the existence of the new phase is robust against the choice of $`V_{\mathrm{pair}}`$. We have found that both $`T_{tr}`$ and $`T_m`$ are rather sensitive to the change of $`V_{\mathrm{pair}}`$. An example of $`T_{tr}`$ is given in Fig. 5(b). However, there is always a melting transition after the I/I’ transition with increasing temperature, suggesting the robustness of the existence of the $`Cmca`$ phase.

###### Acknowledgements.

This work is supported by Chinese National Foundation for Doctoral Education, Scientific Research Foundation for Returned Oversee Chinese Scholars, Chinese National Natural Science Foundation, and Japanese Monbusho Scholarship for Exchange Visiting Researchers.