text:
# The Functional Form of Angular Forces around Transition Metal Ions in Biomolecules

## I Introduction

Atomic-level simulations of proteins interacting with transition-metal ions have the promise of elucidating a wide variety of biophysical phenomena (Kaim, 1994). Such simulations can help establish both the native structure and reaction paths of metalloproteins, and explain the roles of metal ions in protein folding and stability. Similar applications can be seen for simulations of metal ions interacting with DNA and RNA. However, the utility of such simulations depends critically on the availability of accurate but computationally tractable force fields for metals interacting with proteins. There are no quantitatively accurate force fields for these purposes, and it is probably not possible to construct force fields of a simple enough form to be useful that have quantitative accuracy. However, one can hope to generate force fields that include the basic physical effects at a level of accuracy sufficient to ascertain chemical trends as metals or ligands vary. The situation in this regard is best for “simple” metal ions that have no partly occupied shells. Here an ionic picture based on electrostatics, supplemented by empirical repulsive forces, may hope to treat the most important physical effects. Even here though, care should be taken, since the bond is not completely ionic but has some partial covalent character. For transition metals, the situation is much more difficult. The ligand-field splitting of the $`d`$-shell leads to important electronic contributions which cannot be ignored. These are manifested, for example, in the typical square or tetragonal coordination of Cu2+ complexes. Such effects cannot be described by radial interactions, but rather require the introduction of angular forces describing the energetics of the transition-metal $`d`$-shell. Most existing methods for including angular forces in simulations of metal ions have used assumed functional forms (Comba et al., 1995, Allured et al., 1991, Timofeeva et al., 1995, Wiesemann et al., 1994, Sayle et al., 1995, Sayle et al., 1997) based, for example on the observed structures of small complexes. A recently developed method (Landis et al., 1998) treats the environment-dependent energetics of transition metals via a valence-bond approach. This method appears promising for covalently bonded metals. However, the bonding state of transition metals in proteins is may be rather different from the $`sd^n`$ hybridized configuration assumed in that work. In the present analysis, we seek to develop a method suitable for more ionic configurations.

We have recently shown (Carlsson, 1998) how the ligand-field splitting effects of $`\sigma `$-bonded transition metal ions can be described by a force field of a fairly simple form. One starts with an explicitly quantum-mechanical form for the electronic energy, and then then solves the quantum-mechanical problem with systematic approximation methods to extract the real-space description of the electronic bonding energy. A second-order treatment of the hybridization terms between the ligand and metal orbitals is used to generate a ligand-field Hamiltonian for the $`d`$-electrons. The electronic bonding energy of this Hamiltonian is analyzed in terms its moments. The moment analysis of the energy gives a “semiclassical” energy function with a simple trigonometric angular form. It is not precisely an additive function of ligand-ligand interactions, but is nearly as simple computationally. This energy function gives quite accurate energy results for $`\sigma `$-bonded ligands, with less than 10% error in the ligand-field stabilization energy (defined below).

This approach is extended in the present paper, in two ways. First, we develop energy functions which treat $`\pi `$-bonded ligands as well. This results in new and distinct angular forms. Second, we present results for the lowest-energy structures of small clusters, which are a guide to understanding the bonding preferences of transition metals in proteins.

The organization of the remainder of the paper is as follows. Section II develops the formalism underlying the angular forces. Section III gives tests of the functional form of the angular forces by comparing with results from diagonalization of simple cluster Hamiltonians. Section IV gives the small-cluster results for ideal examples, and for an analogue of the copper environment in blue-copper proteins.

## II Model

We model the electronic structure of the environment of a single transition metal ion in terms of orbitals $`|L,\mu `$ localized on the ligands and $`|M,\nu `$ localized on the single metal ion. The orbitals taken to be orthogonal for simplicity of calculation. Here “L” denotes a particular ligand atom and the $`|L,\mu `$ are distinct orbitals on that atom. The Hamiltonian takes the following form:

$`\widehat{H}=`$ $`{\displaystyle \underset{L,\mu }{}}\epsilon _{L,\mu }|L,\mu L,\mu |+{\displaystyle \underset{\nu }{}}\epsilon _{M,\nu }|M,\nu M,\nu |`$ (2)
$`+{\displaystyle \underset{L,\mu ,\nu }{}}\left[h_{LM}^{\mu \nu }|L,\mu M,\nu |+h_{ML}^{\nu \mu }|M,\nu L,\mu |\right].`$
Here the $`\epsilon `$ terms are on-site energies for the orbitals, and the $`h`$ terms are hybridization energies between the ligand orbitals and the metal orbitals. We ignore electron-interaction terms as well as explicit electrostatic effects. The energy associated with this Hamiltonian can then be obtained by diagonalizing its matrix as given in the $`|L,\mu |M,\nu `$ basis, and taking the sum of the energies of the occupied eigenvectors.

To obtain a real-space description of the energetics of this Hamiltonian, we make the simplifying approximation that the dimensionless ratios $`|h_{LM}^{\mu \nu }|/(\epsilon _{M,\nu }\epsilon _{L,\mu })`$ are small and can be used as expansion parameters. This will be the case if the bonding is primarily ionic as opposed to covalent. In the case of transition metal ions, this allows a great simplification by the application of ligand-field theory, which gives an effective Hamiltonian for the $`d`$-shell:

$$\widehat{H}_d=\underset{\mu ,\nu }{}h_d^{\mu \nu }|M,\mu M,\nu |,$$
(3)
where

$$h_d^{\mu \nu }=\underset{L,\eta }{}(\epsilon _d\epsilon _{L,\eta })^1h_{ML}^{\mu \eta }h_{LM}^{\eta \nu }.$$
(4)
Note that we have assumed all of the $`d`$-orbitals on the transition metal ion to have the same energies before ligand-field effects are “turned on”.

In a fully quantum-mechanical ligand-field-theory calculation, one would numerically diagonalize the matrix $`\widehat{H}_d`$. In order to obtain a force field with a nearly classical form, we instead work with the moments of this matrix, defined in terms the traces of its powers, as follows. The first moment is defined as the average energy of the $`d`$-complex, or

$`\overline{\epsilon }`$ $`=`$ $`(1/5)\mathrm{Tr}\widehat{H}_d=(1/5){\displaystyle \underset{L,\eta }{}}(\epsilon _d\epsilon _{L,\eta })^1{\displaystyle \underset{\mu }{}}h_{ML}^{\mu \eta }h_{LM}^{\eta \mu }`$ (5)
$`=`$ $`(1/5){\displaystyle \underset{L,\eta }{}}(\epsilon _d\epsilon _{L,\eta })^1g_{\eta \eta }`$ (6)
where

$$g_{\eta ,\eta ^{}}=\underset{\mu }{}h_{L^{}M}^{\eta ^{}\mu }h_{ML}^{\mu \eta }.$$
(7)
Thus $`\overline{\epsilon }`$ is given as a sum of independent contributions from the ligands. Because the $`d`$-shell by itself has spherical symmetry, the contribution from each ligand is a radial function (no angular dependence) of the metal-ligand distance. Our major interest is in the angular terms resulting from the ligand-field splitting, so we do not consider the $`\overline{\epsilon }`$ term further.

The width $`W`$ of the $`d`$-complex corresponding to the ligand-field splitting can in the simplest picture be described in terms of the second moment or variance $`\delta \epsilon ^2`$ of the eigenvalues of $`\widehat{H}_d`$. We expect that $`W\sqrt{\delta \epsilon ^2}`$. We can obtain a simple real-space form for $`\delta \epsilon ^2`$ as follows: First we note that $`\delta \epsilon ^2=(1/5)_n(\epsilon _n\overline{\epsilon })^2=(1/5)\mathrm{Tr}(\widehat{H}_d\overline{\epsilon }\widehat{I})^2`$ where the $`\epsilon _n`$ are the eigenvalues of $`\widehat{H}_d`$. From Eqs. (3), (4), and (7), we see that

$`\mathrm{Tr}\widehat{H}_d^2`$ $`=`$ $`{\displaystyle \underset{\mu ,\nu }{}}h_d^{\mu \nu }h_d^{\nu \mu }`$ (8)
$`=`$ $`{\displaystyle \underset{L,\eta ,L^{},\eta ^{}}{}}(\epsilon _d\epsilon _{L,\eta })^1(\epsilon _d\epsilon _{L^{},\eta ^{}})^1{\displaystyle \underset{\mu }{}}h_{L^{}M}^{\eta ^{}\mu }h_{ML}^{\mu \eta }{\displaystyle \underset{\nu }{}}h_{LM}^{\eta \nu }h_{ML^{}}^{\nu \eta ^{}}`$ (9)
$`=`$ $`{\displaystyle \underset{L,\eta ,L^{},\eta ^{}}{}}(\epsilon _d\epsilon _{L,\eta })^1(\epsilon _d\epsilon _{L^{},\eta ^{}})^1g_{\eta \eta ^{}}^2.`$ (10)
Then the variance is

$`\delta \epsilon ^2`$ $`=`$ $`(1/5)\mathrm{Tr}\left(\widehat{H}_d\overline{\epsilon }\widehat{I}\right)^2=(1/5)\left(\mathrm{Tr}\widehat{H}_d^25\overline{\epsilon }^2\right)`$ (11)
$`=`$ $`(1/5){\displaystyle \underset{L,\eta ,L^{},\eta ^{}}{}}(\epsilon _d\epsilon _{L,\eta })^1(\epsilon _d\epsilon _{L^{},\eta ^{}})^1\left[g_{\eta \eta ^{}}^2(1/5)g_{\eta \eta }g_{\eta ^{},\eta ^{}}\right].`$ (12)
We will focus on the ligand-field stabilization energy, which is defined as the sum of the energies of the occupied orbitals relative to $`\overline{\epsilon }`$, or $`E_{\mathrm{LFSE}}=_n^{}(\epsilon _n\overline{\epsilon })`$, where the sum is over only the occupied eigenfunctions. We assume that $`E_{\mathrm{LFSE}}`$ is proportional to the ligand-field splitting $`W`$, and also identify $`W`$ with the standard deviation $`\sqrt{\delta \epsilon ^2}`$. Thus our expression for the ligand-field stabilization energy has the form

$$E_{\mathrm{LFSE}}=\left\{\alpha /5\underset{L,\eta ,L^{},\eta ^{}}{}(\epsilon _d\epsilon _{L,\eta })^1(\epsilon _d\epsilon _{L^{},\eta ^{}})^1\left[g_{\eta \eta ^{}}^2(1/5)g_{\eta \eta }g_{\eta ^{},\eta ^{}}\right]\right\}^{1/2},$$
(13)
where $`\alpha `$ is a dimensionless constant which will later be used as a fitting parameter.

The accuracy of this form will be demonstrated in the next section via specific numerical experiments. At present, we will show how this form for $`E_{\mathrm{LFSE}}`$ results in an expression for the energy in terms of simple angular interactions between he ligands. We thus need to develop analytic forms for the $`g_{\eta \eta ^{}}`$. First we note that we can write

$`g_{\eta \eta ^{}}`$ $`=`$ $`L,\eta |\widehat{H}P_d^2\widehat{H}|L^{},\eta ^{}`$ (14)
$`=`$ $`\mathrm{\Psi }|\mathrm{\Psi }^{},`$ (15)
where $`\widehat{P}_d=_\mu |M,\mu M,\mu |`$ is the projection onto the $`d`$-subspace of the transition metal ion Here we define $`|\mathrm{\Psi }=\widehat{P}_dH|L,\eta `$ and we have have used the relation $`\widehat{P}_d^2=\widehat{P}_d`$, which follows from $`\widehat{P}_d`$’s being a projection operator.

Consider first the case where both the $`|L,\eta `$ and $`|L^{},\eta ^{}`$ orbitals have $`\sigma `$-character with respect to the metal ion. Then they couple only to the $`d`$-orbitals $`|M,\sigma `$ and $`|M,\sigma ^{}`$ that have $`\sigma `$-character with respect to the bond axes, so that $`|\mathrm{\Psi }=h_\sigma |M,\sigma `$ and $`|\mathrm{\Psi }^{}=h_\sigma ^{}|M,\sigma ^{}`$, where $`h_\sigma `$ and $`h_\sigma ^{}`$ are the appropriate coupling strengths. Thus $`g_{\eta \eta ^{}}=h_\sigma h_\sigma ^{}M,\sigma |M,\sigma ^{}`$, and the angular dependence is contained in the last inner product. But this is simply a matrix element of a rotation about $`M`$ which carries $`L`$ into $`L^{}`$. We note that $`|M,\sigma `$ is equivalent to $`|M,m=0`$, where $`m=0`$ denotes the angular dependence of the spherical harmonic $`Y_{2m}`$. Choosing our coordinate system so that $`L`$ is along the $`z`$-axis and $`L^{}`$ is in the $`z`$-$`x`$ plane, we find that

$`M,\sigma |M,\sigma ^{}`$ $`=`$ $`D_{00}^{(2)}(0,\theta ,0)`$ (16)
$`=`$ $`(3\mathrm{cos}^2\theta 1)/2`$ (17)
where the $`D`$-term is a matrix element of the $`l=2`$ representation of the rotation group, and the second equality follows from the explicit formulas of the $`D`$-terms given in Ref. (Wigner, 1959). In summary, for $`\sigma `$-bonded ligands,

$$g_{\eta \eta \mathrm{`}}=h_\sigma h_\sigma ^{}(3\mathrm{cos}^2\theta 1)/2.$$
(18)
For the case of a $`\pi `$-bonded orbital $`L`$ and a $`\sigma `$-bonded orbital $`L^{}`$, we consider only the case in which the axis of the orbital lies along the circle connecting $`L`$ and $`L^{}`$; if it is perpendicular to this circle, $`L`$ and $`L^{}`$ have different inversion symmetries so their coupling vanishes. We write $`|\mathrm{\Psi }=h_\pi |M,\pi `$ where $`|M,\pi =(1/\sqrt{2})(|M,m=1+|M,m=1)`$ and $`m`$ the usual aximuthal angular momentum index for the spherical harmonics. (The “$``$” sign comes from the definition of the spherical harmonics). Then, again choosing a coordinate system in which $`L`$ is along the $`z`$-axis and $`L^{}`$ is in the $`z`$-$`x`$ plane, and following reasoning parallel to the $`\sigma `$-bonded case, we see that

$`M,\pi |M,\sigma ^{}`$ $`=`$ $`(1/\sqrt{2})\left[D_{10}^{(2)}(0,\theta ,0)+D_{10}^{(2)}(0,\theta ,0)\right]`$ (19)
$`=`$ $`\sqrt{3}\mathrm{sin}\theta \mathrm{cos}\theta .`$ (20)
and

$$g_{\eta \eta \mathrm{`}}=\sqrt{3}h_\pi h_\sigma ^{}\mathrm{sin}\theta \mathrm{cos}\theta .$$
(21)
Finally we turn to the case of two $`\pi `$-bonded orbitals. To obtain the subsequent results, it is sufficient to consider the case 1) in which the orbitals are parallel to the arc connecting $`L`$ and $`L^{}`$, and the case 2) in which they are perpendicular to it. By reasoning similar to that above, one sees that in case 1)

$`M,\pi |M,\pi ^{}`$ $`=`$ $`(1/2)\left[D_{11}^{(2)}(0,\theta ,0)+D_{11}^{(2)}(0,\theta ,0)D_{11}^{(2)}(0,\theta ,0)D_{11}^{(2)}(0,\theta ,0)\right]`$ (22)
$`=`$ $`\mathrm{cos}2\theta `$ (23)
and

$$g_{\eta \eta \mathrm{`}}=h_\pi h_\pi ^{}\mathrm{cos}2\theta .$$
(24)
In case 2),

$`M,\pi |M,\pi ^{}`$ $`=`$ $`(1/2)\left[D_{11}^{(2)}(0,\theta ,0)+D_{11}^{(2)}(0,\theta ,0)+D_{11}^{(2)}(0,\theta ,0)l+D_{11}^{(2)}(0,\theta ,0)\right]`$ (25)
$`=`$ $`\mathrm{cos}\theta `$ (26)
and

$$g_{\eta \eta \mathrm{`}}=h_\pi h_\pi ^{}\mathrm{cos}\theta .$$
(27)
Thus, combining Eqs. (13), (18), (21), (24), and (27), we find the following explicit form for $`E_{\mathrm{LFSE}}`$ as a sum of ligand-ligand interactions:

$$E_{\mathrm{LFSE}}=\sqrt{\underset{L,\eta ,L^{},\eta ^{}}{}}U(L,\eta ;L^{},\eta ^{})$$
(28)
where the ligand-ligand interaction $`U`$ is given as follows:

For $`\sigma `$-$`\sigma `$ interactions,

$$U(L,\eta ;L^{},\eta ^{})=(\alpha /5)e_\sigma (r)e_\sigma (r^{})u_{\sigma \sigma }(\theta )$$
(29)
for $`\pi `$-$`\sigma `$ interactions,

$$U(L,\eta ;L^{},\eta ^{})=(\alpha /5)e_\pi (r)e_\sigma (r^{})u_{\pi \sigma }(\theta )$$
(30)
and for $`\pi `$-$`\pi `$ interactions,

$$U(L,\eta ;L^{},\eta ^{})=(\alpha /5)e_\pi (r)e_\pi (r^{})u_{\pi \pi }(\theta )$$
(31)
where

$`e_\sigma (r)=`$ $`h_{\sigma }^{}{}_{}{}^{2}/(\epsilon _d\epsilon _{L,\eta }),`$ (32)
$`e_\pi (r)=`$ $`h_{\pi }^{}{}_{}{}^{2}/(\epsilon _d\epsilon _{L,\eta ^{}}),`$ (33)
$`u_{\sigma \sigma }(\theta )=`$ $`(9\mathrm{cos}^4\theta 6\mathrm{cos}^2\theta +1/5)/4,`$ (34)
$`u_{\pi \sigma }(\theta )=`$ $`(3\mathrm{cos}^4\theta +3\mathrm{cos}^2\theta 2/5),`$ (35)
$`u_{\pi \pi }(\theta )=`$ $`(4\mathrm{cos}^4\theta 3\mathrm{cos}^2\theta +1/5).`$ (36)
In each case, $`U`$ is given as a product of radial terms involving the two ligands, and a simple angular function. (Note that although denoted an interaction here, $`U`$ does not have units of energy because of the square root in Eq. (28)). In the calculation of the terms involving $`\pi `$-$`\sigma `$ interactions, each term involves a sum over two $`\pi `$-orbitals, parallel and perpendicular to the arc connecting the two ligands. In the latter case, $`g_{\eta \eta ^{}}`$ in Eq. (13) vanishes, but $`g_{\eta \eta }`$ and $`g_{\eta ^{}\eta ^{}}`$ do not. In the calculation of the terms involving $`\pi `$-$`\pi `$ interactions, one has a similar scenario except that one sums over two pairs of $`\pi `$-orbitals.

These forms are plotted out in Figure 1a. Note that the $`\sigma `$-$`\sigma `$ interactions are fairly similar in form to the $`\pi `$-$`\pi `$ interactions, both having pronounced minima at $`180^{}`$ (as well as the physically irrelevant one at $`0^{}`$), and a shallower minimum at $`90^{}`$. The $`\sigma `$-$`\pi `$ interaction is complementary to these, having minima at $`45^{}`$ and $`135^{}`$. We shall see later that these differences lead to large differences in ground-state structures of small clusters. For comparison, we show in Figure 1b an empirical angular interaction curve (Comba et al., 1995) assumed in some previous calculations of small transition-metal structures. The angular dependence is based on the observed square structure of small complexes of the transition metals of interest, and is proportional to $`\mathrm{sin}^22\theta `$. Since the metals which have square coordination generally have predominantly $`\sigma `$-bonds to their neighbors, the most relevant comparison is to the $`\sigma `$-$`\sigma `$ curve in Figure 1a. We see that the behavior is quite different. The empirical curve has equivalent minima at $`180^{}`$ and $`90^{}`$, while in the theoretical curve the $`180^{}`$ minimum is much deeper.

## III Tests of Functional Form

In order to evaluate the accuracy of the semiclassical form of Eq. (13) for the ligand-field stabilization energy, we have performed explicit tests for small clusters. These clusters consist of a central transition-metal ion with four neighbors placed at random orientations at random distances relative to the central ion. We consider only the minority-spin orbitals, as these determine $`E_{\mathrm{LFSE}}`$ for high-spin late transition metals. Four electrons are placed in these states, corresponding to Cu2+; other band filling values gives similar results. The random distances are taken into account by varying the couplings $`h_\sigma `$ and $`h_\pi `$ uniformly over a finite interval ranging from zero to $`h_{\mathrm{max}}^\sigma `$ or $`h_{\mathrm{max}}^\pi `$. We take $`h_{\mathrm{max}}^\pi =0.5h_{\mathrm{max}}^\sigma `$. The energies are obtained by explicit diagonalization of a tight-binding Hamiltonian for this cluster. All of the ligand orbitals are taken to have the same value of $`\epsilon _{L,\eta }`$. The only dimensionless variable that enters the results is then $`\gamma =|h_{\mathrm{max}}/(\epsilon _d\epsilon _{L,\eta })|`$. For small values of $`\gamma `$, the bonding is primarily ionic, and for larger values it acquires more covalent character. We use $`\gamma =0.1`$ for our results, but even when $`\gamma `$ is significantly larger we find that the accuracy of the semiclassical form is essentially the same.

We have fitted the energies of these clusters to the following semiclassical form for the energy:

$$E_{\mathrm{LFSE}}^2=\underset{L,\eta ,L^{},\eta ^{}}{}U(L,\eta ;L^{},\eta ^{})+\beta ,$$
(37)
involving two parameters $`\alpha `$ and $`\beta `$ (where $`U`$ contains $`\alpha `$). We use a database of 10,000 clusters to determine the parameters, and then test them on a set of 1000 clusters not included in the “training” set. Typical results are shown in Fig. 2, which shows the semiclassical energies vs. exact energies for 1000 clusters. Figure 2a corresponds to three $`\sigma `$-bonded ligands and one $`\pi `$-bonded one, while Fig. 2b corresponds to four $`\pi `$-bonded ligands. The rms errors in $`E_{\mathrm{LFSE}}`$ for these two cases are 9% and 14%, respectively. For clusters with two and three $`\pi `$-bonded ligands, the rms errors are 10% and 11%. This is to be compared with rms errors of 25%–30% that are obtained with empirical force fields (Carlsson, 1998).

## IV Small Cluster Minimum-Energy Structures

In this section, we describe some of the implications of the angular forms developed above for the structure of small model clusters consisting of a metal atom and four ligands. Such cluster calculations cannot treat the protein environment accurately; this would await parametrization and incorporation of the force field into protein codes, which is in progress in our group. However, from the small-cluster calculations it is possible to see the structural preferences of the angular forces by themselves. These are an important factor in the final structure adopted by a protein.

To keep the calculations as simple as possible, we place the ligands at frozen bond lengths from the central metal ion. Here “frozen means” that they are fixed at a given set of values, but these values are not necessarily the same for all of the ligands. The energy terms include the ligand-field stabilization energy as described by our angular terms, as well as a repulsive radial interaction between ligands forbidding close approaches. The latter term is needed because otherwise spurious structures involving $`45^{}`$ bond angles can appear when $`\pi `$-bonding ligands are present. For simplicity, we assume that the values of $`h_\sigma `$ and $`h_\pi `$ at a given distance are the same. To evaluate the magnitudes of these couplings, we assume a value of 100 kJ/mole $`=23.9`$ kcal/mole $`=1.04`$ eV for the ligand-field stabilization energy of a cluster with four $`\sigma `$-ligands in square coordination; this number would correspond to the more strongly stabilized complexes (Cotton, 1972). For clusters with unequal bond lengths, we assume an exponential decay, so that $`h_\sigma (r)=h_\sigma (r_0)\mathrm{exp}[\kappa (rr_0)]`$, where $`r_0`$ is the reference distance and $`\kappa `$ is a decay parameter. Since the most prominent case that we consider of a distant ligand is sulfur, we identify $`\kappa `$ with the spatial decay rate coming from the measured first ionization energy $`E_1`$ of sulfur, 10.4 eV (= 240 kcal/mole), using the formula $`\mathrm{}^2\kappa ^2/2m=E_1`$. This yields $`\kappa =1.65`$Å-1. The ligand-ligand terms contain an exponential term taken as the repulsive part of the van der Waals interactions as given in the “MM2” force field (Sprague et al., 1987). In our model clusters, we use parameters and bond lengths typical for copper or nickel interacting with nitrogen and/or sulfur ligands. This is because copper and nickel have the largest ligand-field energies among the 3d transition metals, and typical ligands for these metals are nitrogen and sulfur. We have examined four simple geometries; the corresponding minimum-energy clusters are shown in Fig. 3.

a) A cluster with all four ligands $`\sigma `$-bonded and placed at equal distances of 2.0Å. We take repulsion parameters appropriate for nitrogen. The result is a square-planar cluster (cf. Fig. 3a), consistent with the observed structures of many copper and nickel complexes. With this parameterization, the square structure is quite strongly favored over the tetrahedral one, by 0.31 eV or 7.0 kcal/mole. The square structure remains stable under increases of the ligand-ligand repulsion strengths. of up to a factor of over three. We note, however, that the stability of the square structure is strongly sensitive to the choice of repulsive terms. For example, we find that when the van der Waals term from the OPLS-II force field (Jorgensen et al., 1988) is used instead of the MM2 form, an increase of the repulsion strengths of only 20% is enough to stabilize the tetrahedral structure. We expect that input from fully quantum-mechanical total-energy calculations will be necessary to precisely pin down the competition between the electronic energy and the repulsion terms in determining the structural energy differences.

b) A cluster similar to that of a), but with one $`\pi `$-bonded ligand. This leads to a substantial deviation from planarity, as shown in Fig. 3b, where the $`\pi `$-bonded ligand is atom 4. As expected from the maximum of the $`\pi `$-$`\sigma `$ potential in Fig. 1, the angle between atom 4 and atoms 1–3 is greater than $`90^{}`$, about $`123^{}`$. The angle between the $`\sigma `$-bonded atoms 1, 2, and 3, is slightly greater than the ideal value of $`90^{}`$, because of the repulsive energy term.

c) A cluster with with two $`\sigma `$-bonded nitrogen ligands 1 and 2, one $`\sigma `$-bonded sulfur ligand 3, and a $`\pi `$-bonded sulfur ligand 4. The distances are 2.04Å for the nitrogen ligands, 2.18Å for the $`\pi `$-bonded sulfur ligand, and 2.64Å for the $`\sigma `$-bonded sulfur ligand. This geometry is motivated by the observed geometry of copper sites in “blue-copper” proteins. In these proteins, the nitrogen ligands and a cysteine sulfur ligand are close to the copper, and it is believed that the cysteine sulfur ligand is primarily $`\pi `$-bonded. An additional $`\sigma `$-bonded methionine sulfur ligand is farther from the copper. The values of the ligand distances here are taken from cluster calculations for blue-protein models (Ryde et al., 1996). The van der Waals parameters are taken from the OPLS parameter set but we do not have a reliable procedure for determining the differences between the ligand-ion interactions of the nitrogen and the sulfur ligands. Since the greater distance of the sulfur at 2.18Å from the copper will be compensated to some extent by the greater size of sulfur relative to nitrogen, we simply assume that the three close ligands have the same coupling strength to the central atom. We have varied the ratio of the sulfur to nitrogen coupling strengths by up to 30% in both directions and found changes of only a few degrees in the bond angles of the cluster. For the far ligand, we use the scaling procedure described above, which leads to a coupling that has 50% of the strength of the close ligands. We have varied this coupling from 0% to 75% of the close-ligand value, and again found bond-angle changes of only a few degrees. The lowest-energy structure for this cluster is indicated in Fig. 3c. The geometry is trigonal, with bond angles of $`125^{}`$ between the cysteine sulfur and the nitrogens; the nitrogen-nitrogen bond angle is $`91^{}`$. By comparison, the optimal values for protein models obtained by quantum-chemical calculations (Ryde et al., 1996) are $`125^{}`$ for the cysteine sulfur-nitrogen bond angle and $`103^{}`$ for the nitrogen-nitrogen bond angle. The results obtained by the present method are do of course not have quantitative accuracy, but the overall structure of the coordination shell is quite similar to that obtained in (Ryde et al., 1996), which is shown in Fig. 3d. The stability of this structure is not predetermined by the bond lengths that we used as input. To demonstrate this, we have considered a cluster with the same bond lengths as the blue-protein fragment, with with all ligands $`\sigma `$-bonded. In this case, the structure is planar, as shown in Fig. 3e. Thus the formation of the trigonal structure is directly related to the special character of the angular interactions associated with $`\pi `$-bonding.

## V Conclusion

We have seen that it is possible to extract the functional form of angular forces around transition metals in biomolecules by using an approximate treatment of the ligand-field Hamiltonian. The approximations made here, involving primarily ionic bonding, are complementary to those used in (Landis et al., 1998). The resulting angular forms have simple trigonometric forms, which are very different for $`\pi `$-bonding as compared to $`\sigma `$-bonding. Ligand-field energies are represented well by these angular forces, and much better than by force fields with assumed functional forms. Future work should aim to include the functional forms derived here in widely used biomolecular simulation packages.

###### Acknowledgements.

This work was supported by the National Science Foundation under Grant DMR-9971476. We appreciate informative conversations with Garland Marshall and Jay Ponder.

REFERENCES

* Allured, V. S., C. M. Kelly, and C. R. Landis. 1991. Shapes Empirical Force Field: New Treatment of Angular Potentials and its Application to Square-Planar Transition-Metal Complexes. J. Am. Chem. Soc. 113:1-12.
* Carlsson, A. E. 1998. Angular Forces Around Transition Metals in Biomolecules. Phys. Rev. Letters 81:477-480.
* Comba, P., T. W. Hambley, and M. Ströhle. 1995. The Directionality of d-Orbitals and Molecular-Mechanics Calculations of Octhedral Transition-Metal Compounds. Helv. Chim. Acta 78:2042-2047.
* Cotton, F. A., and G. Wilkinson. 1972. Advanced Inorganic Chemistry. Interscience, New York.
* Jorgensen, W. L. and J. Tirado-Rives. 1988. The OPLS Potential Functions for Proteins, Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 110:1657-1666.
* Kaim, W., and B. Schwederski. 1994. Bioinorganic Chemistry: Inorganic Elements in the Chemistry of Life. Wiley, New York.
* Landis, C. R., T. Cleveland, and T. K. Firman. 1998. Valence Bond Concepts Applied to the Molecular Mechanics Description of Molecular Shapes. 3. Applications to Transition Metal Alkyls and Hydrides. J. Am. Chem. Soc. 120:2641-2649.
* Ryde, U., M. H. M. Olsson, K. Pierllot, and B. O. Roos. 1996. The Cupric Geometry of Blue Copper Proteins is not Strained. J. Mol. Biol. 261:586-596.
* Sayle, D. C., M. A. Perrin, P. Nortier, and C. R. A. Catlow. 1995. Simulation Study of Copper(I) and Copper(II) Species in ZSM-5 Zeolite. J. Chem. Soc. Chem. Commun. 9:945-947.
* Sayle, D. C., et al., 1997. Computer Modeling of the Active-Site Configurations within the NO Decomposition Catalyst CU-ZSM-5. J. Phys. Chem. 101:3331-3337.
* Sprague, J. T., J. C. Tai, Y. Yuh, and N. L. Allinger. 1987. The MMP2 Calculational Method. J. Comput. Chem. 8:581-603.
* Timofeeva, T. V., J. H. Lii, and N. L. Allinger. 1995. Molecular Mechanics Explanation of the Metallocene Bent Sandwich Structure. J. Am. Chem Soc. 117:7452-7459.
* Wiesemann, F., S. Teipel, B. Krebs, and U. Howeler. 1994. Force-Field Calculations on the Structures of Transition-Metal Complexes: 1. Application to Copper(II) Complexes in Square-Planar Coordination. Inorg. Chem. 33:1891-1898.
* Wigner, E. P. 1959. Group Theory and its Application to the Quantum Mechanics of Atomic Spectra. Academic Press, New York.


format: text

text:
# PUZZLE OF THE CONSTANCY OF FUNDAMENTAL CONSTANTS

## 1 Introduction

Possible variability of fundamental physical constants is one of the topical problems of contemporary physics. The modern theory (Supersymmetric Grand Unification Theory – SUSY GUT, Superstring models, etc.) has established that the coupling constant values which characterize different kinds of interactions (i) are “running” with the energy transfer and (ii) may be different in different regions of the Universe and vary in the course of cosmological evolution (e.g., Ref. ). The energy dependence of the coupling parameters has been reliably confirmed by high-energy experiments (see, e.g., Ref. ), whereas the space-time variability of their low-energy limits so far escapes detection.

Note that a numerical value of any dimensional physical parameter depends on arbitrary choice of physical units. In turn, there is no way to determine the units in a remote space-time region other than through the fundamental constants. Therefore it is meaningless to speak of a variation of a dimensional physical constant without specifying which of the other physical parameters are defined to be invariable. Usually, while speaking of variability of a dimensional physical parameter, one implies that all the other fundamental constants are fixed. So did Milne  and Dirac  in their pioneering papers devoted to a possible change of the gravitational constant $`G`$. More recently, a number of authors considered cosmological theories with a time varying speed of light $`c`$ (e.g., Ref. ). However, if we adopt the standard definition of meter  as the length of path traveled by light in vacuum in 1/299 792 458 s, then $`c=\mathrm{2.997\hspace{0.17em}924\hspace{0.17em}58}\times 10^{10}\mathrm{cm}\mathrm{s}^1`$ identically. Similarly, one cannot speak of variability of the electron mass $`m_e`$ or charge $`e`$ while using the Hartree units ($`\mathrm{}=e=m_e=1`$), most natural in atomic physics. Thus, only dimensionless combinations of the physical parameters are truly fundamental, and only such combinations will be considered hereafter.

At present, the most promising candidate for the theory which is able to unify gravity with all other interactions is the Superstring theory, which treats gravity in a way consistent with quantum mechanics. All versions of the theory predict existence of the dilaton – a scalar partner to the tensorial graviton. Since the dilaton field $`\varphi `$ is generally not constant, the coupling constants and masses of elementary particles, being dependent on $`\varphi `$, should vary in space and time. Thus, the existence of a weakly coupled massless dilaton entails small, but non-zero, observable consequences such as Jordan–Brans–Dicke-type deviations from General Relativity and cosmological variations of the gauge coupling constants . These variations depend on cosmological evolution of the dilaton field and may be non-monotonous as well as different in different space-time regions.

In this paper, we focus on the space-time variability of the low-energy limits of two fundamental constants which are of paramount importance for atomic and molecular spectroscopy:
(i) the fine-structure constant $`\alpha =e^2/\mathrm{}c`$ (Sommerfeld parameter),
(ii) the electron-to-proton mass ratio $`\mu =m_e/m_p`$ (Born–Oppenheimer parameter).

The next section presents a compendium of the basic methods allowing one to obtain restrictions on possible variations of fundamental constants. In Sects. 3 and 4, we consider recent estimates of the values of $`\alpha `$ and $`\mu `$, respectively, at cosmological redshifts $`z=1`$–4 which correspond to epochs $`7`$–13 billion years ago. Conclusions are given in Sect. 5.

## 2 Tests of variability of fundamental constants

Techniques used to investigate time variation of the fundamental constants may be divided into extragalactic and local methods. The latter ones include astronomical methods related to the Galaxy and the Solar system, geophysical methods, and laboratory measurements.

### 2.1 Local tests

#### 2.1.1 Laboratory measurements

Laboratory tests are based on comparison of different frequency standards, depending on different combinations of the fundamental constants. Were these combinations changing differently, the frequency standards would eventually discord with each other. An interest in this possibility has been repeatedly excited since relative frequency drift was observed by several research groups using long term comparisons of different frequency standards. For instance, a comparison of frequencies of He-Ne/CH4 lasers, H masers, and Hg+ clocks with a Cs standard  has revealed relative drifts. Since the considered frequency standards have a different dependence on $`\alpha `$ via relativistic contributions of order $`\alpha ^2`$, the observed drift might be attributed to changing of the fine-structure constant. However, the more modern was the experiment, the smaller was the drift. Taking into account that the drift may be also related to some aging processes in experimental equipment, Prestage et al.  concluded that the current laboratory data provide only an upper limit $`|\dot{\alpha }/\alpha |3.7\times 10^{14}\mathrm{yr}^1`$.

#### 2.1.2 Analysis of the Oklo phenomenon

The most stringent limits to variation of the fine-structure constant $`\alpha `$ and the coupling constant of the strong interaction $`\alpha _s`$ have been originally inferred by Shlyakhter  from results of an analysis of the isotope ratio 149Sm/147Sm in the ore body of the Oklo site in Gabon, West Africa. This ratio turned out to be considerably lower than the standard one, which is believed to have occurred due to operation of the natural uranium fission reactor about $`2\times 10^9`$ yr ago in those ores. One of the nuclear reactions in the fission chain was the resonance capture of neutrons by 149Sm nuclei. Actually, the rate of the neutron capture reaction is sensitive to the energy of the relevant nuclear resonance level $`E_r`$, which depends on the strong and electromagnetic interaction. Since the capture has been efficient $`2\times 10^9`$ yr ago, in means that the position of the resonance has not shifted by more than it width (very narrow) during the elapsed time. At variable $`\alpha `$ and invariable $`\alpha _s`$ (which is just a model assumption), the shift of the resonance level would be determined by changing the difference between the Coulomb energies of the ground-state nucleus 149Sm and the nucleus 150Sm excited to the level $`E_r`$. Unfortunately, there is no experimental data for the Coulomb energy of the excited 150Sm in question. Using order-of-magnitude estimates, Shlyakhter  concluded that $`|\dot{\alpha }/\alpha |10^{17}\mathrm{yr}^1`$. From an opposite model assumption that $`\alpha _s`$ is changing whereas $`\alpha =`$constant, he derived a bound $`|\dot{\alpha }_s/\alpha _s|10^{19}\mathrm{yr}^1`$.

Damour and Dyson  performed a more careful analysis, which resulted in the upper bound $`|\dot{\alpha }/\alpha |7\times 10^{17}\mathrm{yr}^1`$. They have assumed that the Coulomb energy difference between the nuclear states of 149Sm and 150Sm in question is not less than that between the ground states of 149Sm and 150Sm. The latter energy difference has been estimated from isotope shifts and equals $`1`$ MeV. However, it looks unnatural that a weakly bound neutron ($`0.1`$ eV), captured by a 149Sm nucleus to form the highly excited state 150Sm, can so strongly affect the Coulomb energy. Moreover, heavy excited nuclei often have Coulomb energies smaller than those for their ground states (e.g., Ref. ). This indicates the possibility of violation of the basic assumption involved in Ref. , and therefore this method may possess a lower actual sensitivity. Furthermore, a correlation between $`\alpha `$ and $`\alpha _s`$ (which is likely in the frame of modern theory) might lead to considerable softening of the above-mentioned bound, as estimated by Sisterna and Vucetich .

#### 2.1.3 Some other local tests

Geophysical, geochemical, and paleontological data impose constraints on a possible changing of various combinations of fundamental constants over the past history of the Solar system, however most of these constraints are very indirect. A number of other methods are based on stellar and planetary models. The radii of the planets and stars and the reaction rates in them are influenced by values of the fundamental constants, which offers a possibility to check variability of the constants by studying, for example, lunar and Earth’s secular accelerations. This was done using satellite data, tidal records, and ancient eclipses. Another possibility is offered by analyzing the data on binary pulsars and the luminosity of faint stars. Most of these have relatively low sensitivity. Their common weak point is the dependence on a model of a fairly complex phenomenon, involving many physical effects.

An analysis of natural long-lived $`\alpha `$\- and $`\beta `$-decayers in geological minerals and meteorites is much more sensitive. For instance, a strong bound, $`|\dot{\alpha }/\alpha |<5\times 10^{15}\mathrm{yr}^1`$, was obtained by Dyson  from an isotopic analysis of natural $`\alpha `$\- and $`\beta `$-decay products in Earth’s ores and meteorites.

Having critically reviewed the wealth of the local tests, taking into account possible correlated synchronous changes of different physical constants, Sisterna and Vucetich  derived restrictions on possible variation rates of individual physical constants for ages $`t`$ less than a few billion years ago, which correspond to cosmological redshifts $`z0.2`$. In particular, they have arrived at the estimate $`\dot{\alpha }/\alpha =(1.3\pm 6.5)\times 10^{16}\mathrm{yr}^1`$.

All the local methods listed above give estimates for only a narrow space-time region around the Solar system. For example, the epoch of the Oklo reactor ($`1.8\times 10^9`$ years ago) corresponds to the cosmological redshift $`z0.1`$. These tests cannot be extended to earlier evolutionary stages of the Universe, because the possible variation of the fundamental constants is, in general, unknown and may be oscillating . Another investigation is needed for higher cosmological redshifts.

### 2.2 Extragalactic tests

Extragalactic tests, in contrast to the local ones, concern values of the fundamental constants in distant areas of the early Universe. A test which relates to the earliest epoch is based on the standard model of the primordial nucleosynthesis. The amount of 4He produced in the Big Bang is mainly determined by the neutron-to-proton number ratio at the freezing-out of n$``$p reactions. The freezing-out temperature $`T_f`$ is determined by the competition between the expansion rate of the Universe and the $`\beta `$-decay rate. A comparison of the observed primordial helium mass fraction, $`Y_p=0.24\pm 0.01`$, with a theoretical value allows one to obtain restrictions on the difference between the neutron and proton masses at the epoch of the nucleosynthesis and, through it, to estimate relative variation of the curvature radius $`R`$ of extra dimensions in multidimensional Kaluza–Klein-like theories which in turn is related to the $`\alpha `$ value . However, as noted above, different coupling constants might change simultaneously. For example, increasing the constant of the weak interactions $`G_F`$ would cause a weak freezing-out at a lower temperature, hence a decrease in the primordial 4He abundance. This process would compete with the one described above, therefore, it reduces sensitivity of the estimates. Finally, the restrictions would be different for different cosmological models since the expansion rate of the Universe depends on the cosmological constant $`\mathrm{\Lambda }`$.

The most unambiguous estimation of the atomic and molecular constants at early epochs and in distant regions of the Universe can be performed using the extragalactic spectroscopy. Accurate measurements of the wavelengths in spectra of distant objects provide quantitative constraints on the variation rates of the physical constants. This opportunity has been first noted and used by Savedoff , and in recent years exploited by many researchers (see, e.g., Refs.  and references therein). At present, the extragalactic spectroscopy enables one to probe the physical conditions in the Universe up to cosmological redshifts $`z4`$, which correspond, by order of magnitude, to the scales $`10^{10}`$ yr in time and $`10^9`$ parsec in space. In the following sections, we review briefly the studies of the space-time variability of the fine-structure constant $`\alpha `$ and the electron-to-proton mass ratio $`\mu `$, based on the latter method.

## 3 Non-variability of $`\alpha `$

We have already mentioned in Sect. 2.1.1 that several laboratory tests hinted at a tentative time variation of $`\alpha `$, but were later refuted by measurements at a higher level of accuracy. A similar situation has occurred for extragalactic tests at larger space-time scales.

Bahcall and Schmidt  were the first to use spectral observations of distant quasars to set a bound on the variability of the fine-structure constant. They have obtained an estimate $`\mathrm{\Delta }\alpha /\alpha =(2\pm 5)\times 10^2`$ at $`z=1.95`$. Later statistical analyses  of fine-structure doublet lines in quasar spectra appeared to indicate a tentative variation of $`\alpha `$ (of the order of $`0.3\%`$ at the cosmological redshift $`z2`$). However, this tentative variation has been shown to result from a statistical bias .

Another statistical examination of the fine-doublet wavelengths of absorption lines in quasar spectra  indicated a tentative (at the 2–3$`\sigma `$ level) variability of $`\alpha `$ values by $`0.1\%`$ over the celestial sphere (as function of angle) at redshifts $`z2`$–3. However, this result has not been confirmed by a later analysis , which was based on higher-quality spectra and yielded an order of magnitude higher precision.

Quite recently, Webb et al.  have estimated $`\alpha `$ by comparing wavelengths of Fe ii and Mg ii fine-splitted spectral lines in extragalactic spectra and in the laboratory. Their result suggests a time-variation of $`\alpha `$ at the incredibly high accuracy level of $`10^3\%`$: the authors’ estimate reads $`\mathrm{\Delta }\alpha /\alpha =(1.9\pm 0.5)\times 10^5`$ at $`z=1.0`$–1.6. Note, however, two important sources of a possible systematic error which could mimic the effect: (a) Fe ii and Mg ii lines used are situated in different orders of the echelle-spectra, so relative shifts in calibration of the different orders can affect the result of comparison, and (b) if the isotopic composition varies during the evolution of the Universe, then the average doublet separations should vary due to the isotopic shifts. Were the relative abundances of Mg isotopes changing during the cosmological evolution, the Mg ii lines would be subject to an additional $`z`$-dependent shift relative to the Fe ii lines, quite sufficient to simulate the variation of $`\alpha `$ (this shift can be easily estimated from recent laboratory measurements ).

The method based on the fine splitting of a line of the same ion species is not affected by these two uncertainty sources. We have studied the fine splitting of the doublet lines of Si iv, C iv, Mg ii and other ions, observed in spectra of distant quasars. According to quantum electrodynamics, the relative splitting of these lines $`\delta \lambda /\lambda `$ is proportional to $`\alpha ^2`$ (neglecting small relativistic corrections, recently estimated by Dzuba et al. ). We have selected the results of high-resolution observations , most suitable for an analysis of the variation of $`\alpha `$. According to our analysis, presented elsewhere , the most reliable estimate of the possible deviation of the fine-structure constant at $`z=2`$–4 from its present ($`z=0`$) value:

$$\mathrm{\Delta }\alpha /\alpha =(4.6\pm 4.3[\mathrm{stat}]\pm 1.4[\mathrm{syst}])\times 10^5.$$
(1)
Thus, only an upper bound can be derived at present for the long-term variability of $`\alpha `$:

$$|\dot{\alpha }/\alpha |<1.4\times 10^{14}\mathrm{yr}^1$$
(2)
(at the 95% confidence level).

## 4 Non-variability of $`\mu `$

The dimensionless Born–Oppenheimer constant $`\mu =m_e/m_p`$ approximately equals the ratio of the constant of electromagnetic interaction $`\alpha =e^2/\mathrm{}c1/137`$ to the constant of strong interaction $`\alpha _s=g^2/\mathrm{}c14`$, where $`g`$ is the effective coupling constant calculated from the amplitude of $`\pi `$-meson–nucleon scattering at low energy.

An early limit on the possible variation of this constant, $`|\dot{\mu }/\mu |<1.2\times 10^{10}\mathrm{yr}^1`$, has been derived from the concordance of K–Ar and Rb–Sr geochemical ages . The first astrophysical bound , based on the agreement between redshifts of atomic hydrogen and other lines in quasar absorption spectra, turned out to be twice stronger.

Orders-of-magnitude more precise analysis has become possible due to discovery  of a system of H2 absorption lines in the spectrum of quasar PKS 0528$``$250 at $`z=2.811`$. A study of this system yields information about physical conditions and, in particular, the value of $`\mu `$ at this redshift (corresponding to the epoch when the Universe was several times younger than now). A possibility of distinguishing between the cosmological redshift of spectral wavelengths and shifts due to a variation of $`\mu `$ arises from the fact that the electronic, vibrational, and rotational energies of H2 each undergo a different dependence on the reduced mass of the molecule. Hence comparing ratios of wavelengths $`\lambda _i`$ of various H2 electron-vibration-rotational lines in a quasar spectrum at some redshift $`z`$ and in laboratory (at $`z=0`$), we can trace variation of $`\mu `$. We have calculated  sensitivity coefficients $`K_i`$ of the wavelengths $`\lambda _i`$ with respect to possible variation of $`\mu `$ and applied a linear regression analysis to the measured redshifts of individual lines $`z_i`$ as function of $`K_i`$. If the proton mass in the epoch of line formation were different from the present value, the measured $`z_i`$ and $`K_i`$ values would correlate:

$$\frac{z_i}{z_k}=\frac{(\lambda _i/\lambda _k)_z}{(\lambda _i/\lambda _k)_0}1+(K_iK_k)\left(\frac{\mathrm{\Delta }\mu }{\mu }\right).$$
(3)
We have performed a $`z`$-to-$`K`$ regression analysis using a modern high-resolution spectrum of PKS 0528$``$250. Eighty-two of the H2 lines have been identified. The resulting parameter estimate and $`1\sigma `$ uncertainty is

$$\mathrm{\Delta }\mu /\mu =(11.5\pm 7.6[\mathrm{stat}]\pm 1.9[\mathrm{syst}])\times 10^5.$$
(4)
The $`2\sigma `$ confidence bound on $`\mathrm{\Delta }\mu /\mu `$ reads

$$|\mathrm{\Delta }\mu /\mu |<2.0\times 10^4.$$
(5)
Assuming that the age of the Universe is $`1.5\times 10^{10}`$ yr the redshift of the H2 absorption system $`z=2.81080`$ corresponds to the elapsed time $`1.3\times 10^{10}`$ yr (in the standard cosmological model). Therefore we arrive at the restriction

$$|\dot{\mu }/\mu |<1.5\times 10^{14}\mathrm{yr}^1$$
(6)
on the variation rate of $`\mu `$, averaged over 90% of the lifetime of the Universe.

## 5 Conclusions

Despite the theoretical prediction that fundamental constants of Nature should vary, no statistically significant variation of any of the constants has been reliably detected up to date, according to our point of view substantiated above. The upper limits obtained indicate that the constants of electroweak and strong interactions did not significantly change over the last 90% of the history of the Universe. The striking tightness of these limits is really astonishing and has already ruled out some theoretical models (see Refs. ). A more elaborated theory (e.g., Ref. ) cannot be ruled out yet, but its parameters can be severely restricted (e.g., see Ref. ). This shows that more precise measurements and observations and their accurate statistical analyses are required in order to detect the expected variations of the fundamental constants.

Acknowledgments. This work was performed in frames of the Russian State Program “Fundamental Metrology” and supported by the grant RFBR 99-02-18232.


text:
# Still flat after all these years!

## Abstract

The Universe could be spatially flat, positively curved or negatively curved. Each option has been popular at various times, partly affected by an understanding that models tend to evolve away from flatness. The curvature of the Universe is amenable to measurement, through tests such as the determination of the angles of sufficiently large triangles. The angle subtended by the characteristic scale on the Cosmic Microwave sky provides a direct test, which has now been realised through a combination of exquisite results from a number of CMB experiments.

After a long and detailed investigation, with many false clues, it seems that the mystery of the curvature of the Universe is now solved. It’s an open and shut case: the Universe is flat!

Humankind has long pondered whether the world is flat or round – whether boats would fall off the Earth and whether one can reach the edge of Creation. After the coming of General Relativity, the idea of a non-Euclidean geometry in 3 spatial dimensions began to be applied to the entire Universe. We imagine space being curved in an imaginary 4th dimension, and that this curvature can be detected by sufficiently accurate measurements in the 3 dimensions in which we live.

The classic example of a test for non-flatness involves measuring the angles of a large triangle and determining whether they sum to $`180^{}`$, or to some larger or smaller number. We know that in the 2-d analogy, on the surface of the Earth for example, the angles of spherical triangles always come to more than two right angles. And for triangles drawn on the inside of a trumpet horn, an example of a 2-d hyperbolic space, the angles sum to less than $`180^{}`$. The same test can be applied in 3-d, and, in addition, we can also consider surface area or volume tests. The trick is to apply such tests on a sufficiently large scale, since space is locally extremely flat. As we shall see, this very test has now been applied using CMB anisotropies, on a distance scale which is approximately the Hubble radius.

Let us define some terminology. In an isotropic and homogeneous Universe Einstein’s equations imply the following evolution equation for the scale factor $`R(t)`$:

$$H^2\left(\frac{\dot{R}}{R}\right)=\frac{8\pi G\rho }{3}\frac{c^2}{^2}$$
(1)
where $`\rho `$ is the total energy density and $``$ is the radius of curvature of the Universe (we write $`i`$ for open models, and $`\mathrm{}`$ for flatness). Conventionally curvature is measured through the density parameter $`\mathrm{\Omega }`$, which is defined by $`\mathrm{\Omega }8\pi G\rho /3H^2`$. This relates to curvature because $`\mathrm{\Omega }1=c^2/H^2^2`$. If $`\mathrm{\Omega }=1`$ then the Universe is flat, while $`\mathrm{\Omega }>1`$ and $`\mathrm{\Omega }<1`$ imply a closed and open geometry, respectively. There are in principle many ways to measure the total energy density of the Universe, but most methods do this indirectly, for example by measuring the density in different species (obeying different equations of state) separately and then summing to infer the value of $`\mathrm{\Omega }`$.

From the time of Euclid to the time of Einstein, scientists had a firm understanding of flat geometry, and never doubted that this described the Universe. Soon after curved Universes were proposed mathematically, closed models became the most popular, at least theoretically. Although there was never good evidence for believing that $`\mathrm{\Omega }`$ could be significantly bigger than unity, nevertheless there was some allure for a Universe with finite volume that also had a very definite future. $`\mathrm{\Omega }`$ is rather difficult to measure, and it wasn’t until perhaps the 1960s that open Universes gained observational favour. Throughout the 1970s and 1980s there was a split between theoretical and observational cosmologists. Observers saw only evidence for relatively low values of $`\mathrm{\Omega }`$, while theorists preferred to imagine that there might be enough dark matter to preserve the simplicity of $`\mathrm{\Omega }=1`$. Theorists continued to hope for $`\mathrm{\Omega }_{\mathrm{matter}}=1`$ even when there was considerable evidence to the contrary. There was a short period in the mid 1990s when open universes were popular, and even theorists got excited by inflationary open models. Now it appears that after many false trails we are once again back on the road to flatness.

Flatness does of course imply extreme fine tuning in the early Universe , since $`\mathrm{\Omega }=1`$ is unstable. There has been much philosophical debate over whether it is a serious problem that the Universe can have contrived to avoid diverging much from flatness. On the one hand there are those who contend that since $`\mathrm{\Omega }=\mathrm{constant}`$ has only 3 solutions, 0, 1 and $`\mathrm{}`$, then faced with that choice there is only one reasonable option. Others argue that this only makes sense if one invokes a mechanism for fine tuning the Universe to $`\mathrm{\Omega }=1\pm ϵ`$ at early times, and that inflation comes to the rescue. Yet others would argue that the anthropic principle saves us from living in Universes where $`\mathrm{\Omega }`$ departs too far from unity. We prefer to take the pragmatic approach – let us wait to see what the curvature of the Universe turns out to be, and then worry about the explanation. It appears that the wait is now over. $`\mathrm{\Omega }`$ is indeed very close to unity. So now we can worry about what it all means!

The newly applied test, using Cosmic Microwave Background anisotropies, has the appeal of being a much more direct measurement of $`\mathrm{\Omega }`$, regardless of the specific kinds of energy or matter which contribute to the $`\mathrm{\Omega }`$ census. The only loop-hole is that one must first be confident that the test is being applied within the correct cosmological model. Once that has been established, then $`\mathrm{\Omega }`$ can be constrained by looking at the position of the peak in the CMB power spectrum, which is expected to exist in most popular cosmological models. CMB anisotropies are measured by looking at the variation of mean square temperature fluctuation with angular scale – conventionally one uses the power spectrum versus multipole $`\mathrm{}`$, where $`\mathrm{}1/\theta `$. Some examples of power spectra are shown in Fig. 1.

The CMB triangle test involves using a standard ruler on the last scattering surface and then attempting to measure its angular size. The ruler here is the sound horizon at last scattering, corresponding to the size of the fundamental acoustic mode at that epoch. The lengths of the two longest sides of the triangle are the distance to the last scattering surface, which is essentially fixed. With all 3 sides of the triangle having their lengths established, all that remains is to measure the angle subtended by the shortest side, the standard ruler on the last scattering surface, and we have determined whether or not the Universe is flat .

In the most popular models, $`\mathrm{}_{\mathrm{peak}}220`$ in a flat Universe. If $`\mathrm{\Omega }`$ is changed, the main effect is that the location of the peak shifts according to $`\mathrm{}_{peak}\mathrm{\Omega }^{1/2}`$. Of course the height of the peak is also affected, but in a manner which strongly depends on other parameters, while the position of this main acoustic peak is a robust measure of $`\mathrm{\Omega }`$ alone.

In order to apply the CMB flatness test, it is important first to establish the framework in which to perform the test. There are 3 separate classes of model for generating cosmological perturbations:

* adiabatic initial conditions;
* isocurvature initial conditions;
* topological defect sources.

As shown in Fig. 1, all of these ideas for the origin of perturbations can produce at least one peak in the power spectrum. To understand the basic differences, imagine a scenario in which fluctuations are laid down at very early times in a way which is apparently acausal. In other words, modes over all relevant length scales begin with the same phase, irrespective of whether there has apparently been enough time for the establishment of causal contact. Inflation carries this off by having the Universe effectively expand faster than the speed of light  – but any initial condition synchronized over a sufficiently large volume would be equivalent. Such perturbations can then be either of the adiabatic or isocurvature type or, in principle, a combination of the two. ‘Adiabatic’ means that the entropy per species is unperturbed, while ‘isocurvature’ means that the total energy density is unperturbed. CMB anisotropies can be thought of in terms of a driven harmonic oscillator  with the driving force being the gravitational perturbation. For adiabatic modes, we have a driving force which already exists when the modes start to oscillate, shortly after they come inside the horizon, whereas for isocurvature modes the driving force is growing from zero. The phases of the oscillating fluid are captured in a snapshot at the last scattering epoch and hence we find that the power spectra of adiabatic and isocurvature modes are in anti-phase with each other (Fig. 1).

Topological defect models generate density perturbations as the field re-orders itself on roughly the horizon scale. These modes are therefore causal. We can think of a model like cosmic strings as arising from a number of patches on the sky, each of which is like an isocurvature model which chose some random number for the initial phase of the modes. The final power spectrum is the incoherent average of a number of shifted isocurvature power spectra, and thus shows a single broad peak (Fig. 1).

In Fig. 1 we show an example of a power spectrum from each of the three main families of models that have been considered for generating structure. We have specifically chosen models dominated by a cosmological constant (to be consistent with currently popular ideas). The curves have been scaled to have the same peak height. These three generic predictions look quite different. Note that it is possible to argue for other sources of anisotropy at the smallest $`\mathrm{}`$s, and one should focus one’s attention at the higher multipoles.

Evidence has been accumulating from CMB experiments, as well as the relative normalization with the matter fluctuation power spectrum, which indicate that the fluctuations are adiabatic. Some individual CMB data sets give evidence for a localised peak . However, it is the compilation of all available data that gives a convincing demonstration. To this end, we have calculated a binned anisotropy power spectrum, which we plot in Fig. 2 . The resulting power spectrum is flat at large angles, with a gradual rise to a prominent peak a little below $`1^{}`$, with a clear decrease thereafter. This is precisely the shape predicted by inflationary-inspired adiabatic models.

Although we can imagine that small components of isocurvature or defect-produced perturbations may still be possible , it is apparent is that overall the adiabatic models are in very good shape. Therefore, we can assume the adiabatic family, and use the position of the peak in $`\mathrm{}`$ to constrain the geometry of the Universe within the framework of adiabatic models. We do this rather simply by taking a cosmological constant dominated model, of the sort which fits a wide range of present cosmological data, and we rescale the $`\mathrm{}`$ scale by $`\mathrm{\Omega }^{1/2}`$ to find a range of acceptable values of $`\mathrm{\Omega }`$. This provides a 95% confidence region $`0.79<\mathrm{\Omega }<1.17`$, or approximately $`\mathrm{\Omega }=1.0\pm 0.1`$. It thus appears to be settled that the Universe is very close to flat – much closer than would be implied by the amount of visible or dark matter inferred from studies of galaxies, typically $`\mathrm{\Omega }_{\mathrm{matter}}0.3`$.

The next wave of CMB experiments should determine more definitively whether the perturbations are adiabatic, by measuring the power spectrum at angular scales where oscillations are predicted. The position of the first peak should soon be determined more precisely. In addition, other tests, including those from distant supernovae, should help determine the constituents that make up $`\mathrm{\Omega }`$.

It may be that one cosmic mystery has been solved, and it is tempting to regard this as something akin to a proof that the Universe once went through a period of inflation. However, new mysteries have been uncovered: what exactly was the physics which led to these adiabatic perturbations; what is this ‘dark energy’ which makes the Universe flat? and why did the Universe contrive to contain baryons, dark matter and dark energy in similar proportions?

###### Acknowledgements.

This research was supported by the Natural Sciences and Engineering Research Council of Canada, and by the US National Science Foundation.