text:
# Free energy calculations of elemental sulphur crystals via molecular dynamics simulations

## I Introduction

In recent papers we studied the crystalline phases of the molecular compound $`S_8`$, the most abundant natural form of elemental sulphur around ambient pressure and temperature (STP). Three crystalline allotropes of this compound are known: $`\alpha `$-$`S_8`$ phase is orthorhombic with 16 molecules in the non-primitive unit cell meyer1 ; steudel . If several samples of this crystal are slowly heated, some of them will show a phase transition to $`\beta `$-$`S_8`$ at 369 K, which in turn melts at 393 K steudel , but most of them show a metastable melting point at 385.5 K meyer1 ; steudel . Monoclinic $`\beta `$-$`S_8`$ is usually obtained from the melt steudel  and has 6 molecules per primitive unit cell, two of them orientationally disordered above 198 K beta1 . Monoclinic $`\gamma `$-$`S_8`$ has four molecules in a pseudo-hexagonal close packed unit cell, with a density 5.8% higher than $`\alpha `$-$`S_8`$ at STP meyer1 .

Using a very simple molecular model consisting of a cyclic semi-flexible chain of 8 atoms, with constant S-S bond lengths, it is possible to reproduce many features of the complex phase diagram of the $`S_8`$ molecule pap1 ; pap2 . The intra-molecular potential model includes a harmonic bending potential for S-S-S angles and a double well for the torsional angles cardini1 , the intermolecular potential used was a simple *Lennard-Jones* non-bonded atom-atom interaction, the details and potential parameters are given in the Section III. Using this simple interaction model we could reproduce, via a series of classical constant temperature- constant pressure simulations pap1 ; pap2 , the following experimental facts: the crystalline structure, the configurational energy and the lattice, bending and torsional dynamics (as given by the calculated density of vibrational modes) of $`\alpha `$-, $`\beta `$\- and $`\gamma `$-$`S_8`$ for T$``$200 K; the orientational dynamical order - disorder phase transition of $`\beta `$-$`S_8`$ and, finally, the solid-liquid phase transition of a cubic disordered sample was calculated near the experimental value pap1 ; pap2 .

Nevertheless, we found a fact that cannot be reproduced by this simple molecular model pap1 ; pap2 : when the temperature of a $`\alpha `$-$`S_8`$ MD sample of 288 molecules is lowered below 200 $`K`$, our orthorhombic $`\alpha `$-$`S_8`$ sample shows a structural phase transition to a monoclinic phase, with a molecular array similar to that of $`\alpha `$-$`S_8`$ and that we called $`\alpha ^{}`$-$`S_8`$ pap2 . This has not been experimentally observed. This spontaneous change was most probably due to the large fluctuations associated to a relatively small sample. A larger sample of 512 molecules didn’t show this distortion, but its configurational energy and volume per molecule remained, nevertheless, higher than those calculated for $`\alpha ^{}`$-$`S_8`$ pap2 . Fig. 1 shows both structures and the configurational energy vs. temperature, that of $`\alpha `$-$`S_8`$ was calculated by slowly decreasing T in steps of 25 to 50 K. The same study for $`\alpha ^{}`$-$`S_8`$ is performed by increasing T, after a long stabilization and annealing to obtain a totally ordered structure at the starting point of 50 K, this is the reason for the different $`U_{conf}`$ values, in fig. 1, at this point. Experimental measurements of the $`\alpha `$-$`S_8`$ structure at 300 and 100 K clearly disregard a structural phase transition, unless a metastable state has been experimentally measured alfa2 ; exp-alemanes , and this fact shows the limits of applicability of the simple molecular model.

In refs. pap2 ; pap3  we analyzed if an anisotropic non- bonded atom-atom intermolecular potential model was able to reproduce all the experimental facts. Using the program GAMESS gamess , we performed *ab initio* calculations of the electronic density distribution of the $`S_8`$ molecule, and the measured anisotropy of this distribution around each atom due to the atomic lone pairs. The calculated electronic density reproduces the experimental crystalline measurements of Coppens *et. al.* alfa3 , with a lone pair center at 0.7Å of the S location; nevertheless, at the nearest crystalline distances found between non-bonded atoms, the atomic anisotropy is low. The MD simulations performed with this anisotropic potential did not improve the previous result (shown in Fig. 1), unless a quite unrealistic atomic anisotropy is used pap3 ; pap4 .

It has to be pointed out that the search of a reliable molecular model for elemental sulphur molecules is extremelly useful due to the practical impossibility of performing quantum mechanical simulations of the complex phase diagram of these molecular crystals, with a large number of atoms in the primitive cells.

In this paper we decided to implement a reliable calculation of the free energy of $`\alpha `$-$`S_8`$ and $`\alpha ^{}`$-$`S_8`$ crystals in order to review the conclusions obtained with the first molecular model pap1 ; pap2 , based on an estimation of the free energy that takes into account the contribution of the vibrational modes to the zero point energy and entropy of the sample pap1 ; pap2 . Our present calculations are useful to check our previous free energy estimations and, at the same time, for the comparison of both free energy calculation methods. The usefulness of this comparison is due to the fast way in which the free energy can be estimated, within the quasi-harmonic approximation, using the data of a single MD simulation.

In the following sections we give the details of the implemented free energy calculations, the inter- and intra-molecular potential model used, the performed MD simulations, the obtained results and our conclusions.

## II Free energy calculations

In order to obtain reliable free energy values we implemented the thermodynamic integration method of Frenkel and Ladd free-frenkel-ladd1 , including the two correction terms given in refs. free-frenkel-correction ; understanding . From the results of a single MD simulation, we also estimated the free energy in the quasi-harmonic approximation, that is usually employed in lattice dynamics calculations bornhuang . In Section IV we give the details of our performed MD simulations.

### II.1 Thermodynamic integration method

This method proposes a thermodynamic integration along a reversible path between the system we are interested in and another one for which the free energy can be analytically calculated free-frenkel-ladd1 ; free-frenkel-correction ; understanding . For crystalline samples the reference system is an Einstein crystal, for which the internal energy is calculated as:

$$U_{Einstein}=U(\{𝐫_0^N\})+\underset{i=1}{\overset{N}{}}\alpha (𝐫_i𝐫_0)^2,$$

where $`\{𝐫_0^N\}`$ is the set of coordinates of the minimum energy configuration (the equilibrium atomic coordinates) and $`U(\{𝐫_0^N\})`$ its potential energy, the second term is the harmonic potential energy of the N oscillators. The force constant $`\alpha `$ is taken equal for all degrees of freedom and its value, $`\alpha =`$18.0844 kJ/mol, is chosen so as to obtain a mean square displacement equal to the average one in the real crystal: $`\mathrm{\Delta }r^2_{Einstein}\mathrm{\Delta }r^2_{real}`$. The potential energy difference between the real and the reference crystals can be calculated along a reversible artificial pathway that links both systems. At any intermediate point of the path, given by the coordinate $`\lambda `$ $`(0\lambda 1)`$, the local value of $`U`$ is

$$\stackrel{~}{U}(\lambda )=U(\{𝐫_0^N\})+(1\lambda )[U(\{𝐫^N\})U(\{𝐫_0^N\})]+\lambda \underset{i=1}{\overset{N}{}}\alpha (𝐫_i𝐫_0)^2.$$
(1)
So, for $`\lambda =0`$ the potential energy of the real crystal is recovered and, for $`\lambda \text{=1}`$ the potential energy is that of the Einstein crystal:

$$\stackrel{~}{U}(\lambda =0)=U_{sist}$$

and

$$\stackrel{~}{U}(\lambda =1)=U_{Einstein}.$$

The free energy difference between the real and the reference crystal is given by free-frenkel-ladd1 ; understanding :

$$F(\lambda =0)F(\lambda =1)=_0^1𝑑\lambda (\frac{F(\lambda )}{\lambda })=_0^1𝑑\lambda \frac{\stackrel{~}{U}(\lambda )}{\lambda }_\lambda ,$$

and deriving eq. 1:

$$F=F_{Einstein}+_{\lambda =1}^{\lambda =0}𝑑\lambda \underset{i=1}{\overset{N}{}}\alpha (𝐫_i𝐫_0)^2[U(𝐫^N)U(𝐫_0^N)]_\lambda .$$
(2)
The free energy of the Einstein crystal, with a fixed center of mass for the MD sample, is polson-frenkel1 :

$$\frac{F_{Einstein}}{N}=U(𝐫_0^N)\frac{3(N1)}{2N\beta }\mathrm{ln}(\frac{2\pi }{\alpha \beta })\frac{3\mathrm{ln}N}{2N}\frac{\mathrm{ln}V_0}{\beta N}+3\mathrm{ln}(\mathrm{\Lambda }),$$
(3)
where $`\mathrm{\Lambda }=\frac{h}{\sqrt{2\pi mk_BT}}`$ is the de Broglie thermal wavelength , $`\beta =\frac{1}{k_BT}`$, $`N`$ is the number of atoms and $`V_0`$ is the equilibrium volume of the sample.

If $`N_{mol}`$ is the number of molecules of the system, the final expression, including a correction term due to the center of mass constraint of the sample polson-frenkel1 ; free-frenkel-ladd1 ; understanding  and a second correction term due to its finite size free-frenkel-correction , is:

$$\frac{\beta F}{N_{mol}}=3\frac{N}{N_{mol}}\mathrm{ln}(\mathrm{\Lambda })\frac{3(N1)}{2N_{mol}}\mathrm{ln}(\frac{2\pi }{\beta \alpha })\frac{3\mathrm{ln}(N)}{2N_{mol}}\frac{\mathrm{ln}V_0}{N_{mol}}\frac{\beta }{N_{mol}}_0^1U_{ein}U_\lambda .$$
(4)
The accuracy of this calculation is determined by a correct evaluation of the last term in eq. 4, where long MD simulations are needed to measure the integrand for each value of $`\lambda `$.

For all values of $`\lambda `$, the center of mass of the sample is held fixed at the origin because the ‘springs’ of the reference Einstein crystal are tied to the equilibrium atomic positions. This is important to avoid divergences of the integrand in eq. 4 when $`\lambda `$ is close to 0 free-frenkel-ladd1 ; understanding .

The third term in eqs. 3 and 4 includes the change of the sign mentioned in ref. free-frenkel-correction , with respect to the original work of Frenkel and Ladd free-frenkel-ladd1 , due to a new consistent calculation of the partition function of the constrained CM system free-frenkel-correction . This term is negligible in a large system but must be taken into account to compute absolute free energies of a typical MD sample.

### II.2 Quasi-harmonic approximation

Following Born and Huang bornhuang , the free energy of a harmonic crystal with $`N_d`$ degrees of freedom ( i. e. $`N_d`$ independent harmonic oscillators) can be calculated as:

$$F=F_i=U_{conf}+\frac{1}{2}\underset{i=1}{\overset{N_d}{}}h\nu _i+k_BT\underset{i=1}{\overset{N_d}{}}\mathrm{ln}(1e^{\frac{h\nu _i}{k_BT}})\text{,}$$
(5)
where $`U_{conf}`$ is the potential energy and the $`\nu _i`$ are the vibrational frequencies of a system with $`N_d`$ degrees of freedom.

We have implemented this method, as in ref. pap1 ; pap2 , in two steps. First we calculate the atomic velocities self correlation function $`C(t)`$, using the MD data stored in a free trajectory in the phase space of the sample:

$$C(t)=\frac{v(t).v(0)}{v(0).v(0)}\text{,}$$

where $`<>`$ imply averages over atoms and over different initial times.

In a second step we calculate the vibrational density of states $`D(\nu )`$ *via* a Fourier transform of $`C(t)`$, and replace the sums over frequencies in eq. 5 by an integral on frequencies, weighted by the calculated density of vibrational modes. We have to take into account that the normalization factor for $`D(\nu )`$ is such that $`D(\nu )𝑑\nu =N_d.`$ In our case is $`N_d=16N_{mol}`$, because we have 8 bond constraints per molecule.

We can then calculate the internal energy, the entropy and free energy as bornhuang :

$$F=U_{conf}+\frac{1}{2}𝑑\nu D(\nu )h\nu +k_BT𝑑\nu D(\nu )\mathrm{ln}(1e^{\frac{h\nu }{k_BT}}),$$
(6)
$$S=(\frac{F}{T})_V=\frac{e^{\frac{h\nu }{k_BT}}}{1e^{\frac{h\nu }{k_BT}}}\frac{h\nu }{T}D(\nu )𝑑\nu ,$$

$$U=FTS=U_{conf}+\frac{1}{2}𝑑\nu D(\nu )h\nu +\frac{e^{\frac{h\nu }{k_BT}}}{1e^{\frac{h\nu }{k_BT}}}\frac{h\nu }{T}D(\nu )𝑑\nu .$$

The approximations involved in our quasi-harmonic calculations, with the results of a single MD simulation, are the following:

* $`U_{conf}`$ is calculated by performing averages, over time and molecules, of the potential energy per molecule calculated in the MD run, $`U_{conf}=<U>`$, it is not the value of the minimum potential energy $`U_0`$ of the Einstein crystal in equation 3.
* The calculated vibrational density of states $`D(\nu )`$ of the MD crystalline sample, due to its finite size ( a few number of primitive unit cells are included in it) does not give an accurate measurement of the density given by the real dispersion curves, the frequencies are only measured in a few points of the reciprocal space.
* The anharmonic frequencies of the sample can be accurately measured with a MD simulation, but eq. 6 is called the ’quasi-harmonic’ approximation because eq. 5 is exact only for harmonic potentials.

## III Intra- and intermolecular potential model

The potential model is that of refs. pap1 ; pap2 . The flexible molecular model includes all low frequency molecular distortions that mix with lattice modes and can therefore be of relevance in a possible structural phase transition. The S-S bond distances are kept constant ($`d_{SS}=2.0601`$ Å) because the stretching modes are well above in energies ($`\nu >400cm^1`$) than the rest of the vibrational modes ($`\nu <250cm^1`$).

We must stress here that bond constraints can have nonnegligible contributions to the free energy calculations. Such contributions have been studied elsewhere , with different methods and techniques constraints-estima1 ; constraint-estima2 ; constraints-estima3 . The main conclusions of those explicit calculations, using both MD constraint-estima2  and MCTI constraints-estima1 , are that the influence of bond constraints is actually negligible when the bond length is not changed from the initial to the final state. This is the case in our $`S_8`$ phases, in which the primitive cell will change and the molecules will distort from one phase to the other but the molecules bond length is held fixed within the same value for every calculation presented in this work. If changes in bond length between initial and final states of the free energy calculation were involved, an explicit constraint contribution to the free energy must be considered as explained in refs. constraints-estima1 ; constraints-estima3 . For the sake of completeness, a simple estimation of the change in the probability density between the standard $`NVT`$ and the ‘bond-constrained’ ensemble of $`S_8`$ is given in the Appendix.

The bending intramolecular potential for S-S-S angles is harmonic cardini1 ,

$$V(\beta )=\frac{1}{2}C_\beta (\beta \beta _0)^2$$

with a force constant of Cβ=25725 $`k_B/`$rad2 and $`\beta _0`$= 108 deg. The intramolecular potential for torsion angles is a double well cardini1 ,

$$V(\tau )=A_\tau +B_\tau \mathrm{cos}(\tau )+C_\tau \mathrm{cos}^2(\tau )+D_\tau \mathrm{cos}^3(\tau )\text{,}$$

with $`A_\tau `$=57.192 kB, $`B_\tau `$=738.415 kB, $`C_\tau `$=2297.880 kB and $`D_\tau `$=557.255 kB. These parameters describe a double well with minima at $`\tau `$=180 $`{}_{}{}^{+}{}_{}{}^{}`$ 98.8 deg., and a barrier height of about 9 kJ/mol at $`\tau `$=180 deg. The barrier height at $`\tau `$=0 deg. is of 30 kJ/mol, out of the range of energies explored in these simulations.

The intermolecular potential is of the non-bonded atom-atom Lennard-Jones type, with parameters $`\epsilon `$=1.70 kJ/mol and $`\sigma `$=3.39 Å cardini1 ; pap1 . The cut-off radius of the interactions is 15 Å and correction terms, to the energy and pressure, due to this finite cut-off are taken into account by integrating the contribution of an uniform distribution of atoms for distances larger than our cut-off of 15 Å.

## IV Molecular dynamics simulations

The MD simulations and samples are entirely identical to those previously performed in refs. pap1 ; pap2 . The bond constraints are held fixed with the SHAKE algorithm tildesley  and the temperature is maintained using the Nosé-Hoover chains method nose-hoover-chains , although the same behavior was found using the standard Nosé method nose-hoover1 ; nose-hoover2  in our NPT ensemble simulations.

Taking the equilibrated samples of $`\alpha `$\- and $`\alpha `$’-S$`_\text{8}`$ from our simulations in refs. pap1 ; pap2  as starting point, we first did a careful measurement of their crystalline structures (i. e. the lattice parameters and volume $`V_0`$) by averaging the calculated values over long runs performed in the NPT ensemble. The time step of the simulations was of 0.01 ps, the thermalization of the samples was of 80000 time steps (800 ps) and they were measured in the following 30000 (300 ps) time steps. The lattice parameters of each one of the four samples, together with their averaged atomic positions, were taken as the equilibrium atomic locations for the corresponding Einstein crystals.

The next series of MD simulations, as a function of the $`\lambda `$ parameter $`(0<\lambda <1)`$, were performed at constant volume, with a standard time step of 0.005 ps and runs of 160000 time steps (800 ps) of thermalization followed by a 40000 time steps (200 ps) of a free trajectory in the phase space, that were used to measure the systems.

For $`\lambda `$ values close to 1, we had to reduce the time step up to a minimum of 0.0001 ps in order to obtain the same accuracy in the total energy as for the other values of $`\lambda `$. Accordingly, the total number of time steps were increased so as to accumulate nearly the same total time of thermalization and storage of the other cases.

As in refs. polson-frenkel1 ; polson-frenkel2 , ten values of the parameter $`\lambda `$ are defined by a ten point Gauss-Legendre quadrature method, used to resolve the thermodynamic integration on the $`\lambda `$ coordinate (last term en eq. 4).

## V Results

### V.1 Thermodynamic integration method

Fig. 2 shows the values of $`U_{free}(\lambda )=U_{Einstein}U_\lambda `$, obtained with averages over the lasts 200 ps of each run, for $`\alpha `$\- and $`\alpha `$’-$`S_8`$ at 300 K. Similar curves $`U_{free}(\lambda )`$ are calculated for $`\alpha `$ and $`\alpha ^{}`$-$`S_8`$ at 100K, except that for $`\lambda `$1 $`<U_{free}>`$ is about $`260`$ kJ/mol. For both temperatures the values for the $`\alpha ^{}`$-$`S_8`$ sample are systematically lower than those of $`\alpha `$-$`S_8`$, the difference increases for increasing values of $`\lambda `$.

Note that the values $`U_{free}(\lambda )`$ near $`\lambda `$1 are those that weight more in the calculation of the integral in eq. 4. $`\lambda `$1 implies strong springs of Einstein crystals and weak interactions of the original system, in this case the slope of $`<U_{free}(\lambda 1)>`$ has a large increment because each atom is subjected almost only to the harmonic force of the Einstein crystal, allowing molecular configurations very different from those of the real system. Each atom shows a mean square displacement that depends only on temperature and not in the forces of the real system. At higher temperatures the interactions are more repulsive (because of the shorter mean interatomic distances) and $`<U_{free}>`$ takes more negative values ($`U`$ has a negative sign in the integrand of eq. 4). This can be verified noting that for $`T=100`$ K the minimum $`<U_{free}>`$ is approximately $`260`$ kJ/mol, almost half the value at $`300`$ K (-510 kJ/mol, see figure 2).

The contributions of the different terms to the total free energy of eq. 4, in adimensional units, are included in Table 1. The first three rows present the same values for all the cases, the first two of them are large terms that almost cancel each other. The next three rows show the terms that are calculated different in each case, and are obtained from the MD simulation data. The main contributions are from the $`U_0`$ static term and from the thermodynamic integration term. In the last row the adimensional free energy values for all the cases are shown.

The integral for thermodynamic integration in eq. 4 was also calculated using numerical integration for the Akima interpolation curve (see figure 2) and the difference with the Gauss-Legendre quadrature is within 6.4 %. The error of this term, included in table 1 was obtained with the same Gauss-Legendre quadrature using the values $`U_{free}(\lambda )+\mathrm{\Delta }U(\lambda )`$ and $`U_{free}(\lambda )\mathrm{\Delta }U(\lambda )`$, where $`\mathrm{\Delta }U`$ is the statistical error of each MD simulation.

There is a further comment about this free energy calculation method and its use in systems with multiple constraints, as are the very commonly used semiflexible molecular models with bond length constraints. In our simulations the constraints were not only applied to the real systems but also to the corresponding Einstein crystals. In this way we computed free energy differences between two constrained systems (the real and the reference) for each one of the four studied cases. This is irrelevant as long as we compare the free energy between the $`\alpha `$\- and $`\alpha `$’-$`S_8`$ crystals or between them and their reference systems. However, when computing the absolute free energy of each crystal, it has to be taken into account that the bond length constraints affect its MD trajectory in the phase space . The correction terms for systems with multiple types of constraints were recently discussed by Otter and Briels briels-fe-const , they found that the bond constraint contribution is not negligible only when comparing systems of different bond lengths constraint-estima2 ; constraints-estima1 ; constraints-estima3 .

### V.2 Quasi-harmonic approximation

The absolute values of the free energies of the four samples were also calculated *via* the contribution of the density of vibrational modes $`D(\nu )`$ to the entropy and zero point energy of the sample. With the measurement of an unique run in the $`NPT`$ ensemble, for each one of the four real systems, we obtained the values included in Table 2. The samples are equilibrated and measured at a pressure of $`0._{}^+0.03`$ kbar and the $`pV`$ term is not included to allow a direct comparison with the thermodynamic integration method. The volume per molecule $`(V_m)`$, configurational energy ($`U_{conf}`$), zero point energy ($`E_{Q_0}`$), entropy ($`S`$) and free energy ($`F`$) are shown. The free energy calculated with this approximation for various temperatures between 100 and 300 K, not included here, showed a systematically lower values for $`\alpha `$’-$`S_8`$.

### V.3 Comparison of both methods

Table 3 includes the calculated absolute values of free energy, in kJ/mol units, for $`\alpha `$ and $`\alpha ^{}S_8`$ at 100 and 300 K, using both methods.

We expect higher accuracy in the absolute free energy values for the thermodynamic integration method, mainly at higher temperatures when the harmonic approximation is worst than in the case of 100 K, due to the fact that anharmonic effects are expected to increase at higher temperatures.

With the thermodynamic integration method the difference of free energies is calculated about $`\mathrm{\Delta }F_{\alpha \alpha ^{}}(T=300K)2\frac{kJ}{mol}`$ , being $`\alpha `$’-$`S_8`$ the more stable phase. Although the absolute free energy values are not similar, the difference between both phases is also of about 2 kJ/mol for the quasi-harmonic method. This difference can be experimentally measured and closely follows that between the corresponding $`U_0`$ and $`U_{conf}`$ included in table 4.

## VI Conclusions

In this paper we calculated the free energy of orthorhombic $`\alpha `$-$`S_8`$ crystals and compared it with that of monoclinic $`\alpha `$’-$`S_8`$ pap1 ; pap2 ; pap3 ; pap4 , a previously proposed phase that is calculated with lower mean volume per molecule and lower potential energy than $`\alpha `$-$`S_8`$ crystals. The monoclinic $`\alpha `$’-$`S_8`$ crystal was obtained in the course of MD simulations, using simple inter- and intramolecular potential models, but has not been experimentally observed alfa2 ; exp-alemanes . Therefore, in this work we performed an accurate measurement of the free energy difference between both crystals at 100 and 300 K (above and below the spontaneous transition found in ref. pap1 ), as given by the simple molecular model of pap1 . The thermodynamic integration method free-frenkel-ladd1 ; free-frenkel-correction  is implemented, taking as reference system the corresponding Einstein crystals. For the sake of comparison we also included calculations of the free energy in the quasi-harmonic approximation, a fast estimation method that can be easily implemented in MD simulations.

The accurate free energy calculations performed here still show that monoclinic $`\alpha `$’-$`S_8`$ is calculated more stable than the orthorhombic experimentally observed $`\alpha `$-$`S_8`$, in accord with our previous calculations pap1 ; pap2  performed with the same simple molecular model and in the quasi-harmonic approximation.

As regards the accuracy of the free energy calculations, we found that for the thermodynamic integration method the results are extremely sensitive to the accurate calculation of the values $`U(\lambda )`$. In particular for $`\lambda 1`$ long MD simulations are needed. Although the absolute free energy values differ for both calculation methods, both measure a difference of about 2 kJ/mol between $`\alpha `$\- and $`\alpha `$’-$`S_8`$ crystals, being $`\alpha `$’-$`S_8`$ the more stable.

The approximated method for calculating the free energy *via* the contribution of the vibrational modes to the energy and entropy of the system turns out to be very useful to determine relative values between polymorphic phases or between different temperatures of a given sample. This result is very important, taking into account that this method is, at least, ten times faster than the thermodynamic integration method, that requires at least ten lengthly MD simulations (one MD run for each value of the $`\lambda `$ parameter), in order to calculate the free energy of one sample at a given temperature.

###### Acknowledgements.

C. P. thanks Ignacio Urrutia for fruitful and encouraging discussions as regards bond constrains in free energy calculations and FOSDIC for a partial support. This work was partially supported by the grant PIP 0859/98 of CONICET and Fundación Balseiro-2001.

## Appendix A APPENDIX

We estimate here the difference in the calculated statistical averages for MD simulations performed with and without molecular bond length constraints. This topic has been extensively treated in literature and is not of minor importance fe-const2 ; fe-const3 ; fe-const4 ; understanding . Indeed, it has been widely studied in the case of free energy calculations that involve a reaction coordinate fe-const3 ; fe-const4 .

We recall the results of chapter 10 from ref. understanding  and the appendix from ref. fe-const2 . The expression for the ratio of molecular position probability for constrained and unconstrained system is:

$$\frac{\rho (𝐪_1,\mathrm{},𝐪_{3N})}{\rho (𝐪_1,\mathrm{},𝐪_{3Nl,}𝐪_1^c=\sigma _1,\mathrm{},𝐪_l^c=\sigma _l)}=c\sqrt{\left|Z\right|},$$
(7)
where we explicitly note the constrained degrees of freedom with $`𝐪_i^c=\sigma _i`$, being the standard bond length constraint,

$$\sigma _i=\sqrt{m_i}(\sqrt{𝐫_{ij}^2}d0),$$
(8)
and $`l`$ constrained bonds are assumed. The $`l\times l`$ $`Z`$ matrix, in the right hand side of 7 is defined as understanding ; fe-const2 :

$$Z_{mn}=\underset{i=1}{\overset{N}{}}\frac{1}{m_i}(\frac{\sigma _m}{𝐫_i}).(\frac{\sigma _n}{𝐫_i}),$$
(9)
and $`\left|\right|`$ stands for the determinant. The origin of $`Z`$ is from the fact that not only $`\sigma (𝐫_i)=0`$ must be required to constraint the bonds in the simulations, but also $`\dot{\sigma }(𝐫_i)=0`$ at all times fe-const4 . It should be noted that $`𝐪`$ in eq. 7 is used for generalized coordinates while $`𝐫_i`$ for cartesian ones.

For the $`S_8`$ molecule we have 24 degrees of freedom, from which 8 have been constrained in our model (see sec. III). Therefore we have 16 degrees of freedom with 8 bond constraints following eq. 8, with $`d=`$2.06 Å. Matrix $`Z`$ can be explicitly calculated for a $`S_8`$ molecule from eq. 9:

$$\left(\begin{array}{cccccccc}2& \mathrm{cos}(\beta _{12})& 0& 0& 0& 0& 0& \mathrm{cos}(\beta _{18})\\ \mathrm{cos}(\beta _{12})& 2& \mathrm{cos}(\beta _{23})& 0& 0& 0& 0& 0\\ 0& \mathrm{cos}(\beta _{23})& 2& \mathrm{cos}(\beta _{34})& 0& 0& 0& 0\\ 0& 0& \mathrm{cos}(\beta _{34})& 2& \mathrm{cos}(\beta _{45})& 0& 0& 0\\ 0& 0& 0& \mathrm{cos}(\beta _{45})& 2& \mathrm{cos}(\beta _{56})& 0& 0\\ 0& 0& 0& 0& \mathrm{cos}(\beta _{56})& 2& \mathrm{cos}(\beta _{67})& 0\\ 0& 0& 0& 0& 0& \mathrm{cos}(\beta _{67})& 2& \mathrm{cos}(\beta _{78})\\ \mathrm{cos}(\beta _{18})& 0& 0& 0& 0& 0& \mathrm{cos}(\beta _{78})& 2\end{array}\right),$$

where $`\beta _{ij}`$ is the bending bond angle between consecutive angles $`i,j`$ and $`j+1`$.

If we take a mean $`\beta `$ value for all the bending angles and evaluate the determinant in eq. 7, we get:

$$\left|Z\right|=256576\mathrm{cos}(\beta )^2+400\mathrm{cos}(\beta )^488\mathrm{cos}(\beta )^6.$$
(10)
We can now take a statistical mean value over the NVT ensemble, averaging eq. 10 over a MD run and normalizing with the maximum value. This provides an estimation of the effect of constraints over statistical averages such as used in eq. 4, to carry out the thermodynamic integration. We got the following main difference between the constrained and unconstrained probabilities at $`T=300K`$:

$$\frac{\rho }{\rho _{constrained}}\sqrt{\left|Z\right|}0.954.$$

The closeness to 1 indicates that the effect of bond constrains is practically negligible for simulations in the studied range of thermodynamic variables. From this estimation we can conclude that the points obtained in figure 2 are almost identical to those that can be obtained in a similar simulation without bond constrains, that would have been quite more expensive in CPU time. This simple calculation, that we present here for close ring-molecules, is complementary to that already obtained in ref. const-fixman  for open chain molecules.

### Contribution of stretching modes to free energy

We include here an estimation of the free energy difference between samples with and without molecular bond constrains, under the framework of quasi-harmonic approximation commented in sec. II.2. If we think in the normal modes as harmonic oscillators of frequency $`\nu _i`$, we have a ‘trivial’ partition function $`Z_i`$ for each oscillator of the formbornhuang :

$$Z_i=e^{i\frac{1}{2}h\nu _i}\underset{s=0}{\overset{\mathrm{}}{}}e^{sh\nu _i/kT}=\frac{e^{i\frac{1}{2}h\nu _i}}{1e^{ih\nu _i}}$$

where $`h`$ is the Planck constant, $`k`$ the Boltzmann factor and $`s`$ the positive integer quantum number of the oscillator of frequency $`\nu _i`$. For a system of such independent oscillators we have a partition function:

$$Q_{simul}=_0^{\mathrm{}}D(\nu )(\frac{e^{\frac{1}{2}\beta h\nu }}{1e^{\beta h\nu }})𝑑\nu ,$$
(11)
where $`D(\nu )`$ is the density of vibrational states already defined in sec. II.2and $`Q_{simul}`$ is calculated from our bond constrained MD simulation. The partition function of the unconstrained system $`Q`$ can be estimated in the following way: taking into account that for the $`S_8`$ molecule the stretching vibrational modes are well above bending and torsional modes in the $`\nu `$ scale (see sec. III), we have separable contributions for the total partition function of the molecule and

$$Q=_0^{\nu _{max}}D_{simul}(\nu )(\frac{e^{\frac{1}{2}\beta h\nu }}{1e^{\beta h\nu }})𝑑\nu +_{_{\nu _{max}}}^{\mathrm{}}D_{stretch}(\nu )(\frac{e^{\frac{1}{2}\beta h\nu }}{1e^{\beta h\nu }})𝑑\nu Q_{simul}+Q_{stretch},$$
(12)
where we have assumed that the interval $`(0,\nu _{max})`$ accounts for the bending, torsional (dihedral) and lattice vibration modes and the high frequencies of bond stretching are ‘isolated’ in the interval $`(\nu _{max},\mathrm{})`$. The first term in eq. 12 is obtained from our simulation and the second one can be simply approximated by:

$$Q_{stretch}=N^{stretch}(\frac{e^{\frac{1}{2}\beta h\nu _{stretch}}}{1e^{\beta h\nu _{stretch}}}),$$

here $`N^{stretch}`$ are the 8 stretching degrees of freedom (constrained in the MD simulations) and $`\nu _{stretch}`$ is the mean stretching frequency. This is equivalent to use the relationship $`D_{stretch}(\nu )=N^{stretch}\times \delta (\nu \nu _{stretch})`$. The free energy difference between considering or not the stretching degrees of freedom can be estimated, under the quasi-harmonic approximation by:

$$\beta \mathrm{\Delta }F=\beta (F_{tot}F_{simul})=\mathrm{ln}(\frac{Q_{simul}}{Q_{simul}+Q_{stretch}})$$

For $`T=300K`$ and the sulphur mean stretching frequency $`\nu _{stretch}=453.8cm^1`$, we get $`\beta \mathrm{\Delta }F=0.075`$ which is a relatively and also an absolute small number when compared with the $`\beta F`$ values of table 1.

This simple calculation can be considered complementary to those of ref. constraints-estima1 ; constraints-estima3 ; constraint-estima2 , that provided a thorough discussion of the role of bond constrains in free energy calculations.


format: text

text:
# Time-reversal violating Schiff moment of 225Ra

## I Introduction

Experiments with K and B mesons indicate that time-reversal invariance (T) is violated through phases in the Cabibbo-Kobayashi-Maskawa matrix that affect weak interactions Kowalewski . The suspicion that extra-Standard-Model physics, e.g. supersymmetry, also violates T has motivated a different kind of experiment: measuring the electric dipole moments (EDMs) of the neutron and of atoms. Because any such dipole moment must be proportional to the expectation value of the T-odd spin operator, it can only exist when T (and parity) is violated Sachs (1987); Khriplovich and Lamoreaux (1997). So far the experiments have seen no dipole moments, but they continue to improve and even null results are useful, since they seriously constrain new physics. Whatever the experimental situation in the future, therefore, it is important to determine theoretically what the presence or absence of EDMs at a given level implies about T-violating interactions at elementary-particle scales. Our focus here is atoms, which for some sources of T violation currently provide limits as good or better than the neutron Romalis et al. (2001).

One way an atom can develop an EDM is through T and P violation in its nucleus. Let us assume that given a fundamental source of the broken symmetry one can use effective-field theory and QCD to calculate the strength of the resulting T-violating nucleon-pion interaction. One then needs to connect the strength of that interaction to the resulting nuclear “Schiff moment”, which, because the nuclear EDM is screened Schiff (1963), is the quantity responsible for inducing an EDM in electrons orbiting the nucleus. The Schiff moment is defined classically as a kind of radially weighted dipole moment:

$$𝑺=\frac{1}{10}\mathrm{d}^3r\rho _{\mathrm{ch}}(𝒓)\left(r^2\frac{5}{3}\overline{r_{\mathrm{ch}}^2}\right)𝒓,$$
(1)
where $`\rho _{\mathrm{ch}}`$ is the nuclear charge density and $`\overline{r_{\mathrm{ch}}^2}`$ is the mean-square charge radius. Recent papers Spevak et al. (1997); Auerbach et al. (1996) have argued that because of their asymmetric shapes, octupole-deformed nuclei in the light-actinide region should have collective Schiff moments that are 100 to 1000 times larger than the Schiff moment in 199Hg, the system with the best experimental limit on its atomic EDM Romalis et al. (2001). Ref. Engel et al. (1999) suggested that certain many-body effects may make the enhancement a bit less than that. The degree of enhancement is important because several experiments in the light actinides are contemplated, planned, or underway Chupp ; Holt . They may see nonzero EDMs, and even if they don’t we will need to be able to compare their limits on fundamental sources of T violation to those of experiments in other isotopes.

Perhaps the most attractive octupole-deformed nucleus for an experiment is 225Ra. Though radioactive, it has a ground-state angular momentum $`J=1/2`$, which minimizes the effect of stray quadrupole electric fields in an experiment to measure a dipole moment111The statement that the nucleus has octupole and quadrupole deformation really refers to its intrinsic state, a concept we elaborate on below, and does not contradict its insensitivity to applied electric fields with multipolarity greater than one.. In addition, the associated atom has close-lying electronic levels of opposite parity and is relatively easy to trap and manipulate. As a result, at least one group is at work on a measurement in 225Ra Holt . Here we calculate its Schiff moment, attempting to incorporate the effects discussed in Ref. Engel et al. (1999) through a symmetry-unrestricted mean-field calculation. We begin in the next section by describing the physics of the Schiff moment in octupole-deformed nuclei, briefly reviewing prior work in the process. In Section III we test our mean-field approach by calculating properties of even Ra isotopes. In Section IV we discuss issues peculiar to mean-field calculations in odd nuclei and then present our results for the Schiff moment of 225Ra, focusing particularly on the degree of enhancement. Section V is a brief conclusion.

## II Enhancement of Schiff Moments in Octupole-Deformed Nuclei – Previous Work

In analogy with dipole moments in atoms, static Schiff moments in nuclei can exist only if T is broken. Because T-violating forces are much weaker than the strong interaction, the Schiff moment can be accurately evaluated through first-order perturbation theory as

$$S\mathrm{\Psi }_0|\widehat{S}_z|\mathrm{\Psi }_0=\underset{i0}{}\frac{\mathrm{\Psi }_0|\widehat{S}_z|\mathrm{\Psi }_i\mathrm{\Psi }_i|\widehat{V}_{PT}|\mathrm{\Psi }_0}{E_0E_i}+\text{c.c.},$$
(2)
where $`|\mathrm{\Psi }_0`$ is the member of the ground-state multiplet with $`J_z=J0`$, the sum is over excited states, and $`\widehat{S}_z`$ is the operator

$$\widehat{S}_z=\frac{e}{10}\underset{p}{}\left(r_p^2\frac{5}{3}\overline{r_{\mathrm{ch}}^2}\right)z_p,$$
(3)
with the sum here over protons. The operator $`\widehat{V}_{PT}`$ is the T- (and parity-) violating nucleon-nucleon interaction mediated by the pion Haxton and Henley (1983); Herczeg (1988a) (shown to be more important than other mesons in Ref. Towner and Hayes (1994)):

$`\widehat{V}_{PT}(𝒓_1𝒓_2)`$ $`={\displaystyle \frac{gm_\pi ^2}{8\pi m_N}}`$ $`\{(𝝈_1𝝈_2)(𝒓_1𝒓_2)[\overline{g}_0\stackrel{}{\tau }_1\stackrel{}{\tau }_2{\displaystyle \frac{\overline{g}_1}{2}}(\tau _{1z}+\tau _{2z})+\overline{g}_2(3\tau _{1z}\tau _{2z}\stackrel{}{\tau }_1\stackrel{}{\tau }_2)]`$ (4)
$`{\displaystyle \frac{\overline{g}_1}{2}}(𝝈_1+𝝈_2)(𝒓_1𝒓_2)(\tau _{1z}\tau _{2z})\left\}{\displaystyle \frac{\mathrm{exp}(m_\pi |𝒓_1𝒓_2|)}{m_\pi |𝒓_1𝒓_2|^2}}\right[1+{\displaystyle \frac{1}{m_\pi |𝒓_1𝒓_2|}}],`$
where arrows denote isovector operators, $`\tau _z`$ is +1 for neutrons, $`m_N`$ is the nucleon mass, and we are using the convention $`\mathrm{}=c=1`$. The $`\overline{g}`$’s are the unknown isoscalar, isovector, and isotensor T-violating pion-nucleon couplings, and $`g`$ is the usual strong $`\pi NN`$ coupling.

In a nucleus such as 199Hg, with no intrinsic octupole deformation, many intermediate states contribute to the sum in Eq. (2). By contrast, the asymmetric shape of 225Ra implies the existence of a very low-energy $`|1/2^{}`$ state, in this case 55 keV above the ground state $`|\mathrm{\Psi }_0|1/2^+`$, that dominates the sum because of the corresponding small denominator. To very good approximation, then,

$$S\frac{1/2^+|\widehat{S}_z|1/2^{}1/2^{}|\widehat{V}_{PT}|1/2^+}{\mathrm{\Delta }E}+\text{c.c.},$$
(5)
where $`\mathrm{\Delta }E`$ = 55 keV. The small denominator is part of the reason for the enhancement of the Schiff moment. The other part is the matrix element of the Schiff operator in Eq. (5). In the limit that the deformation is rigid, the ground state and first excited state in octupole-deformed nuclei are partners in a parity doublet, i.e., projections onto good parity and angular momentum of the same “intrinsic state” that represents the wave function of the nucleus in its own body-fixed frame. The matrix elements in Eq. (5) are then proportional (again, in the limit of rigid deformation) to intrinsic-state expectation values, so that Spevak et al. (1997)

$$S2\frac{J}{J+1}\frac{\widehat{S}_z\widehat{V}_{PT}}{\mathrm{\Delta }E},$$
(6)
where $`J`$ is the ground-state angular momentum, equal to 1/2 for 225Ra, and the brackets indicate expectation values in the intrinsic state. The intrinsic-state expectation value $`\widehat{S}_z`$ is generated by the collective quadrupole and octupole deformation of the entire nucleus; it is much larger than a typical matrix element in a spherical or symmetrically deformed nucleus. Together with the small energy denominator, this large matrix element is responsible for the enhancement of laboratory-frame Schiff moments in nuclei such as 225Ra.

The amount of the enhancement is not easy to calculate accurately, however. The reason is that the matrix element of the two-body spin-dependent operator $`\widehat{V}_{PT}`$ in Eq. (5) depends sensitively on the behavior of a few valence particles, which carry most of the spin. In the approximation that particles (or quasiparticles) move in independent orbits generated by a mean field, the potential can be written as an effective density-dependent one-body operator that we will denote $`\widehat{U}_{PT}`$, defined implicitly by

$$a|\widehat{U}_{PT}|b=\underset{c<F}{}ac|\widehat{V}_{PT}|bc,$$
(7)
where $`|a`$, $`|b`$, and $`|c`$ are eigenstates of the mean field and the matrix elements of $`\widehat{V}_{PT}`$ are antisymmetrized. With the further approximation that the mass of the pion is very large, $`\widehat{U}_{PT}`$ can be written as a local operator, in a form we display in the Section IV. Evaluating its matrix element is tricky.

The authors of Refs. Spevak et al. (1997); Auerbach et al. (1996) used a version of the particle-rotor model Leander and Sheline (1984) to represent the odd-$`A`$ nucleus. In this model, all but one of the nucleons are treated as a rigid core, and the last valence nucleon occupies a deformed single-particle orbit, obtained by solving a Schrödinger equation for a Nilsson or deformed Wood-Saxon potential. The model implies that the core carries no intrinsic spin whatever, that the neutron and proton densities are proportional, and that the exchange terms on the right-hand side of Eq. (7) are negligible. Under these assumptions, $`\widehat{U}_{PT}`$, which now acts only on the single valence nucleon, reduces to Sushkov et al. (1984)

$$\widehat{U}_{PT}(𝒓)\eta \frac{G}{2m_N\sqrt{2}}𝝈\mathbf{}\rho _0(𝒓).$$
(8)
where $`G`$ is the Fermi constant, inserted to follow convention, and $`\rho _0`$ is the total nuclear mass density. The dimensionless parameter $`\eta `$ is then a function of the couplings $`\overline{g}_i`$ and the isospin of the nucleus.

Ref. Engel et al. (1999) confirmed the collectivity of the intrinsic Schiff moments obtained in Refs. Spevak et al. (1997); Auerbach et al. (1996), but questioned the accuracy of some of the assumptions used to evaluate the matrix element of $`\widehat{V}_{PT}`$, suggesting that either core-spin polarization or self-consistency in the nuclear wave function might reduce laboratory Schiff moments. The zero-range approximation and the neglect of exchange in $`\widehat{U}_{PT}`$ are also open to question. As a result, it is not clear whether the Schiff moment of 225Ra is 1000 times that of 199Hg or 100 times, or even less. In what follows, we provide a (tentative) answer by moving beyond the particle-rotor model. Our calculation is not the final word on Schiff moments in octupole-deformed nuclei — we only do mean-field theory, neglecting in particular to project onto states with good parity, and do not fully account for the pion’s nonzero range — but is a major step forward.

## III Mean-field calculations for other Ra isotopes

### III.1 Mean-field calculations

Self-consistent mean-field theory is widely used for describing bulk properties of nuclei Bender et al. (2003). In the guise of density-functional theory, it is also used throughout atomic and molecular physics. The approach is more “microscopic” — nucleons are the only degrees of freedom — and far less phenomenological than the collective particle-rotor model. Self-consistency connects the single-particle states and the actual density distribution. The variational principle that determines the single-particle wave functions thus optimizes all multipole moments not fixed by global symmetries. The density distributions of neutrons and protons are not proportional to each other; they have slightly different deformations and radial profiles. In odd-$`A`$ nuclei, self-consistent calculations include rearrangement due to the unpaired particle. Rearrangement causes polarization of the even-$`A`$ core through orbital-current and spin-density terms in the effective interaction. Core polarization is one of the effects on the Schiff moment of 225Ra that we investigate below.

Our approach is nonrelativistic and employs Skyrme interactions. To get an idea of the range of results this kind of calculation can produce, we use four different parameterizations of the Skyrme energy functional, i.e., four different Skyrme forces. The four give similar results for many observables near stability, but still have significant differences. Our favorite interaction, for reasons explained below, is SkO’ Bender et al. (2002); Reinhard et al. (1999), but we also show results for the commonly used forces SIII Beiner et al. (1975), SkM Bartel et al. (1982), and SLy4 Chabanat et al. (1998).

### III.2 Related Observables in Even Isotopes

Intrinsic-parity breaking in even radium isotopes is the subject of several theoretical analyses; see the review in Ref. Butler and Nazarewicz (1996) and the more recent studies in Refs. Garrote et al. (1997); Tsvetkov et al. (2002). To assess the ability of the Skyrme interactions to handle it, we perform a series of Hartree-Fock (HF) + BCS calculations for even Radium isotopes. We use the Skyrme-HF+BCS code from Ref. Bender (1998); it represents single-particle wave functions on an axially symmetric mesh, and uses Fourier definitions of the derivative, $`1/r_{}`$, and $`1/r_{}^2`$ operators. We choose 75 grid points in the $`z`$ direction, and 27 in the $`r_{}`$ (perpendicular) direction, with 0.8 fm between them. The code uses a density-independent zero-range pairing interaction with a self-adjusting cutoff as described in Ref. Bender et al. (2000). For each Skyrme force we adjust the pairing strength separately for protons and neutrons Bender et al. (2000). We should note that other self-consistent mean-field models, namely HF+BCS with the nonrelativistic Gogny force Egido and Robledo (1989) and the relativistic mean-field model Rutz et al. (1995), yield results that are similar to those we describe now.

Figure 1 illustrates the calculated evolution of intrinsic deformation with increasing neutron number in the Radium isotopes. It plots the intrinsic ground-state mass-density contours predicted by SkO’. The mean-field ground states go from having a spherical shape at the magic number $`N=126`$ to a quadrupole deformed (reflection-symmetric) shape at $`N=130`$, then to quadrupole+octupole deformed (reflection-asymmetric) shapes for $`132N140`$, and finally back to quadrupole deformed shapes at higher $`N`$. Because the ground states are obtained from a variational principle, all shape moments higher than octupole are also optimized (the isoscalar dipole moment is constrained to be zero). The nucleus 225Ra, with $`N=137`$, will clearly be well deformed in both the quadrupole and octupole coordinates. The structures at small radii visible for $`N132`$ reflect small oscillations of the density distribution around the saturation value (for a given neutron excess) caused by shell effects.

We must note that the octupole-deformed minima are not equally pronounced for all forces and isotopes. In addition, in some of the isotopes with reflection-symmetric minima, some of the Skyrme forces predict an excited octupole-deformed minimum separated by a small barrier from the ground-state minimum. Furthermore, in the transitional nuclei, which have soft potentials in the octupole direction, all parity-breaking intrinsic deformations are subject to collective correlations as discussed in Ref. Egido and Robledo (1989). The influence of correlations will be smallest for the nuclides with the most pronounced octupole-deformed minima, usually 222Ra and 224Ra. This fact supports our belief that our mean-field calculations supply a good approximation to the intrinsic structure in 225Ra.

Figure 2 shows the relative error in the predicted binding energies $`\delta E=(E_{\text{calc}}E_{\text{expt}})/E_{\text{expt}}`$ for all four forces, and the predicted two-neutron separation energies, along with the measured values. All the forces do a good job with binding, which is not surprising given the way their parameters were fit. The fact that the error in binding for SkO’ is nearly constant with $`N`$ for $`N>130`$ is reflected in the near perfect agreement in the bottom panel with the measured two-neutron separation energies $`S_{2n}`$. The errors in predicted values of $`S_{2n}`$ around $`N=128`$ probably reflect the deficiencies of mean-field models in transitional nuclei.

Figure 3 shows three parity-violating intrinsic quantities. In the top panel is the ground-state octupole deformation $`\beta _3=4\pi r^3Y_{30}/(3AR^3)`$ (where $`R=1.2A^{1/3}`$), as a function of neutron number. The trend mirrors that in the density profiles shown earlier. At $`N=136`$, one less neutron than in 225Ra, the forces all predict almost identical octupole deformation, a result we like. Experimental data for octupole moments are still sparse in this region; we are only aware of $`\beta _3=0.105(4)`$ for $`N=138`$, a value that can be deduced from the $`B(E3;0_1^+3_1^{})`$ given in Ref. Spear (1989). (In 224Ra and 226Ra, by the way, we agree fairly well with the quadrupole moments obtained from $`B(E2)`$’s in Ref. Raman et al. (2001). For example, SkO’ gives $`\beta _2=0.184`$ in 224Ra and experiment gives $`\beta _2=0.179(4)`$.)

The second panel in the figure shows the absolute values of intrinsic dipole moments $`D_0=e_pz_p`$, along with experimental data extracted from $`E1`$ transition probabilities Butler and Nazarewicz (1996). The calculated values for $`D_0`$ change sign from positive to negative between $`N=134`$ and $`N=138`$, reflecting a small change in the location of the center of charge from the “top” half of the pear-shaped nucleus to the “bottom” half. This predicted sign change is consistent with the near-zero experimental value for $`N=136`$. None of the forces precisely reproduces the trend through all the isotopes, but the comparison has to be taken with a grain of salt because “data” derive from transitions between excited rotational states, and therefore are not necessarily identical to the ground-state dipole moments. Cranked Skyrme-HF calculations without pairing correlations Tsvetkov et al. (2002) and cranked HFB calculations with the Gogny force Garrote et al. (1997) predict that for most Ra isotopes $`D_0`$ changes significantly with angular momentum. In any event, as thoroughly discussed in Ref. Butler and Nazarewicz (1996), the intrinsic dipole moment is a small and delicate quantity.

The intrinsic Schiff moment $`S_z`$, the quantity we’re really interested in, is more collective and under better control, as the bottom panel of the figure shows. The various predictions are usually within 20$`\%`$ of one another and large, confirming the predictions originally made in Refs. Spevak et al. (1997); Auerbach et al. (1996). The octupole deformation and intrinsic dipole moment have been shown to change only slightly with parity projection from the intrinsic states Garrote et al. (1997), and the same is probably true of the intrinsic Schiff moment.

By turning the pairing force off, we are able to see whether the parity-violating quantities in Fig. 3 are affected by pairing correlations. In 224Ra, for example, SkO’ gives $`\beta _3=0.141`$, $`D_0=0.103`$ $`e`$ fm, and $`S_z=34.4`$ $`e`$ fm3 without pairing, and $`\beta _3=0.143`$, $`D_0=0.093`$ $`e`$ fm, and $`S_z=34.3`$ $`e`$ fm3 when pairing is included. In this nucleus uncertainties related to pairing are very small.

Finally, in Fig. 4, we show the predicted proton and neutron single-particle spectra generated by the ground-state mean-field in 224Ra. The combination of quadrupole, octupole, and higher deformations reduces the level density around the Fermi surface for both kinds of nucleon, leading to significant deformed $`Z=88`$ and $`N=132`$ shell closures for all interactions, and a somewhat weaker $`N=136`$ subshell closure for SIII, SkM\* and SkO’. The small level density around the Fermi surface might explain the insensitivity of the deformation to pairing correlations mentioned above. For all the forces except SkM, the first empty neutron level clearly has $`j_z=1/2`$, implying that in 225Ra the ground-state parity-doublet bands will be built on $`J^\pi =1/2^\pm `$ states. For SkM the situation is less clear because the $`j_z=1/2`$ and 3/2 states are nearly degenerate, and it is necessary to carry out the calculation in 225Ra itself to see which becomes the lowest.

## IV Calculating the Schiff moment of 225Ra

### IV.1 Odd-$`A`$ Nuclei and Schiff Moments in Mean-Field Approximation

Fully self-consistent calculations in odd-$`A`$ nuclei are possible but seldom performed. For many physical observables it is enough to neglect correlations between the odd particle and the core, which amounts to dropping a valence particle into the field generated by that core. Other quantities, however, are sensitive to the interaction between the last particle and the core. The interaction can change the deformation and pairing strength, and produce various kinds of core polarization. In fact, a self-consistent odd-$`A`$ calculation is equivalent to first performing the calculation in the even-even nucleus with one less neutron, then placing the last nucleon in the first empty orbit and treating its polarizing effect on the core in RPA Brown (1970b); Blaizot and Ripka (1986). Fully self-consistent calculations take all kinds of polarization into account simultaneously. For us, the spin polarization is of particular interest.

To proceed here we use the public-domain Skyrme-Hartree-Fock code HFODD (v2.04b) Dobaczewski and Dudek (1997, 2000); Olbratowski et al. , which allows the mean-field and the associated Slater determinant to simultaneously break invariance under rotation, reflection through three planes, and time reversal (the code we used in the last section does not). Breaking the first is necessary to describe a deformed nucleus. Breaking all three reflections (and not only parity) is necessary to represent axial octupole deformation with the spin aligned along the symmetry axis. Breaking the last induces spin polarization in the core, which because of Kramers degeneracy cannot otherwise occur. Incorporating spin polarization is important because it has the potential to significantly alter the matrix element of $`\widehat{U}_{PT}`$ in Eq. (8) from its value in the particle-rotor model, where the spin is carried entirely by one valence particle. The code HFODD cannot yet treat pairing when it allows T to be broken, but pairing in T-odd channels is poorly understood. No existing codes can do more than HFODD in odd-$`A`$ octupole-deformed nuclei.

As above, we use the Skyrme interactions SIII, SkM, SLy4, and SkO’. The reason SkO’ is our favorite has to do with the part of the energy functional composed of T-odd spin densities (which, following common practice, we refer to as the “T-odd functional”, even though the entire functional must be even under T). The T-odd functional plays no role in the mean-field ground states of even nuclei, but can be important in any state with nonzero angular momentum. Of the forces above, only SkO’ has been seriously investigated in T-odd channels. In Ref. Bender et al. (2002), the T-odd part of the functional was adjusted to reproduce Gamow-Teller resonances, resulting in an effective Landau parameter $`g_0^{}=1.2`$. In the isoscalar channel, the force was adjusted to reproduce the commonly used value $`g_0=0.4`$ Osterfeld (1992). Although there are not enough data to constrain other relevant parameters in the functional, and although a very recent calculation starting from a realistic nucleon-nucleon interaction Zuo et al. (2003), while confirming the value $`g_0^{}=1.2`$, finds $`g_0=0.85`$, SkO’ is clearly the best available Skyrme interaction for describing spin-spin interactions. The corresponding T-odd terms in the functional are precisely those that will polarize the spin in the core. There are other terms that could be added to the standard Skyrme interaction and do the same thing —- the tensor force for example — but they are almost never used and their effects still need to be investigated.

The parameters of the other three forces SIII, SkM, and SLy4 were adjusted entirely to ground-state properties in even-even nuclei, and so the Landau parameters $`g_0`$ and $`g_0^{}`$ were not fit. Here we set them to zero by treating the T-odd and T-even terms in the Skyrme functional independently, as described in Ref. Bender et al. (2002). Only orbital terms, which are fixed by gauge invariance Dobaczewski and Dudek (1995) (a generalization of Galilean invariance), appear in the T-odd parts of these forces.

We rely on the crudest forms of projection. For parity, that means none at all, and for angular momentum it means inserting the rigid-rotor factor $`J/(J+1)=1/3`$ in front of the intrinsic Schiff moment, as described above. In other words, we use Eq. (6), with the intrinsic state taken to be the Hartree-Fock ground state produced by HFODD. Just as in the particle-rotor model, the intrinsic Schiff moment is given by the classical expression Eq. (1), but with $`\rho _{\mathrm{ch}}`$ equal to $`e`$ times the Hartree-Fock ground-state proton density. As already mentioned, the Hartree-Fock approximation allows us, by summing over occupied orbits, to write the intrinsic matrix element of the two-body potential $`\widehat{V}_{PT}`$ as the expectation value of an effective one-body operator $`\widehat{U}_{PT}`$. Because we now have a microscopic version of the “core”, this effective potential is more complicated than in Eq. (8), and it now acts on all the nucleons:

$$\widehat{U}_{PT}=\frac{g}{2m_\pi ^2m_N}\underset{i=1}{\overset{A}{}}𝝈_i\tau _{z,i}\mathbf{}\mathrm{d}^3r^{}\left(\frac{m_\pi ^2e^{m_\pi |𝒓𝒓^{}|}}{4\pi |𝒓𝒓^{}|}\right)\left[(\overline{g}_0+2\overline{g}_2)\rho _1(𝒓^{})\overline{g}_1\rho _0(𝒓^{})\right]+\text{exch.}$$
(9)
Here $`\rho _0(𝒓)\rho _n(𝒓)+\rho _p(𝒓)`$ and $`\rho _1(𝒓)\rho _n(𝒓)\rho _p(𝒓)`$ are the isoscalar and isovector densities. The piece coming from exchange terms in the original two-body interaction $`\widehat{V}_{PT}`$ is nonlocal, just as in the usual Hartree-Fock mean field, and we have not written it explicitly here (though we do below).

The code HFODD at present cannot evaluate the expectation value of a folded potential like that above, which is due to the finite pion range. Nevertheless, even in the zero-range approximation we can avoid several of the assumptions — proportionality of neutron and proton densities, negligibility of exchange terms, and absence of core spin — leading to the extremely simplified potential in Eq. (8). The zero-range approximation is equivalent to assuming the pion is very heavy, so that the term involving the pion mass in Eq. (9) becomes a delta function. Under that assumption, but none others, the exchange terms become local and $`\widehat{U}_{PT}`$ takes the form:

$`\widehat{U}_{PT}(𝒓)`$ $`{\displaystyle \frac{g}{2m_\pi ^2m_N}}`$ $`\{{\displaystyle \underset{i=1}{\overset{A}{}}}𝝈_i\tau _{z,i}[(\overline{g}_0+2\overline{g}_2)\mathbf{}\rho _1(𝒓)\overline{g}_1\mathbf{}\rho _0(𝒓)]`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{A}{}}}𝝈_i[(3\overline{g}_0+\overline{g}_1\tau _{z,i})𝑱_0(𝒓)+(\overline{g}_1+\overline{g}_0\tau _{z,i}4\overline{g}_2\tau _{z,i})𝑱_1(𝒓)]\}.`$
Here $`𝑱(𝒓)`$ is the “spin-orbit” current, defined, e.g., in Ref. Bender et al. (2002) and references therein, and the subscripts 0 and 1 refer to isoscalar and isovector combinations as they do for the density. The terms in $`\widehat{U}_{PT}`$ that contain $`𝑱`$ are the exchange terms omitted above. We will evaluate them, but argue later that their effects are probably small when the finite range is restored. The terms containing the density $`\rho `$ all come from the direct part of $`\widehat{V}_{PT}`$. We do not simplify things further to obtain something like Eq. (8) because $`\rho _p`$ is not really proportional to $`\rho _n`$ and the core nucleons do carry some spin. We will manage nevertheless, to compare our results with those of Ref. Spevak et al. (1997). We will also estimate the effect of a finite pion range on the direct terms, though our inability to do so more precisely at present is the most significant shortcoming of this work.

HFODD works by diagonalizing the interaction in the eigenbasis of an optimal anisotropic three-dimensional harmonic oscillator. For 225Ra, algorithms developed in Ref. Dobaczewski and Dudek (1997) give oscillator frequencies of $`\mathrm{}\omega _z`$=7.0625 and $`\mathrm{}\omega _{}`$=8.6765 MeV in the directions parallel and perpendicular to the elongation axis. The matrix element of $`\widehat{U}_{PT}`$ converges only slowly as we increase the number of levels in the basis. When the interaction polarizes the core, it takes 2500 or more single-particle basis states to get convergence. The basis then contains up to $`N_z`$=26 and $`N_{}`$=21 oscillator quanta.

### IV.2 Laboratory Schiff Moment of 225Ra

We turn finally to results in 225Ra itself. For SkO’, our HFODD calculations yield $`\beta _2=0.190`$, $`\beta _3=0.146`$, and $`\beta _4=0.136`$ for the usual first order approximation to the deformation parameters determined from mass multipole moments Ring and Schuck (1980). The laboratory Schiff moment, Eq. (6), is proportional to the product of the intrinsic Schiff moment $`\widehat{S}_z`$ and the expectation value $`\widehat{V}_{PT}`$. Table 1 shows the intrinsic moments and the expectation values of the 6 operators that enter the zero-range approximation to $`\widehat{V}_{PT}`$ in Eq. (IV.1). Before commenting on the entries, we mention what is in each of the forces and calculations.

For all the forces, terms in the functional that are proportional to Laplacians of spin densities ($`𝒔\mathrm{\Delta }𝒔`$) and density-dependent spin-spin terms ($`f(\rho )𝒔𝒔`$), cf. Ref. Dobaczewski and Dudek (1995); Bender et al. (2002), which enter through the T-odd part of the Skyrme functional, have been turned off. For the first three lines in Table 1 \[forces labeled with (0)\], the spin-spin terms have also been turned off, so that the only nonzero terms in the T-odd functional (as noted above) are those required by gauge invariance. For the fourth line \[SkO’(—)\], all T-odd terms in the functional have been turned off, so that aside from the self-consistency in the wave functions the calculation resembles one with a phenomenological (non-self-consistent) potential, for which T-odd mean-fields are never considered. We include this result so that we can distinguish the role played by core polarization. The results in the line labeled SkO’ include the time-odd channels, adjusted as discussed above Bender et al. (2002). This is the force in which we have the most confidence. The last entry is the result of Ref. Spevak et al. (1997), with the implicit assumption that the neutron and proton densities are proportional.

In our calculations, the intrinsic Schiff moments are close to one another, and all are less than twice the estimate of Ref. Spevak et al. (1997). The agreement reflects the collective nature of these intrinsic moments; they are even larger than the particle-rotor estimates. But the matrix elements of $`\widehat{V}_{PT}`$, the other ingredient in Eq. (6) for the laboratory Schiff moment, are a bit more delicate. Our results show the exchange terms on the right side of the table to be comparable to the direct terms, a result that is surprising because for a spin-saturated core (or in the particle-rotor model) the exchange terms vanish exactly. We think, however, that the ratio of exchange to direct terms would become small were the finite range of the interaction reintroduced and short-range NN correlations inserted.

Though unable to include either effect here, we did so in a Nilsson model for 225Ra. We took nucleons there to occupy independent single-particle levels generated by a deformed oscillator potential with $`\beta _2=0.138`$, $`\beta _3=0.104`$, and $`\beta _4=0.078`$, values taken from Ref. Spevak et al. (1997). We then evaluated the ground-state expectation value of the full two-body interaction $`\widehat{V}_{PT}`$, with and without the zero-range approximation (and in the latter case, with short-range correlations included à la Ref. Miller and Spencer (1976)). In this simple model, the valence nucleon carries all the spin, and only the neutron-proton and neutron-neutron parts of $`\widehat{V}_{PT}`$ contribute. The direct $`np`$ term shrank by a factor of 1.5, while the corresponding exchange term shrank by a factor of 1400 (both independently of the $`\overline{g}`$’s in Eq. (4), it turns out) when the range of the interaction was set to its proper value. The results in the $`nn`$ channel were less dramatic: the direct part again shrank by 1.5 and the exchange part by a factor of 5. When we moved the valence neutron to higher orbits, these numbers changed some — the direct terms sometimes were not suppressed at all and other times shrank by factors of up to 6, but the ratios of the exchange to direct contributions almost always ended up small. Similar behavior was found for parity-violating forces in Ref. Adelberger and Haxton (1985), where it was traced in part to the different average momenta carried by the pion in direct and exchange graphs. So that we can compare our results with those of Ref. Spevak et al. (1997), we will neglect the exchange terms from now on, though we caution that this step should eventually be justified more rigorously, e.g., by actually calculating them with the finite-range force in the full mean-field model. The reduction we see in the direct terms is in line with the results of Ref. Griffiths and Vogel (1991), though we find it more variable222We performed another test, using the direct part of Eq. (IV.1) with the valence wave function taken from the Nilsson model just described, but with the neutron and protons densities assumed to have more realistic Wood-Saxon forms. The direct terms were again suppressed by factors of 1.5 to almost 10 that depended significantly on the valence orbit..

Though we cannot yet be more quantitative about finite-range effects, we do quantify the core polarization in Table 1. For the first three lines of the table, where the forces are labeled (0), the spin-spin terms are absent from the energy functional, and the protons in the core develop only a tiny spin density from the T-odd terms required by gauge invariance. For the fourth line, SkO’(—), all T-odd terms are absent and the protons can have no spin at all. This means that the operators $`f(𝒓)𝝈`$ and $`f(𝒓)𝝈\tau `$ have either the same or almost the same expectation value for any $`f(𝒓)`$ so that columns 4 and 6 ($`𝝈𝑱_0`$ and $`𝝈\tau 𝑱_0`$) have identical or nearly identical entries for these forces, and so do columns 5 and 7 ($`𝝈𝑱_1`$ and $`𝝈\tau 𝑱_1`$). The fifth line of the table contains the effects of spin polarization, which are primarily to alter the neutron-spin density; the equalities between the columns are not badly broken, so the protons do not develop much spin. The same is true of the terms involving $`\rho `$, though that is not obvious from the table because we display only the two terms that appear in Eq. (IV.1).

These near equalities and the probable irrelevance of the exchange terms when the finite range is taken into account imply that only the quantities $`𝝈_n\mathbf{}\rho _n`$ and $`𝝈_n\mathbf{}\rho _p`$ are ultimately important. We display them in Table 2. Except for SIII, the neutron-density distribution affects the matrix element much more than the that of protons. By comparing the fourth and fifth lines, however, we see that spin correlations increase the role of the protons, while reducing that of the neutrons slightly. Thus, while the spin-spin interactions do not cause the protons to develop much net spin, they do correlate the neutron spin with the proton density.

There is not too much scatter in our results. The entries in the second column ($`𝝈\tau \mathbf{}\rho _0`$) of Table 1 differ by factors of two or three, and the entries in the third ($`𝝈\tau \mathbf{}\rho _1`$) by a little more, though they are all smaller than those in the second column (which is not surprising — the third column subtracts the neutron and proton densities while the second adds them). In the neutron-proton scheme (table 2) all of our numbers are smaller than those of Ref. Spevak et al. (1997), a result that was anticipated in Ref. Engel et al. (1999). The difference from the earlier estimate for the larger matrix elements ranges from factors of two to four, though the isovector combination — the third column in table 1 — is sometimes actually enhanced a little.

What, at last, have we to say about the real laboratory Schiff moment $`S`$? The lab moment is given by the product of the matrix elements just discussed, the intrinsic Schiff moments, and the unknown coefficients $`\overline{g}_i`$. Our intrinsic Schiff moments are about 1.5 times larger than those of Ref. Spevak et al. (1997), while our $`\widehat{V}_{PT}`$ matrix elements, in the zero-range approximation, are smaller than theirs, usually by a somewhat larger amount. Overall, our lab moments will usually be smaller by factors between about 1.5 and 3 than the estimates of Ref. Spevak et al. (1997) (an exception can occur if for some reason $`\overline{g}_1`$ is considerably less than the other two coefficients).

How big are our moments compared to that of 199Hg? The most comprehensive calculation in that nucleus, which appeared very recently Dmitriev and Sen’kov , improved on the work of Ref. Flambaum et al. (1986) by including the effects of the residual strong interaction and the full finite-range form for $`\widehat{V}_{PT}`$. The new results are smaller than that of ref. Flambaum et al. (1986), only slightly so for the isovector part of $`\widehat{V}_{PT}`$, but by a considerably amount in the isoscalar and isotensor channels. The authors write their results in terms of the pion-nucleon couplings as

$$S_{\mathrm{Hg}}=.0004g\overline{g}_0+.055g\overline{g}_1+.009g\overline{g}_2[e\mathrm{fm}^3].$$
(11)
Our result for radium, with the zero-range approximation and exchange terms neglected, translates to

$$S_{\mathrm{Ra}}^{\mathrm{zero}\mathrm{range}}=5.06g\overline{g}_0+10.4g\overline{g}_110.1g\overline{g}_2[e\mathrm{fm}^3].$$
(12)
If the three $`\overline{g}`$’s are comparable, therefore, our Schiff moment is several hundred times larger than that of Ref. Dmitriev and Sen’kov , in part because the isoscalar and isotensor interactions are more effective in Ra than in Hg. \[If $`\overline{g}_1`$ is larger than the other two couplings, as in left-right symmetric models Herczeg (1988b), our result is less than 200 times bigger than the latest one in 199Hg. The very small coefficient of $`\overline{g}_0`$ for 199Hg in Eq. (11), by the way, has significant consequences Herczeg (1988b) for the limit on the QCD T-violating parameter $`\overline{\theta }`$ that can be inferred from the experimental limit in Ref. Romalis et al. (2001).\] Accepting the work of Ref. Dzuba et al. (2002) on atomic physics in Ra and Hg, the enhancement of the atomic EDM of 225Ra is about about three times that of the Schiff moment, i.e. potentially more than 1000. We again caution, however, that we have yet to include the full finite-range version of $`\widehat{V}_{PT}`$ and that our preliminary investigations suggest that doing so will decrease our Schiff moment at least a little. Ironically, Ref. Dmitriev and Sen’kov  finds that including the finite range actually increases the matrix element in 199Hg, though only slightly.

We hope to make other improvements in our calculation as well. Projection onto states of good parity will change the results a bit, as will angular-momentum projection. Our conclusions about the size of spin-polarization corrections could be modified by two terms in the Skyrme functional we’ve set to zero, or by a better determined value of the Landau parameter $`g_0`$. Whatever the result of such corrections, however, it is clear that the atomic EDM of 225Ra will always be significanlty larger than that of 199Hg.

## V Conclusions

We have calulated the Schiff moment in 225Ra in the approximation that the T-violating interaction has zero range. Our calculations, which are self-consistent and include core polarization, give results that are generally just a few times smaller than earlier estimates based on the particle-rotor model. Accepting the very recent results of Ref. Dmitriev and Sen’kov , we currently find the Schiff moment of 225Ra to be (generically) several hundred times that of 199Hg, a result that strengthens the case for an atomic-EDM experiment in Ra, though the enhancement factor depends significantly on the source of T violation, and we expect it to decrease at least a little when we use the finite-range force. Work towards including a finite range in HFODD is in progress. We also plan to apply the self-consistent methods used here to other light actinides, as well as to 199Hg, where we suspect octupole correlations may play some role Engel et al. (1999). Maintaining self consistency in 199Hg should automatically control the spurious Schiff strength encountered in Ref. Dmitriev and Sen’kov . The source of the insensitivity of the Schiff moment to T violation in the isoscalar channel in that work should be checked and understood.

After many years of neglect, the question of which isotopes are best for EDM measurements is now being rapidly addressed.

###### Acknowledgements.

This work was supported in part by the U.S. Department of Energy under grant DE-FG02-02ER41216, by the Polish Committee for Scientific Research (KBN) under Contract No. 5 P03B 014 21, and by computational grants from the *Regionales Hochschulrechenzentrum Kaiserslautern* (RHRK) Germany and from the Interdisciplinary Centre for Mathematical and Computational Modeling (ICM) of the Warsaw University. M. B. acknowledges support through a European Community Marie Curie Fellowship, and J. H. J. acknowledges partial support through a Doctoral Fellowship from the Portuguese Foundation for Science and Technology.