text:
# References

The low energy dynamics of the globally supersymmetric extension of four dimensional quantum Yang-Mills theory with classical gauge group, $`𝒩=2`$ supercharges, and $`N_f`$ flavours of fundamental matter hypermultiplets, at weak coupling, has been found to be exactly solvable, using methods based on the pioneering work of Seiberg and Witten in . However, the forms of the exact solutions proposed \[3–8\] for theories with gauge group $`SU(N)`$ differ from each other. Furthermore, the results for the low energy effective actions derived from these exact solutions were found not to agree with the instanton predictions in these theories for $`N_f=2N`$ \[14–16, 19\]. There are also discrepancies in the expressions for the quantum moduli for $`N<N_f<2N`$, as reported in .
Recently a matching prescription to resolve these discrepancies was proposed in , generalising other proposals following from one-instanton and two-instanton checks of the conjectured exact results \[10–19\] for the matching between the two sets of predictions.
By considering the permissible non-perturbative redefinitions of the physical quantities involved, the scheme eliminates the reported ambiguities between the proposed exact solutions themselves and between the results from these and their instanton counterparts. Hence the exact results and instanton predictions are in agreement modulo reparameterisations. The constants in the matching scheme can only be fixed by comparison with instanton calculations, however, and cannot be derived from the exact results themselves.
In this letter we perform a non-trivial test of the matching scheme in  for the entire quantum moduli space in $`SU(N)`$ theories with $`N_f=2N`$ massless fundamental hypermultiplets by determining the constants which give agreement between the exact results and the one-instanton predictions for general $`N`$. Previously this check had been performed for a single special point of the moduli space in these theories , where agreement was found within the scheme.
We begin by briefly reviewing the matching scheme  and the conjectured method of exact solution for $`𝒩=2`$ supersymmetric $`SU(N)`$ gauge theories.
The coupling parameter (valid near weak coupling) used throughout is $`q\mathrm{exp}(2\pi i\tau )`$, where

$$\tau =\frac{\theta }{2\pi }+\frac{4\pi i}{g^2}$$
(1)
is the complex gauge coupling constant, for coupling $`g`$. In the weak coupling regime, $`q0`$. The scheme for matching the conjectured exact solutions for the low energy theories to the instanton results is derived by considering the most general mapping between the parameters and vacuum expectation values (VEVs) of each set of results which are consistent at weak coupling and obey the constraints imposed by supersymmetry . It agrees with dimensional analysis and also ensures that the notion of the moduli space is preserved.
For $`𝒩=2`$ supersymmetric $`SU(N)`$ theories, the low energy effective action is a function of the bare masses $`\{m_n\}`$, $`n=1\mathrm{}2N`$ and the coupling parameter $`q`$, and the set of VEVs is $`\{u_i\}`$, $`i=1\mathrm{}N`$, all of which assume complex values. We denote the parameters and VEVs appearing in the proposed exact solutions of these theories $`\{\stackrel{~}{q},\stackrel{~}{m}_n,\stackrel{~}{u}_i\}`$, and those appearing in the instanton predictions $`\{q,m_n,u_i\}`$. The constants in the matching scheme are denoted $`\{C_s,B_s,A_s^{(i;\{i_m\})}\}`$. Holomorphy and the asymptotic behaviour at weak coupling imply that the general map between $`q`$ and $`\stackrel{~}{q}`$ is

$$\stackrel{~}{q}=\underset{s=0}{\overset{\mathrm{}}{}}C_sq^{s+1}.$$
(2)
Dimensionless ratios of the masses cannot enter into Eq. (2) due to the matching that must hold at very weak coupling. The masses $`\{m_n\}`$ and $`\{\stackrel{~}{m}_n\}`$ are related via

$$\stackrel{~}{m}_n=\left(1+\underset{s=1}{\overset{\mathrm{}}{}}B_sq^s\right)m_n.$$
(3)
Finally, the matching relation between the VEVs $`\{u_i\}`$ and $`\{\stackrel{~}{u}_i\}`$ is given by

$$\stackrel{~}{u}_i=\underset{s=0}{\overset{\mathrm{}}{}}A_s^{(i;\{i_m\})}q^s\underset{m=0}{\overset{N}{}}u_{i_m},$$
(4)
in which we define $`u_0m`$ .
We shall employ the relations in Eqs. (2, 4) to one-instanton level, or, equivalently, to $`𝒪(q)`$, in the test of the matching scheme. The matching in Eq. (3) between the masses is not required here since all hypermultiplet matter has zero mass.
The defining quantity in $`𝒩=2`$ supersymmetric $`SU(N)`$ theories is the prepotential, $``$, a function of the superfields, which determines the low-energy effective (Wilsonian) action of the theory.
The prepotential can be decomposed into a classical part ($`_{\mathrm{cl}}`$) and its perturbative corrections ($`_{1\mathrm{loop}}`$), which are one-loop exact in this case (due to nonrenormalisation theorems), and non-perturbative corrections, ($`_{\mathrm{inst}}`$), containing instanton effects of all orders:

$$=_{\mathrm{cl}}+_{1\mathrm{loop}}+_{\mathrm{inst}}.$$
(5)
The instanton contributions have the following form:

$$_{\mathrm{inst}}=\underset{n=0}{\overset{\mathrm{}}{}}q^n_n,$$
(6)
where $`_n=_n(a_i)`$ are functions of the quantumVEVs, $`\{a_i\}`$, of the scalar superfield $`\varphi `$ in the adjoint representation (i.e., the Higgs field), and $`n`$ is the instanton number. The superfield $`\varphi `$ is a member of the vector multiplet of the theory. In general, the prepotential $``$ is a holomorphic function of the VEVs $`\{a_i\}`$. In the theories considered here, the perturbative beta function vanishes, and so $`_{1\mathrm{loop}}=0`$. The VEVs of the electric-magnetic dual of $`\varphi `$ are $`\{a_{D,i}\}`$, and are related to the prepotential via

$$a_{D,i}=\frac{}{a_i}.$$
(7)
The gauge coupling matrix of the theory is given by

$$\tau _{ij}=\frac{^2}{a_ia_j}.$$
(8)
The scalar potential of the $`𝒩=2`$ supersymmetric Lagrangian describing the vector multiplet is a function of $`\varphi `$. At weak coupling, the VEV of $`\varphi `$ is the matrix

$$<\varphi >=\mathrm{diag}[a_1,\mathrm{},a_N].$$
(9)
The region of the quantum moduli space corresponding to the VEVs of scalar fields from the vector multiplet of the theory is referred to as the Coulomb branch, hence we are concerned only with the Coulomb branch (or phase) of these theories. The VEVs, or quantum moduli, $`\{u_i\}`$, of the entire theory can be written in terms of the VEVs $`\{a_i\}`$ of $`\varphi `$:

$$u_n=<\mathrm{tr}(\varphi ^n)>=(1)^n\underset{i_1<\mathrm{}<i_N}{}a_{i_1}\mathrm{}a_{i_N}.$$
(10)
For $`SU(N)`$ theories, classically one has $`\mathrm{tr}(\varphi )=_{i=1}^Na_i=0`$, and $`u_1=0`$, by definition. In particular, the first non-zero classical VEV is given by

$$u_2^{\mathrm{cl}}=<\mathrm{tr}(\varphi ^2)>=\frac{1}{2}\underset{i=1}{\overset{N}{}}a_i^2.$$
(11)
The classical moduli space of these theories does not receive quantum corrections, as it is protected by $`𝒩=2`$ supersymmetry, but the metric $`(ds^2=\mathrm{Im}\tau _{ij}\mathrm{d}a_i\mathrm{d}\overline{a}_j)`$ on it does. Thus the VEVs $`\{a_i\}`$ receive quantum corrections. Following the methods introduced in , one identifies the $`(N1)`$-dimensional quantum moduli space with the moduli space of a genus $`(N1)`$ compact Riemann surface. Then the functions $`\{a_i,a_{D,i}\}`$ can be calculated as the periods (about certain cycles) of the Riemann surface, and the gauge coupling matrix $`\tau _{ij}`$ (Eq. (8)) is the period matrix of such a surface. For comprehensive reviews of the exact results, see, for example, references \[20–22\].
A standard result of the theory of algebraic curves  is that any compact Riemann surface can be specified completely by a class of elliptic ($`N12`$) or hyperelliptic ($`N1>2`$) curves. For the moduli spaces considered here, these curves have the form

$$y^2=F^2(x)G(x),$$
(12)
where $`F(x)`$ is a polynomial of degree $`(N1)`$ in the dummy variable $`x`$, whose coefficients are functions of the set of VEVs $`\{u_i\}`$. The roots of $`F(x)=0`$ are the exact VEVs $`\{e_i\}`$, which parameterise the moduli space, and are related to the VEVs $`\{a_i\}`$ by

$$a_n=(1)^n\underset{i_1<\mathrm{}<i_N}{}e_{i_1}\mathrm{}e_{i_N}.$$
(13)
These parameters obey $`_{i=1}^Ne_i=0`$.
The first non-zero quantum VEV in the exact results is then

$$\stackrel{~}{u}_2=\frac{1}{2}\underset{i=1}{\overset{N}{}}e_i^2.$$
(14)
In our conventions the hyperelliptic curve associated with the matching prescription for $`SU(N)`$ theories with $`2N`$ massless fundamental matter hypermultiplets is :

$$y^2=\left(x^N\underset{k=1}{\overset{N1}{}}\stackrel{~}{u}_{k+1}x^{Nk1}\right)^2\stackrel{~}{q}x^{2N}.$$
(15)
The functions $`\{a_i\}`$ and $`\{a_{D,i}\}`$ can be determined by evaluating the meromorphic one-form

$$\lambda =\frac{xdx}{2\pi iy}\left[\frac{F(x)G^{}(x)}{2G(x)}F^{}(x)\right],$$
(16)
where the prime denotes differentiation with respect to $`x`$, over the canonical basis of homology one-cycles $`\{\alpha _i,\beta _i\}`$ of the Riemann surface:

$`a_i`$ $`=`$ $`{\displaystyle _{\alpha _i}}\lambda ,`$ (17)
$`a_{D,i}`$ $`=`$ $`{\displaystyle _{\beta _i}}\lambda .`$ (18)
Given the curve defining the moduli space of the theory, one can then determine the VEVs $`\{u_i\}`$, and hence the prepotential $``$, exactly via Eqs. (7, 16, 17, 18).
One can perform the integration of the meromorphic one-form in Eq. (16) exactly for $`SU(N)`$ theories . Using the curve in Eq. (15), to order $`𝒪(\stackrel{~}{q})`$ we have:

$$a_i=_{\alpha _i}\frac{dx}{2\pi i}\left(N\frac{xF^{}(x)}{F(x)}+\stackrel{~}{q}\frac{x^{2N}(NFxF^{})}{2F(x)^3}+𝒪(\stackrel{~}{q}^2)\right).$$
(19)
At weak coupling, the homology one-cycles $`\{\alpha _i\}`$ coincide with the exact VEVs $`\{e_i\}`$. (We note that the electric-magnetic duality ambiguity  in the case considered here is trivial.) Discarding a total derivative and performing the integration in Eq. (19) yields the $`(N1)`$ equations

$$a_i=e_i+\stackrel{~}{q}\frac{e_i^{2N1}}{2\mathrm{\Delta }_i(e_i)}\left(N\underset{ij}{}\frac{e_i}{(e_ie_j)}\right),$$
(20)
where $`\mathrm{\Delta }(e_i)=_{il}(e_ie_l)`$, $`l=1\mathrm{}N1`$. In obtaining this result, the reverse direction was taken in performing the period integral Eq. (17), so that $`a_i`$ remains positive. Solving Eq. (20) at leading order in $`\stackrel{~}{q}1`$ gives

$$e_i=a_i\stackrel{~}{q}\frac{a_i^{2N1}}{2\mathrm{\Delta }_i(a_i)}\left(N\underset{ij}{}\frac{a_i}{(a_ia_j)}\right).$$
(21)
Equations (10, 13, 21) show that at the classical level the expected results are reproduced.
We now write $`\mathrm{\Delta }_i\mathrm{\Delta }(a_i)`$ for simplicity. Using the definition in Eq. (14), the exact result $`\stackrel{~}{u}_2`$ is given by

$$\stackrel{~}{u}_2=u_2^{\mathrm{cl}}+\stackrel{~}{q}\stackrel{~}{u}_2^{\mathrm{inst}}=\frac{1}{2}\underset{i=1}{\overset{N}{}}a_i^2\stackrel{~}{q}\underset{i=1}{\overset{N}{}}\frac{a_i^{2N}}{2\mathrm{\Delta }_i^2}\left(1\underset{ij}{}\frac{a_j}{(a_ia_j)}\right).$$
(22)
Equation (22) is the main result of this letter. Our expression for the quantum modulus $`\stackrel{~}{u}_2`$ holds for general values of $`N`$ and agrees up to regular terms with previous results . The Matone relation  enables one to relate $`\stackrel{~}{u}_2`$ to the one-instanton prepotential $`_1`$ derived from the exact results, and we shall employ this in the test of the scheme. The scheme  for the matching of the exact results and the instanton results accounts for the most general mapping which can connect the parameters and the VEVs for either set of results. The parameters and VEVs are $`\{\stackrel{~}{q},\stackrel{~}{u_i}\}`$ for the exact results, and $`\{q,u_i\}`$ for the instanton results, since the hypermultiplet masses are set to zero here. The matching prescription to $`𝒪(q)`$ for the coupling parameter $`q`$ and the VEV $`u_2`$ is:

$`\stackrel{~}{q}`$ $`=`$ $`C_0q,`$ (23)
$`\stackrel{~}{u}_2`$ $`=`$ $`(1+A_1^{(1;1)}q)u_2.`$ (24)
In this matching, no modular invariance (S-duality) is assumed for the space of couplings in the reparameterisation in Eqs. (2, 23), as it is not necessarily a physical attribute of the theory. Writing Eqs. (23, 24) in terms of perturbative and non-perturbative parts gives

$$C_0\stackrel{~}{u}_2^{\mathrm{inst}}=u_2^{\mathrm{inst}}+A_1^{(1;1)}u_2^{\mathrm{cl}}.$$
(25)
The classical VEV $`u_2^{\mathrm{cl}}`$ in the matching above (Eq. (25)) constitutes a regular term; the other terms will then have the same singularity structure .
The Matone relation  (see  for instanton based derivations of this relation) for the quantities found using instanton calculus is

$$2\pi i_n=u_2^n,$$
(26)
and the result for the one-instanton prepotential $`_1`$ for $`𝒩=2`$ supersymmetric $`SU(N)`$ gauge theories with $`2N`$ massless fundamental matter hypermultiplets  is

$$_1=\frac{iC_1^{}\pi ^{2N1}}{2^{2N+2}}\underset{ij}{}\frac{(a_i+a_j)^{2N}}{(a_ia_j)^2\gamma _i\gamma _j},$$
(27)
where $`\mathrm{\Delta }_i=(a_ia_j)\gamma _i`$ and $`\gamma _i=_{ik,kj}(a_ia_k)`$, $`k=1\mathrm{}N2`$. Inserting Eqs. (26, 27) and the exact result for $`\stackrel{~}{u}_2`$, from Eq. (22), into Eq. (25), one has

$$C_0\underset{i=1}{\overset{N}{}}\frac{a_i^{2N}}{2\mathrm{\Delta }_i^2}\left(1\underset{ij}{}\frac{a_j}{a_ia_j}\right)=\frac{C_1^{}\pi ^{2N}}{2^{2N+1}}\underset{i=1}{\overset{N}{}}\underset{ij}{}\frac{(a_i+a_j)^{2N}}{(a_ia_j)^2\gamma _i\gamma _j}+\frac{1}{2}A_1^{(1;1)}\underset{i=1}{\overset{N}{}}a_i^2.$$
(28)
The constant $`C_1^{}`$ is the renormalisation scheme-dependent one-instanton factor; we follow  and use

$$C_1^{}=2^{N+2}\pi ^{2N}.$$
(29)
To extract the constant $`C_0`$, one observes that manipulating Eq. (28) so that both sides of the equality have the same denominator enables one to take the previously singular limit $`a_ia_j`$ and to compare the coefficients of the non-vanishing terms. This gives

$$C_0=2^{N+2}.$$
(30)
The form of $`C_0`$ found here (Eq. (30)) is in exact agreement with the form of $`C_0`$ determined in  for a single point of the moduli space in $`𝒩=2`$ supersymmetric $`SU(N)`$ gauge theories. Using the expression for $`C_0`$, one can explicitly determine the constant $`A^{(1;1)}`$ by expanding Eq. (28) and comparing the coefficients of the leading order terms. Expanding the following objects as

$`\mathrm{\Delta }_i`$ $``$ $`a_i^{N1}+\mathrm{},`$ (31)
$`\gamma _i`$ $``$ $`a_i^{N2}+\mathrm{},`$ (32)
$`(a_i+a_j)^{2N}`$ $`=`$ $`{\displaystyle \underset{r=0}{\overset{2N}{}}}\left(\begin{array}{c}2N\\ r\end{array}\right)a_i^{2Nr}a_j^r,`$ (35)
and comparing the coefficients of the terms of highest order in $`a_i`$ in Eq. (28) implies

$$A_1^{(1;1)}=2^{N+2}+2^{2N}\left(\begin{array}{c}2N\\ N1\end{array}\right).$$
(36)
The formulae Eqs. (30, 36) are valid for general values of $`N`$; this has been checked by an inductive argument using Eq. (28), which we omit here.
Given that the form of the matching between the exact results and the instanton predictions is correct at the one-instanton level, as the preceding calculation and subsequent agreement shows, we now comment on its relevance to the string theoretic derivation of the class of hyperelliptic curves corresponding to those found in the exact solutions . The classical brane configuration of $`N`$ D4-branes suspended between two parallel NS5-branes in Type-IIA string theory corresponds to the vacua of classical $`𝒩=2`$ supersymmetric Yang-Mills theory, and the vacua of the quantum theory corresponds to the supersymmetric configurations of an $`M`$-theory $`M`$5-brane with a particular world volume . To incorporate $`N_f`$ matter hypermultiplets into the system, one attaches $`N_f`$ semi-infinite 4-branes to the NS5-branes. The class of curves corresponding to those appearing in the exact results follows from this brane configuration. For detailed reviews of this construction, see, for example, references .
The dictionary  set up between the parameters of the brane configuration and the parameters of the field theory is only valid at extremely weak coupling. Beyond extremely weak coupling, quantum corrections will in general modify this dictionary, and it will contain ambiguities manifest as the freedom to make non-perturbative redefinitions of the parameters. Hence, the brane–field theory correspondence is valid, but the quantitative dictionary connecting them is ambiguous.
The matching scheme in  uses the most general permissible redefinitions of the parameters and VEVs of these field theories. It is natural to propose that the ambiguities in the $`M`$-theoretic derivation of the curves which exactly solve the low-energy effective actions of the same field theories are resolved by the same matching scheme. That is, the equivalence class of curves derived from $`M`$-theory should coincide precisely with the equivalence class of curves derived from the exact solutions. Then the mappings between the elements of the equivalence class of $`M`$-theoretic curves will be the same mappings between elements of the equivalence class of exact solution curves. Hence the same matching scheme for comparing the exact results to instanton results should also hold for comparing the $`M`$-theoretic results to instanton results. A precise test of this conjecture would validate the matching scheme  for future $`M`$-theory predictions.

In conclusion, the exact results and the instanton predictions can be matched to one-instanton level for $`𝒩=2`$ supersymmetric gauge theories with gauge group $`SU(N)`$ and $`2N`$ massless fundamental matter hypermultiplets, for general $`N`$, via the scheme described, with constants $`C_0`$ and $`A_1^{(1;1)}`$ given by Eqs. (30, 36). It has been shown that this matching can be achieved for the complete quantum moduli space of these theories at the one-instanton level, extending the previous result  for a moduli subspace. The case where the $`2N`$ hypermultiplets have non-zero masses could also be investigated; one expects that the constants $`C_0`$ and $`A_1^{(1;1)}`$ should be the same as in the massless case due to renormalisation group flow arguments. Two-instanton level checks of the matching between the exact results and the instanton predictions in these theories have been performed  only for the case $`N=2`$, where agreement was found. Tests of the scheme at the two-instanton level for general $`N`$ would be desirable since these would provide a physical check of the matching beyond the above constraints on the values of the constants involved.

Acknowledgement

The author is indebted to V. V. Khoze for advice, numerous discussions and for introducing him to this problem. The author acknowledges receipt of a PPARC studentship.


format: text

text:
# Magnetic Phase Diagram of Ca2-xSrxRuO4 Governed by Structural Distortions

\[

## Abstract

We constructed, by the first-principles calculations, a magnetic phase diagram of Sr2RuO4 in the space spanned by structural distortions. Our phase diagram can qualitatively explain the experimental one for Ca2-xSrxRuO4. We found that the rotation and the tilting of RuO6 octahedron are responsible for the ferro- and antiferro-magnetism, respectively, while the flattening of RuO6 is the key factor to stabilize those magnetic ground states. Our results imply that the magnetic and the structural instabilities in Sr2RuO4 are closely correlated cooperatively rather than competitively.

\]

Both the magnetic and the structural instabilities are essential issues for the unconventional superconductivity in Sr2RuO4 , which is the only example of a noncuprate layered perovskite superconductor. It was first suggested that the Sr2RuO4 is close to the ferromagnetic (FM) instability  with strong FM spin fluctuations, which may naturally lead to a spin-triplet $`p`$-wave pairing mechanism . However, the recent observation  of sizable antiferromagnetic (AF) incommensurate spin fluctuation, due to the Fermi surface nesting , indicates that more careful studies are needed. As for the structural aspect, it was pointed out by experiments that Sr2RuO4 is very close to the structural instability with respect to the RuO6 rotation . With such a situation, one may consider that three kinds of instabilities, superconducting, magnetic and structural ones, may compete. Nevertheless, the correlation among those instabilities has not been fully discussed. It was found recently  that the cleaved surface of this material is reconstructed to form the c(2x2) structure which can be regarded as the frozen RuO6 rotation mentioned above. Furthermore, the density-functional calculation predicts that the surface ferromagnetism is strongly stabilized by the structural reconstruction . This prediction suggests that the structural and magnetic instabilities cooperate rather than compete, although the surface ferromagnetism has not been experimentally confirmed up to now.

On the other hand, the recent studies on the Ca2-xSrxRuO4 suggest strongly the cooperative feature of the structural and magnetic instabilities in the bulk. Moreover, the system shows a very rich phase diagram and provides us with an opportunity to analyze the correlation between the magnetism and the structure more extensively. Below is a brief description of the experimental observation for Ca2-xSrxRuO4 by Nakatsuji et al. With the Ca substitution for Sr, the system is successively driven from the non-magnetic (NM) 2-dimensional (2D) Fermi liquid ($`x2.0`$) to a nearly FM metal ($`x0.5`$), an antiferromagnetically correlated metal ($`0.2<x<0.5`$), and finally an AF insulator ($`x<0.2`$). Since the substitution is isovalent, the dominant effects are the structural modifications due to the reduced ionic size of Ca compared with Sr. Evidence  has been presented by neutron scattering that the structural distortions caused by the Ca substitution correlate with the changes in the magnetic and the electronic properties.

The main aim of this letter is to study how and why the magnetism of Ca2-xSrxRuO4 is affected by structural distortions. In order to extract essential aspects, we assume that for a given crystal structure, the electronic structure is not affected by the relative content of Ca and Sr. Therefore, in the following, we study the stable magnetic phases of Sr2RuO4 for given structural distortions. Three types of structure distortions, i.e. RuO6 octahedron rotation about the $`c`$-axis, RuO6 tilting around an axis parallel to the edge of octahedron basal plane and the flattening of RuO6 along the $`c`$-axis are identified from experiments . Our phase diagram can qualitatively explain the experimental phase diagram of Ca2-xSrxRuO4, demonstrating the crucial roles of structural distortions for the tuning of electronic and magnetic properties, and further supporting our previous prediction for the surface . In particular, we found that the RuO6 rotation can enhance the FM instability significantly, while the combination of tilting and rotation of RuO6 is responsible for the enhancement of AF instability. Furthermore, we point out that the flattening of RuO6 is a key factor to stabilize the magnetic (both FM and AF) ground states. The basic physics governing the phase diagram can be understood in terms of the strong coupling between the lattice and the magnetism through the orbital degrees of freedom. Our results strongly suggest that, in Sr2RuO4, the magnetic fluctuations can be significantly enhanced by the structural fluctuations, implying the necessity of reconsidering the coupling mechanism in the bulk superconductivity.

The calculations were performed with the first-principles plane-wave basis pseudopotential method based on the local density approximation (LDA). The validity of LDA treatment for ruthenates was demonstrated in Ref. . The $`2p`$ states of oxygen and $`4d`$ states of Ru are treated by the Vanderbilt ultrasoft pseudopotential , while the norm-conserving scheme  is used for other states. The cutoff energy for the wave function expansion is 30 Ry. The $`k`$-point sampling of the Brillouin zone (BZ) was well checked to provide enough precision in the calculated total energies. The theoretically optimized lattice parameters $`a=3.84`$ Å and $`c=12.70`$ Å for the bulk Sr2RuO4 are in good agreement with the experimental ones $`a=3.86`$ Å and $`c=12.73`$ Å. The degree of flattening of RuO6 octahedron $`\lambda `$ is defined by $`\lambda =d_c/d_{ab}`$ with $`d_c`$ ($`d_{ab}`$) denoting the Ru-O bond length along the $`c`$-axis (in the $`ab`$-plane) with the RuO6 volume fixed. Rotation and tilting of RuO6 octahedron are operated with the Ru-O bond lengths fixed. In order to construct the magnetic phase diagram, the lowest energy magnetic phase for each crystal structure is searched for among different (NM, FM and AF) phases. In the present work, we focus our attention only on phases described within the c(2x2) unit cell. The soft phonon mode of $`\mathrm{\Sigma }_3`$ at the zone boundary in Sr2RuO4  and the AF state of Ca2RuO4  are in this category.

Figure 1 shows the calculated phase diagram , while the Table 1 summarizes the calculated total energies and magnetic moments for some particular points in the phase diagram. Hereafter, $`\varphi `$ and $`\theta `$ denote the rotation angle and the tilting angle, respectively. The apical oxygen and the oxygen in the $`ab`$-plane are called O(2) and O(1) respectively. From right to left of the phase diagram, first the RuO6 starts to rotate along the $`c`$-axis by up to 12, and then with the 12 rotation being fixed, the RuO6 starts to tilt up to 12. The structural analysis by the neutron scattering  allows us to make a one-to-one correspondence between the structural changes, i.e. the horizontal axis of our phase diagram, and the doping level $`x`$ in Ca2-xSrxRuO4. For $`x=2.0`$ (Sr2RuO4), the system has $`I4/mmm`$ symmetry with $`\varphi =\theta =0^{}`$, corresponding to the right end of our phase diagram. With reduction of $`x`$, RuO6 starts to rotate and the symmetry is reduced to $`I4_1/acd`$ until $`x=0.5`$ (Ca1.5Sr0.5RuO4), where $`\varphi =12.78^{}`$ and $`\theta =0^{}`$ at 10 K. With further reduction of $`x`$, RuO6 starts to tilt and the symmetry is further reduced to $`Pbca`$ until $`x=0.0`$ (Ca2RuO4) where $`\varphi =11.93^{}`$ and $`\theta 12^{}`$ at low temperature, corresponding to the left end of our phase diagram. It was also pointed out by the experiment  that, from $`x=2.0`$ to $`x=0.5`$, the degree of flattening $`\lambda `$ remains almost constant ($`1.07`$), while from $`x=0.5`$ to $`x=0.0`$, the rotation angle $`\varphi `$ is almost unchanged ($`12^{}`$). Three representative experimental points are shown in our phase diagram by triangles. Now, the basic tendency suggested by our phase diagram is that the RuO6 rotation will drive the system from a NM state to a FM state, while the subsequent tilting plus the flattening of RuO6 will push the system to an AF region. This general tendency is quite consistent with the experimental results. It should be noted that the rich experimental phase diagram can be simply understood in terms of the close coupling between structural distortions and magnetism. Although the real long range FM ordering in Ca1.5Sr0.5RuO4 is not confirmed yet, the significant enhancement of spin susceptibility in this doping level undoubtedly implies the strengthening of FM instability. Another important aspect in our phase diagram is that the flattening of RuO6 is so important not only for the AF state but also for the FM state. This suggests that simply by uniaxial pressure, the Sr2RuO4 can be driven from the NM state to a FM state.

The basic questions concerning our phase diagram are: 1) Why are the RuO6 rotation and tilting correlated with the tendency to the FM and AF states? 2) Why is the RuO6 flattening so important for the magnetic solutions? Before answering these questions, let us discuss the role of each $`4d`$ orbital in the electronic properties of Sr2RuO4, which is essential to our later discussions. The three Ru $`t_{2g}`$ orbitals ($`d_{xy}`$, $`d_{yz}`$, $`d_{zx}`$) hybridize with each other only very weakly in tetragonal Sr2RuO4. Therefore each orbital plays distinct roles. The projected density of states (DOS) shown in Fig.2(a) indicates that the $`d_{xy}`$ orbital contributes dominantly to the well-known van Hove singularity (VHS) just above the Fermi level. The $`\gamma `$ Fermi surface has the character of $`d_{xy}`$. It is mostly responsible to FM spin fluctuation due to the high DOS around the Fermi level. On the other hand, $`d_{yz}`$ and $`d_{zx}`$ orbitals contribute to the $`\alpha `$ and $`\beta `$ Fermi surfaces and produce the incommensurate spin fluctuation coming from the strong nesting effect due to the quasi-one-dimensional nature of those states. The calculated bare spin susceptibility shown in Fig. 3(a) for the undistorted compound has the incommensurate peak at Q=($`2\pi /3a,2\pi /3a`$), being consistent with previous calculations . In real materials, the two factors, i.e. the FM instability due to the high DOS at the Fermi level combined with the $`q`$ dependent Stoner factor  and the AF instability due to the nesting effect , will compete.

The RuO6 rotation couples mostly with the $`d_{xy}`$ orbital but not with the $`d_{yz}`$, $`d_{zx}`$ orbitals because the $`pd\pi `$ type hybridization between the O(1)-$`2p`$ and the $`d_{xy}`$ states will be significantly reduced by the RuO6 rotation, but those between the O-$`2p`$ and the $`d_{yz}`$, $`d_{zx}`$ states are not affected so much. The direct results of this reduced $`pd\pi `$ type hybridization between the O(1)-$`2p`$ and the $`d_{xy}`$ states are, first the narrowing of $`d_{xy}`$ band width and second the downward shift of $`d_{xy}`$ band, as shown in Fig. 2(b) (about 0.4 eV narrowing and 0.1 eV downward shift of $`d_{xy}`$ band for $`\varphi =12^{},\lambda =1.07`$). As the latter brings the VHS closer to the Fermi level, both of the two results will enhance the DOS at the Fermi level. Another effect coming from the downward shift of the $`d_{xy}`$ states is the population reduction in the $`d_{yz}`$, $`d_{zx}`$ states, which may shift the Fermi surface nesting vector closer to the zone boundary. However, the increase of the DOS at the Fermi level is the dominant effect and the tendency towards FM instability is enhanced by the RuO6 rotation. Once tilting is additionally introduced, all of the $`t_{2g}`$ bands will become narrower. This will enhance the nesting effect and enhance the AF instability. The discussion so far can answer the first question.

Now let us discuss the effects of RuO6 flattening. There are two factors also. First, with the flattening of RuO6 octahedron, the Ru-O(1) bond length will increase, while the Ru-O(1)-Ru angle remains $`180^{}`$. The increased bond length will reduce all the $`pd\pi `$ type hybridizations between the O(1)-$`2p`$ and the $`d_{xy}`$, $`d_{yz}`$, $`d_{zx}`$ states. Therefore, width of all the bands of three Ru-$`t_{2g}`$ states is reduced as shown in Fig.2(c) (about 0.4 eV for $`d_{yz,zx}`$ bands and 0.3 eV for $`d_{xy}`$ band for $`\lambda =0.96`$), making the DOS at the Fermi level higher. This will favor the FM solution. Another very important results of flattening is the orbital polarization. It is obvious that the tetragonal distortion by the flattening will populate the $`d_{xy}`$ states and depopulate the $`d_{yz}`$, $`d_{zx}`$ states (about 0.2 eV downward shift of $`d_{xy}`$ band for $`\lambda =0.96`$). The effect is similar to the RuO6 rotation. The orbital polarization due to flattening will also shift the nesting vector to the zone boundary as shown in the susceptibility calculations (Fig.3). This will favor the commensurate AF state of the system. The net effect by rotation, tilting and flattening of the RuO6 will depend on the competition among them and the phase diagram of Fig.1 demonstrates the situation in a space spanned by those distortion modes. Figure 2(f) shows the calculated DOS for the AF state with $`\varphi =\theta =12^{}`$, i.e. almost the experimental structure of Ca2RuO4. It is clear in this case that the occupied minority spin states mostly come from the $`d_{xy}`$ orbital due to flattening of RuO6. Therefore, the strong superexchange interaction between the occupied majority-spin and unoccupied minority-spin $`d_{yz}`$, $`d_{zx}`$ orbitals will stabilize the AF ground state .

In summary, by constructing a phase diagram of Sr2RuO4 with structural distortions, we find the strong coupling between the lattice and the magnetism. Our phase diagram can qualitatively explain the experimental phase diagram of Ca2-xSrxRuO4. We demonstrate that the RuO6 rotation will enhance the FM instability in the system, while the tilting plus the flattening of RuO6 make the system AF. We pointed out that the flattening of RuO6 is so important not only for the AF state but also for the FM state. An important implication of our results is that the magnetic and the structural instabilities in Sr2RuO4 should be strongly correlated. The structure fluctuation and the magnetic fluctuation cooperate. Actually the phonon mode corresponding to the RuO6 rotation is quite soft  in the bulk, and this rotation will enhance the FM instability. All these results imply the necessity of reconsidering the coupling mechanism for the unconventional superconductivity. In this context, we propose a possible way to identify experimentally the relationship between the FM fluctuation and the superconducting state. As the uniaxial compression of Sr2RuO4 will enhance the FM fluctuation without introducing the disorder, the variation of superconducting transition temperature against uniaxial compression may provide important information.

The authors thank Dr. S. Nakatsuji, Prof. E. W. Plummer, Prof. Y. Tokura, Dr. R. Matzdorf, and Dr. J. D. Zhang for valuable discussions and for providing us with their experimental data. One of the anthors (Z. Fang) also thanks Dr. David Singh for valuable discussions. The present work is partly supported by NEDO.


text:
# Localization and diffusion in Ising-type quantum networks

## I Introduction

Maturing of the techniques for controlling single quantum systems such as atoms, photons, or quantum dots has made it possible to design more complicated structures of these basic elements: quantum networks. Interesting examples of such micro- and nanostructures are for instance thin wires used as atom guides , trapped chains of atoms and ions , or structures grown on a substrate such as quantum dots  or arrays of nanomagnets. These are among the first steps of laboratory demonstrations leading to quantum network engineering and finally to matured and scalable quantum technology. Due to these possible applications investigations of quantum networks have become topical.

In this paper we consider a quantum network configuration which is perhaps the simplest non-trivial one: a network which realizes nearest-neighbour interactions by coupling the individual quantum systems pairwise. As a possible implementation one can think of light modes coupled together by an optical element, for example a beam splitter or a mode coupler. The network is schematically presented in Fig.1. Another common form of nearest neighbour interaction is simulataneous interaction with many neighbours, like electrons interacting with Coulomb interaction. Here we, however, consider pairwise interaction which is typical for light modes, and for many quantum networks in general.

The network in Fig.1 has an interesting and potentially useful connection to the Ising-model . Interference effects are typical in this kind of networks and have been investigated by many authors . The Hamiltonian structure of the network was studied in . In this paper we consider the effect of randomness in such a network. Randomness in certain quantum systems is know to lead to intriguing effects such as Anderson localization and dynamical localization . Also classical wave behaviour combined with randomness can lead to localization, c.f. localization of light

We can view the considered network also as a interferometer. The effects demonstrated can be understood as the result of interferences between paths available for the excitations in the network. However, interferences are known to be very fragile towards phase perturbations. We show that the typical interference effects are modified or completely destroyed as randomness is introduced in the network.

The way in which the randomness enters the network has a profound effect on the behaviour. Therefore we consider two basic variants of the network. First, a 2-D network of the type shown in Fig.1, with randomness all over the network. This means that along each connecting line between two interaction nodes, the phase (free time evolution) is changed randomly within the possible range of $`0,2\pi `$. A corresponding physical system could be for example a network of optical elements with phase fluctuations arising from different optical lengths between the elements. The fluctuations could be produced at will with phase shifters to make systematic studies of the effect of randomness. Alternatively, our calculations can model undesired phase fluctuations in optical networks and can be used to study the stability of interference effects in a complicated interferometer. We show that in the random 2-D network classical stochastic diffusion is typical.

In the second configuration considered, the network models a 1-D system. Now the basic unit of the network contains fixed randomness. The basic unit is formed by the two columns of $`A`$ and $`B`$ (”the motif” element for building the network) in Fig.1, and along the connecting lines between the columns the phase is adjusted arbitrarily, i.e., within the range of $`0,2\pi `$. In such a way we define a random motif matrix. Then this very same basic unit is precisely repeated $`M`$ times. The horizontal dimension of the network can now be understood as time and the network actually models time evolution in one dimension (the vertical dimension is the real space or a mode space). Note that the configuration also corresponds to a thin layer of material inside a cavity, with a light pulse crossing it several times: the 1-D propagation is now the spreading of the light along the material layer after many cavity crossings. We show that randomness in the 1-D type configuration leads to quantum localization.

We explain the different behaviour of the two configurations by relating our observations to the standard theory of localization. The theory of Anderson localization and dynamical localization is formulated considering the properties of the Hamiltonian for the system of interest. This is not always the best possible approach in the quantum engineering context, where the individual components of a more complicated network are described by a unitary transform. The transfer matrix for the total network is easily obtained by simple multiplications of the transfer matrices corresponding to the individual elements, but to find out the Hamiltonian corresponding to the complex system is non-trivial if not impossible . Therefore we find it of interest to consider here, for the first time to our knowledge, the quantum localization problem by starting the theoretical description from the transfer matrices of the system.

Note that in the mathematical level, the unitary transforms describing the network can also be understood as operations performed in the computational space of a quantum computer. Thus our investigations also relate to decoherence affecting certain types of computational operations. Nearest-neighbour coupling configurations have been considered in the quantum information context recently .

In section 2 we introduce the transfer matrix which is the common starting point for both network configurations of interest. The first configuration, leading to classical diffusion, is considered in section 3, the second, showing localization, in section 4. Differences in the quantum behaviour in these two cases are explained with the theory of Anderson/dynamical localization in section 5. The conclusions are presented in section 6.

## II The transfer matrix

The network we consider realizes the nearest-neighbour structure given in Fig.1. Each element (box in the figure) represents a $`2\times 2`$ coupler, e.g. a beam-splitter in case of optical networks. The whole network is build up by repeating the motif structure consisting of a row of components A and B, i.e., we use the outputs of the motif as the input for the next layer. We assume periodic boundary conditions .

The transfer matrix describing the motif unit has the form

$`U=\left[\begin{array}{cccc}\text{A}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{c}\mathrm{}\\ \mathrm{}\end{array}& \\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \text{A}& & \\ \begin{array}{cc}\mathrm{}& \mathrm{}\end{array}& & \mathrm{}& \\ & & & \text{A}\end{array}\right]\left[\begin{array}{ccccc}B_{22}& \begin{array}{cc}0& 0\end{array}& \mathrm{}& & B_{21}\\ \begin{array}{c}0\\ 0\end{array}& \text{B}& & & \\ \mathrm{}& & \text{B}& & \\ & & & \mathrm{}& \\ B_{12}& & & & B_{11}\end{array}\right],`$ (20)
where

$`\text{A}=\left[\begin{array}{cc}\mathrm{cos}\theta & \mathrm{sin}\theta \\ \mathrm{sin}\theta & \mathrm{cos}\theta \end{array}\right]`$ ; $`\text{B}=\left[\begin{array}{cc}\mathrm{cos}\varphi & \mathrm{sin}\varphi \\ \mathrm{sin}\varphi & \mathrm{cos}\varphi \end{array}\right].`$ (25)
When the parameters $`\theta `$ and $`\varphi `$ are real, (25) describes a SU(2) network. Purely imaginary $`\theta `$ and $`\varphi `$ correspond to a SU(1,1) network. Above, $`\mathrm{sin}^2x`$ ($`x`$ is $`\theta `$ or $`\varphi `$) gives the coupling strengths of the elementary units used and $`N`$ is the number of the units in one layer. We will consider only the case $`\mathrm{cos}^2x=0.5`$. In optics, this corresponds to so called balanced beam splitters (or 3-dB couplers).

The transfer matrix $`W`$ of the whole network, in the case all motif sturctures are the same, is given as a power of $`U`$:

$$W(M)=U^M.$$
(26)
The other possibility is that all the motif structures $`U_i`$, $`i[1,M]`$ are different due to overall randomness, then

$$W(M)=U_1D_1U_2D_2\mathrm{}U_M.$$
(27)
The matrices $`D_i`$ represent diagonal matrices giving the phase uncertainities between the random motifs. Here $`M`$ can be understood as the size ($`2N\times M`$) of the network in the case of a 2-D network, or as time when the network models 1-D propagation.

An example of the typical interference behaviour in complete absence of randomness is presented in Fig.2

In the following we will show how this interference pattern is modified by different forms of randomness.

## III Random 2-D network

We consider first a two-dimensional quantum network with randomness. This corresponds to the total transfer matrix

$$W(M)=U_1D_1U_2D_2\mathrm{}U_M,$$
(28)
where $`U_i,D_i`$ are all different. The randomness enters the motif stuctures either via the individual quantum components $`A`$ and $`B`$, or by having an extra random phase shift matrix between the matrices containing the $`N`$ components $`A`$ and $`B`$. For optical systems, the first case would correspond to having fluctuations in the transmittivities of the beam splitters, the latter to random phase shifts occuring in the optical paths. Here we have analyzed phase randomness. The phases between all the connecting points change randomly within the $`0,2\pi `$ interval.

We study the behaviour of the network by following the propagation of an input state initially located at the input port which is at the center of the network. To illustrate our results we fixed the number of network input ports to $`40`$ which corresponds to $`20`$ elements in one layer. As the initial state we consider a single photon guided into the network from the $`20`$-th input port

$$|\psi (0)=|0_1\mathrm{}|0_{19}|1_{20}\mathrm{}|0_{40}.$$
(29)
The spread of a photon given by this initial state in the network is shown in Fig.3

The plot shows the logarithm of the probability for the photon to appear at the given output (given by the dots). Three choices for $`M`$ are shown to illustrate how the distribution changes when the photon propagates in the network. The largest value of $`M`$ equals ten; after ten layers the photon hits the border of the network, then traversing just another motif layer the photon would be recycled back into the system due to periodic boundary conditions. In calculating the probability we have averaged over $`1000`$ runs. The solid line gives the least square fit to the calculated values to show that the distribution follows a quadratic dependence.

As shown in Fig.3, the quantum state in this network behaves like a classical distribution, i.e. it spreads into a broad Gaussian at the outputs. Classical and quantum behaviour in this kind of networks, in the absence of randomness, is studied in . Gaussian outputs were found in the classical case  whereas complicated interference and recurrence phenomena are typical in the quantum case. Our results show that introducing randomness acts as decoherence which reduces the quantum behaviour to the classical one. In section 5 we explain the results further in the light of random matrix theory of localization.

## IV Propagation in 1-D; a thin layer of material

Next we consider the situation where the transfer matrix of one motif structure of the network has the form $`UD`$, where $`D`$ is a diagonal matrix with $`\mathrm{𝑓𝑖𝑥𝑒𝑑}`$ randomness, and the total transfer matrix is simply

$$W(M)=(UD)^M.$$
(30)
Physically this corresponds to quantum propagation in one dimension, or for example a thin layer of material inside a cavity through which light makes several passages. The situation is illustrated in Fig.4

We consider the behaviour of the same input as given in the previous section. Fig.4 shows that the distribution becomes exponentialy peaked at the input position $`20`$. To see the emergence of this exponential localization we present the evolution of the initial photon distribution in dependence on the number of motifs traversed ($`M`$).

We have found that the complicated interference pattern typical for the ideal network turns gradually into quantum localization when the amount of randomness is increased. This is demonstrated in Fig.5.

We compare the distributions at $`M=10`$ for weak and strong randomness. The strength of randomness is controlled by the extend to which the phases in the random diagonal matrices can change, for instance for weak randomness only within the interval $`0,\pi /16`$ and for moderate randomness the interval would be larger $`0,\pi /4`$. The main result is that we observe exponential localization for strong enough randomness. This remarkable quantum feature will be explained in section 5 by making a connection to the standard theory of Anderson and dynamical localization.

### A Intermediate case: extra noise in 1-D propagation or semi-random 2-D network

As an intermediate case we consider a network with the transfer matrix

$$W(M)=UDD_1UDD_2UDD_3\mathrm{},$$
(31)
where $`D_i`$ are random and different from each other. The matrix $`D`$ is also random. Note that with all $`D_i=1`$, $`W(M)`$ would be the same as in the previous section, that is, expected to lead to localization. This intermediate case models a network with phase fluctuations inbetween exactly similar motif structures. Examples of this are an optical network where the nearest neighbour coupling could be essentially perfect but the path lengths inbetween two nearest neighbour coupling layers fluctuate, or fluctuating cavity length in the 1-D case of a layer of material in a cavity. In modelling 1-D propagation, this is also a test of numerical stability of the localization phenomenon when very small extra randomness is introduced.

We now fix the amount of randomness in the matrix $`D`$ that appears repeatedly in $`W(M)`$ (phases chosen between $`0,2\pi `$), and vary the amount of randomness in $`D_i`$. We observe that for weak randomness in $`D_i`$ (phases chosen between $`0,0.1\pi `$), the localization is preserved (note that in this case we probably have approximately the situation $`W(UD)^M`$). For strong randomness (phases chosen between $`0,2\pi `$) Gaussian diffusion is observed. This transition from localization into diffusion is illustrated in Fig.6.

## V Localization — Hamiltonian versus transfer matrix approach

The theory of Anderson and dynamical localization tells that exponential localization of the system eigenvectors, and consequently, of quantum wave packets, occurs when the system Hamiltonian is a random band-diagonal matrix . Band-diagonality means that only a few neighbouring sites (thinking in terms of a lattice) are coupled — in the Lloyd model , which is the simplest model showing localization, only the nearest neighbours are connected. Randomness usually means that the energies of the sites, that is, the diagonal terms of the Hamiltonian, are randomly chosen according to a certain distribution. But randomness can equally well appear in the off-diagonal couplings.

We describe the network of interest by a transfer matrix rather than a Hamiltonian. Connection between these two is non-trivial to obtain for large networks, c.f . But some overall characteristics are easy to sketch. One motif structure of our network clearly couples nearest neighbours, which corresponds to a Hamiltonian $`H_{motif}`$ with only nearest neighbour coupling. In the case of 1-D propagation, or a thin layer of material in a cavity, the total transfer matrix is

$$W(M)=U^M=e^{iH_{motif}M}.$$
(32)
The network size $`M`$ can be interpreted as time, and the Hamiltonian of the system is the one describing one motif structure: random and band diagonal. Therefore it is in accordance with the theory of localization to observe an exponentially localized distribution at the output. The case with overall randomness in a 2-D network is completely different. The total transfer matrix

$$W(M)=U_1D_1U_2D_2\mathrm{}U_M$$
(33)
can no longer be expressed as an $`M`$th power of a certain band-diagonal Hamiltonian. By diagonalizing $`W(M)`$ one can obtain the corresponding Hamiltonian, but it is clear that it will not any more be band-diagonal, but instead couples $`M`$ neighbours. For large $`M`$, the Hamiltonian approaches a full random matrix. It is known that for full random matrices there is no localization . Instead, we observe classical-type diffusion.

## VI Conclusions

We have considered randomness and localization in a specific quantum network which is characterized by nearest-neighbour structure and has a relation to models of statistical mechanics with nearest-neighbour interaction. The physical systems that can be modelled by such a network are for instance beam splitter networks , micro- and nanostructures and networks of intersecting energy levels . We showed that in a 2-D network with overall randomness, the quantum behaviour was reduced to the classical one. In contrast, in 1-D propagation, or a system like a layer of material in a cavity, we observe exponential localization. An intermediate case was showing a smooth transition between these two completely different behaviours.

We relate our findings to the standard theory of localization. Our discussion emphasises the new features which arise when the systems of interest are described by transfer matrices rather than Hamiltonians. This would be the natural approach in case of engineered quantum networks.

Acknowledgements

P.T. acknowledges the support by the Academy of Finland (projects 42588, 47140 and 48845 and Research Centre for Computational Science and Engineering, project 44892, Finnish Centre of Excellence Programme 2000-2005). The financial support by the Alexander von Humboldt foundation, MSM 210000018 (Czech Republic), GACR 202/01/0318 and EU IST-1999-13021 for I.J. is greatefully acknowledged.


text:
# Coriolis Force and Sagnac Effect

## Abstract

We consider the optical Sagnac effect, when the fictitious gravitational field simulates the reflections from the mirrors. It is shown that no contradiction exists between the conclusions of the laboratory and rotated observers. Because of acting of gravity-like Coriolis force the trajectories of co- and anti-rotating photons have different radii in the rotating reference frame, while in the case of the equal radius the effective gravitational potentials for the photons have to be different.

PACS numbers: 04.80.Cc, 04.20.Cv

Since its discovery at the beginning of the XX century the Sagnac effect  has play an important role in the understanding and development of fundamental physics (for a review see ). The Sagnac effect is the dependence of the interference pattern of the rotating interferometer on the direction and speed of rotation. This phenomenon is universal and manifested for any kind of waves, including the matter waves and has found a variety of applications for the practical purposes and in the fundamental physics . It has also a direct experimental verification for large distances in the experiments of clocks transported around Earth . For the laboratory observer the Sagnac effect seems to have simple explanation - because of the rotation the round-trip distance traveled by the waves co-rotating with the platform is grater than for the anti-rotating waves.

Some misunderstandings appear when considering the optical Sagnac effect from the viewpoint of an observer on the rotating disk (see for example ), where the both light beams travel the same distance. Since the round-trip time along a closed path on the disk is different for clockwise and anti-clockwise photons the non-inertial observer can conclude that the speed of the light measured locally on the disk depends on the direction and speed of rotation. On the other hand, in accordance with General Relativity consistently defined speed of light for any observer turns out to be exactly the same, just as in the case of an inertial reference frame .

In the context of the Sagnac effects the null result of the Michelson-Morley experiment is also not clear. Applying the same logic to Sun centered rotating frame in which Earth is fixed, one would expect different light speeds as seen from Earth. Another problem is explanation of the experiments where the effects of placing of transparent media on the paths of the beams were analyzed . No interference fringe change was observed between the light beams in air and in other media configurations. However, the change was observed when the media was stationary and the platform was rotated.

The Sagnac effect have been attempted by several authors to explain: using the General Relativity ; from the special relativistic Doppler effect at the mirrors ; by coupling the momentum of interfering particles to rotation ; introducing various types of relativistic transformations for rotating frames ; using the restricted formula of the 3-space element ; considering the non-invariance of the light speed for the rotating observers .

To make a round-trip on the disk the photon must be reflected by several mirrors and its path has the shape of a polygon. For example, in the classical Sagnac experiment  the light path was quadratic. The form of the trajectory is close to the circle when the number of sides of the polygon is large. In this paper we would like to introduce the fictitious gravitational potential and write the effective metric for the photons in the form

$$ds^2=c^2(12A/r)dt^2dr^2r^2d\varphi ^2dz^2.$$
(1)
The ”gravitational mass” $`A`$ is a constant simulating the reflection of photons. In order to have a closed trajectory of photons in (1) we do not change the spatial part of the metric.

The geodesic equation $`dk^\alpha /d\lambda +\mathrm{\Gamma }_{\mu \nu }^\alpha k^\mu k^\nu =0`$ for the photons moving on the circle ($`z=0,r=\rho =const`$) reduces to

$$\frac{dt}{d\lambda }=const,\frac{d\varphi }{d\lambda }=const,\frac{c^2A}{r^2}\left(\frac{dt}{d\lambda }\right)^2=r\left(\frac{d\varphi }{d\lambda }\right)^2,$$
(2)
where $`k^\nu =dx^\nu /d\lambda `$ is a wave 4-vector and $`\lambda `$ is some parameter varying along the ray.

This system and the Hamilton-Jacobi equation for the photons $`g_{\mu \nu }k^\mu k^\nu =0`$, which is equivalent to $`ds^2=0`$, gives the value of the radius of the trajectory $`\rho =3A`$, and the expression for the angular velocity $`d\varphi /dt=(c^2A/\rho ^3)^{1/2}`$ as well. Finally we find the expected value of the proper time $`\tau =\sqrt{g_{tt}}t=t/\sqrt{3}=2\pi \rho /c`$ which is needed by the photons to travel the circle of the radius $`\rho `$.

For the case of disks uniformly rotating in different directions with the angular velocities $`\pm \omega `$ the metric (1) transforms to

$$ds_\pm ^2=c^2(12A/rr^2\omega ^2/c^2)dt^2dr^2r^2d\varphi ^2\pm 2r^2\omega d\varphi dtdz^2.$$
(3)
Transformation of (1) to (3) is in the full accordance with numerous cyclotron experiments supporting the contention that (3) is the correct form of the metric of a rotating disk, but only for the inertial observer (despite of the fact that it corresponds to the Galilean composition law of velocities ).

The geodesic equations (2) now have the form

$$\frac{dt}{d\lambda }=const,\frac{d\varphi }{d\lambda }=const,\left(\frac{c^2A}{r^2}r^2\omega ^2\right)\left(\frac{dt}{d\lambda }\right)^2=r\left(\frac{d\varphi }{d\lambda }\right)^22r\omega \frac{dt}{d\lambda }\frac{d\varphi }{d\lambda }.$$
(4)
This system and the Hamilton-Jacobi equation, yields $`\rho =3A`$ and $`d\varphi /dt=(c^2A/\rho ^3)^{1/2}\pm \omega `$. In the first approximation in $`(\omega \rho /c)`$ the last relation gives the travel times $`t\tau 2\pi \rho ^2\omega /c^2`$ of co- and anti-rotated photons along the circle of the radius $`\rho `$. For the time lag we recover the well-known Sagnac formula

$$\delta \tau 4S\omega /c^2,$$
(5)
where $`S=\pi \rho ^2`$ is area of the projection of the trajectory spanned by the photons.

As in every stationary field, the clocks on the rotating body cannot be uniquely synchronized at all points and the metric (3) can not be used for a rotating observer. In the rotating frame two fictitious gravity-like forces appear, namely the centrifugal and Coriolis forces. This is an illustration of the equivalence principle, which asserts that gravity and the accelerated motion are locally indistinguishable. The centrifugal force does not affect photons. However, it causes the deformations of the mirrors and shifts the relaxation times for the electrons in the conductivity zone, which are forming the mirrors. Thus the centrifugal force changes the constant $`A`$. In the coordinates of the observer on the disk (1) can be written in the form

$$ds^2=c^2(12B/R)dT^2dR^2R^2d\mathrm{\Phi }^2dZ^2.$$
(6)
The ”mass” $`B`$ simulates the motion of the photons along the circle of the radius $`P`$, just as in the inertial case.

The motion of the photons under the influence of the Coriolis force  can be described by the following equation

$$d\stackrel{}{k}/dT=2[\stackrel{}{k}\stackrel{}{\mathrm{\Omega }}],$$
(7)
where $`\stackrel{}{\mathrm{\Omega }}`$ is some 3-vector which appears in a rotated reference system. The Coriolis force acting on the tangential photons co-rotating (or anti-rotating) with the disk is directed radially outward away from (or inward towards) the centre. It means that the photons are deflected out (or to) the centre as if they are gravitationally repelled (or attracted) by the disk. The result is that the light beam reflected from one mirror on the disk to the next one falls on the point that located before (or next) the one it was falling in the non-rotating case. Thus for a time interval measured by the clock at the mirror the photons cover less (or more) distances along the disk circumference, in the full agreement with the inertial observer.

For the circular motion the Coriolis force has only the radial component and (7) yields

$$dk^\mathrm{\Phi }/dT=\pm 2Rk^\mathrm{\Phi }\mathrm{\Omega }.$$
(8)
Then the geodesic equations of photons moving along the circular orbits ($`Z=0,R=P=const`$) have to be written in the form different from (2)

$$\frac{dT}{d\mathrm{\Lambda }}=const,\frac{d\mathrm{\Phi }}{d\mathrm{\Lambda }}=const,\frac{c^2B}{R^2}\left(\frac{dT}{d\mathrm{\Lambda }}\right)^2=R\left(\frac{d\mathrm{\Phi }}{d\mathrm{\Lambda }}\right)^2\pm 2R\mathrm{\Omega }\frac{d\mathrm{\Phi }}{d\mathrm{\Lambda }}\frac{dT}{d\mathrm{\Lambda }}.$$
(9)
This system and the Hamilton-Jacobi equation leads to

$$3B=P\pm 2\mathrm{\Omega }P^2/c,d\mathrm{\Phi }/dT=c/Pc/3B\pm \mathrm{\Omega }.$$
(10)
¿From the first equation we see that for the observer on the disk the radii of the trajectories for clockwise and anti-clockwise photons are different

$$P_\pm =\pm \frac{c}{4|\mathrm{\Omega }|}\left(\sqrt{1\pm \frac{24B|\mathrm{\Omega }|}{c}}1\right)3B\left(1\frac{3B\mathrm{\Omega }}{c}\right).$$
(11)
¿From (10) it follows that time intervals $`T_\pm =2\pi 3B/c2\pi (3B)^2\mathrm{\Omega }/c^2`$ needed by the beams to travel around the disk are different and give the time lag (5), where now $`S=\pi P^2`$. In order to have the same radius for the beams (this corresponds to $`P_+=P_{}`$ in (11)) in the metric (6) different ”gravitational masses” $`B_\pm `$ have to be introduced. From (11) one obtains $`B_+B_{}6B^2\mathrm{\Omega }/c`$ and the second equation in (10) gives (5) again.

In the framework of the tower experiment  the action of the Coriolis force on the light is equivalent to reception of co-rotated photons by a mirror from a higher gravitational potential (from the upstairs) and of anti-rotated photons from the lower gravitational potential (from the downstairs), since after the round reflection the cumulate angle factor disappears. For the mirror this effect looks like a redshift of the anti-rotated and a blueshift of the co-rotated photons as compared with the inertial ones.

The similar situation takes place for the case of the gravitational clock experiment  where two identical clocks were flown around Earth in different directions. The gravitational blue shift and the special relativistic slowing of the moving clocks are background affects and are the same for both of flying clocks. Only the Coriolis force distinguishes between the co- and anti-rotating senses of motion. For the co-rotated case the ’repulsive’ Coriolis force tend to offset the gravitational attraction effect and the clock is in a lower gravitational potential as seen by the Earth observer. For the anti-rotating case the gravitational and Coriolis forces have the same directions. It is well known that the clocks in a higher gravitational potential tick faster. As the result the anti-rotating clock will be ahead relative to the co-rotating one with the time shift (5). The detailed description of gravitational clock effect for the inertial observer is done in .

In the Michelson-Morley experiment the light makes a round trip along Earth surface, but along one half of the way it co-rotates and then anti-rotates with Earth. Therefore, the effects of ’repulsive’ and ’attractive’ Coriolis forces offset each other. One can also understand the effect of the transparent media on the light beams. The co-rotating medium just changes the ”gravitational mass” $`A`$ in (3) and does not affect (5), while the inertial medium changes the Coriolis force and the time lags as well.

Acknowledgements: Author would like to acknowledge the hospitality extended during his visits at the Theoretical Divisions of CERN and of Abdus Salam International Centre for Theoretical Physics where this work was done.


text:
# Bayesian Methods for Cosmological Parameter Estimation from Cosmic Microwave Background Measurements

## 1 Introduction

Recent determinations of the angular power spectrum of the anisotropy of the cosmic microwave background (CMB)  have created much excitement. These data may be used to estimate cosmological parameters. The maximum in the angular power spectrum around the multipole value $`l200`$ (’the first peak’) is consistent with inflation and vanishing mean spatial curvature . The lack of curvature adds compelling evidence to the case for the existence of some smoothly–distributed negative–pressure component which is possibly a cosmological constant (e.g., ). The constraints on the amplitude of the predicted second peak lead to lower limits on the baryon density  which are in conflict with the standard results from big bang nucleosynthesis (e.g., ); new data are eagerly awaited to provide a resolution. Fortunately, future satellite missions will provide angular power spectrum determinations with greatly improved precision .

The extraction of cosmological parameters from the CMB anisotropy data requires a long and computer–intensive analysis chain. Here we are interested in the last step in the chain which is a derivation of cosmological parameters from the inferred angular power spectrum, $`C_l`$.

This last step is computer–intensive as well and promises to be more so as the number of parameters in ones model increases. Current analysis exercises have included up to ten parameters; the scalar quadrupole and gravity wave perturbation normalizations ($`A_s`$ and $`A_t`$ ), the scalar and tensor power–law indices for primordial perturbations ($`n_s`$ and $`n_t`$), the reionization optical depth $`\tau `$, the spatial curvature $`\mathrm{\Omega }_k`$, the energy densities for baryonic matter ($`\omega _b`$), cold dark matter ($`\omega _{cdm}`$), neutrinos ($`\mathrm{\Omega }_\nu `$) and the vacuum ($`\mathrm{\Omega }_\mathrm{\Lambda }`$). Other logical cosmological parameters could include the number of neutrino families and their masses. Potentially infinite degrees of freedom lie in the primordial spectrum of perturbations.

The estimation of parameters can be thought of as an exercise in Bayesian inference. One starts with the likelihood function, namely the conditional probability distribution function (PDF) of the observation $`𝐳`$ *given* the unknown parameters, $`p(𝐳|𝜽)`$. For studies pertaining to the CMB, the likelihood used is an approximation to the probability of the power spectrum given the data, e.g. . If the likelihood function for a certain parameter vector $`𝜽_1`$ is higher than for a parameter vector $`𝜽_2`$, the observations are more likely to have occurred under $`𝜽_1`$ than under $`𝜽_2`$. Thus, the observations give higher plausibility to $`𝜽_1`$ and finding the parameter value that maximizes the likelihood seems to have some compelling logic. However, the only consistent way to quantify one’s uncertainty about an *unknown* parameter $`𝜽`$ is by specifying a probability distribution or equivalently a PDF for the parameter. After observing the data, this distribution should be updated to the PDF of the parameter *given* the data by incorporating the gained information. We are not interested in the PDF of the data given the  *unknown* $`𝜽`$ but in the PDF of $`𝜽`$ given the *known* data; we want the “inverse probability”, not $`p(𝐳|𝜽)`$ but $`p(𝜽|𝐳)`$. The only coherent approach to update a *prior* probability distribution with experimental information consists of calculating the *posterior* PDF via Bayes’ theorem:

$$p(𝜽|𝐳)=\frac{p(𝜽)p(𝐳|𝜽)}{m(𝐳)}$$
(1)
where $`m(𝐳)=p(𝐳|𝜽)f(𝜽)𝑑𝜽`$ is the marginal PDF of $`𝐳`$. The denominator can be regarded as a normalization constant since it is independent of $`𝜽`$. The posterior PDF thus combines prior with likelihood information. A Bayesian point estimate of the unknown parameter would be the posterior mean, median, or mode.

The usual difficulties of Bayesian inference apply here: namely, the challenge of high–dimensional integration. One needs to be able to perform high-dimensional integrations when calculating the normalization constant $`m(𝐳)`$, for instance, or when calculating the marginal PDF of a single component of the parameter vector $`𝜽`$ by integrating out all the other components, or when calculating posterior means. One of the main objectives of this paper is to demonstrate the potential of what statisticians refer to as simulation–based integration techniques for Bayesian posterior computation of CMB model parameters. These are called simulation–based, not because there is any simulation of the data, but because the posterior distribution is simulated. The particular simulation–based method we use here generates a set of parameter vectors whose distribution simulates the posterior distribution.

Analyses of cosmological parameter constraints from the most recent angular power spectrum determinations start with the calculation of the likelihood on an $`n_p`$–dimensional grid, where $`n_p`$ is the number of parameters . The next step is then the reduction of this large amount of data to a series of one– or two–dimensional probability distributions functions which can easily be plotted. This is achieved by marginalizing over the other variables; i.e., integrating over them. Note that some times (and in all applications for at least some of the variables) this marginalization is approximated by a maximization over the other variables.

The chief drawback to the grid–based approach is the exponential increase in computing times and storage requirements with increasing number of parameters. Even with $`n_p=10`$ (plus a handful of instrumental parameters) overly coarse grids and approximate treatment of some of the marginalizations are necessary to make the problem tractable. Including more parameters, for example ones that describe instrumental effects, foreground contaminants, or departures of the primordial power spectrum from a strict power–law, will make the grid–based approaches even more difficult if not impossible.

In this paper we implement the Bayesian approach to estimation of cosmological parameters using computer–intensive Markov chain Monte Carlo (MCMC) methods. Impediments to multidimensional integrations have been overcome by the progress made within the last decade in Bayesian computational technology via MCMC methods . Since its initial application in digital signal analysis  MCMC methods have revolutionized many areas of applied statistics and we expect there to be an impact on cosmological parameter estimation from CMB measurements as well. A distinct advantage of the MCMC approach is that computational time does not grow exponentially with parameter number, as it does for other methods . MCMC techniques have been applied in numerous areas, from science to economics . Applications of state–space modeling in finance, e.g. stochastic volatility models applied to time series of daily exchange rates or returns of stock exchange indices, easily have 1000–5000 parameters . Specially tailored MCMC algorithms can markedly improve the calculational speed . Hence, the MCMC approach to cosmological parameter estimation may provide the best strategy when testing complex models with numerous parameters.

We test our MCMC parameter estimation routine using simulated data. The likelihood is produced via a prototype fast calculator. The toy model depends on four parameters. We demonstrate the validity of the technique with this example. The goal of our research is to apply MCMC methods to real CMB data, and we are currently optimizing our routines for this multiparameter (~10) analysis. We are optimistic since experience has generally shown scaling of computing times with the number of parameters to be slower than exponential .

In section 2 we review methods that have been used to estimate cosmological parameters from CMB measurements, plus techniques for calculating likelihoods. In section 3 we describe the Bayesian approach to statistical inference and its implementation via Markov chain Monte Carlo methods. In section 4 we describe a method for applying MCMC methods to cosmological parameter estimation with CMB data. We test our method with a four parameter example, and our results are presented in Section 5. Section 6 presents our conclusions.

## 2 Current Statistical Methods with CMB Data

### 2.1 Parameter Estimation

Parameter estimation from CMB data is usually performed with a multi–step process, the last two steps of which are determining the power–spectrum from a map (or otherwise–pixelized data) , and then determining the cosmological parameters from the power spectrum. To date it has been impractical to go straight from the map to the parameters because of the computational expense of evaluating the probability of the map given the parameters thousands of times. Fortunately, the structure of the probability of the map given the power spectrum (the likelihood of the power spectrum) is much simpler than that of the map given the parameters. The simple structure allows for the mode to be found in a small number of iterations of a Levenberg–Marquardt–type search algorithm .

This “radical” compression of the information in the map to information in the power spectrum is possible because we assume that the signal in the maps is Gaussian and statistically isotropic. However, the uncertainties in the resulting power spectrum are not normally–distributed and the above search procedures do not allow one to completely characterize the distribution. Fortunately there are analytic approximations to the complete distribution whose parameters one can calculate with minor adaptations of the power–spectrum mode search algorithms .

The most recent attempts at determining cosmological parameters () are in fact attempts to determine the posterior probability distribution of the cosmological parameters under various prior assumptions. The first step is to evaluate the likelihood on a grid of cosmological parameters. To have constraints of reduced–dimensionality (one or two) suitable for plotting, one marginalizes over the other parameters. Sometimes marginalization is approximated by simply maximizing over the remaining parameters. This approximation is exact if the marginalized parameters are normally distributed.

Direct grid–based evaluations of the likelihood have computation time and storage requirements which rise exponentially with the number of parameters. Such approaches will be difficult to implement for models of greater complexity than have been studied so far. Further, the maximization approximation to marginalization can sometimes lead to spurious results. Numerical maximization techniques, such as the Levenberg–Marquardt method, are only guaranteed to find a local maximum. Once they have reached a local maximum they might get stuck in their search and not reach the global maximum. For this reason statisticians have applied *simulated annealing*, which is a technique for global optimization. Simulated annealing has been attempted for cosmological parameter estimation . Simulated annealing is related to MCMC as its core component is the Metropolis Hastings algorithm. In this method the parameter space is searched in a random way. A new parameter space point is reached with a probability that depends on the likelihood and an effective temperature term. In the limit where the temperature approaches zero the thermodynamics of this parameter space search finds the system approaching the maximum of the likelihood. Although these methods are applicable to high–dimensional problems, they can be very inefficient in certain situations. Furthermore, there is no guarantee that the global maximum will be reached in a finite time. The efficiency depends very much on specifying a good cooling schedule which involves the arbitrary and skillful choice of various cooling parameters.

### 2.2 Calculating the Likelihood

Implicit in the problem of parameter estimation, whether from a frequentist or Bayesian perspective, is the calculation of the likelihood of the observed power spectrum (determined from a map) given some cosmological model. Thus likelihood calculation requires us to be able to calculate the angular power spectrum for a given model (with its associated cosmological parameters). This calculational task is accomplished with computer codes such as CMBfast  or CAMB . These codes accept the cosmological parameters as input, and return the angular power spectrum of the CMB anisotropies, $`C_l`$. These software packages serve as the work–horses of current CMB parameter–determination efforts. For example, the likelihood has been calculated 30,311,820 times in order to cover a region in a ten–dimensional cosmological parameter space . In other studies the likelihood was evaluated as needed within the calculation .

## 3 Bayesian Posterior Computation via MCMC

Parameter estimation can be comprehensively described within the language of Bayesian inference. Application of Bayes’ theorem is well–suited to astrophysical observations . In Bayesian data analysis the model consists of a joint distribution over all unobserved (parameters) and observed (data) quantities. One conditions on the data to obtain the posterior distribution of the parameters. The starting point of the Bayesian approach to statistical inference is setting up a full probability model that consists of the *joint* probability distribution of all observables, denoted by $`𝐳=(z_1,\mathrm{},z_n)`$ and unobservable quantities, denoted by $`𝜽=(\theta _1,\mathrm{},\theta _d)`$. Using the notion of conditional probability, this joint PDF, $`p(𝐳,𝜽)`$can be decomposed into the product of the PDF of all unobservables, $`p(𝜽)`$, referred to as the *prior* PDF of $`𝜽`$, and the conditional PDF of the observables given the unobservables, $`p(𝐳|𝜽)`$, referred to as the sampling distribution or *likelihood*, i.e.

$$p(𝐳,𝜽)=p(𝜽)p(𝐳|𝜽).$$
(2)
The prior PDF contains all the information about the unobservables that is known from substantive knowledge and expert opinion *before* observing the data. All the information about the $`𝜽`$ that stems from the experiment is contained in the likelihood. In the light of the data, the Bayesian paradigm then updates the prior knowledge about $`𝜽`$, $`p(𝜽)`$, to the  *posterior* PDF of $`𝜽`$, $`p(𝜽|𝐳)`$. This is done via an application of Bayes’ theorem through conditioning on the observations

$$p(𝜽|𝐳)=\frac{p(𝜽,𝐳)}{m(𝐳)}p(𝜽)p(𝐳|𝜽)$$
(3)
where $`m(𝐳)=p(𝐳|𝜽)p(𝜽)𝑑𝜽`$ is the marginal PDF of $`𝐳`$ which can be regarded as a normalizing constant as it is independent of $`𝜽`$. The Bayesian approach is based on the likelihood function but also quantifies the uncertainty about the parameters through a joint prior distribution that summarizes the available information about the parameters before observing the data. In the light of the observations, the information about the unknown parameters is then updated via Bayes’ theorem to the posterior distribution which is proportional to the product of likelihood and prior density .

As already mentioned in the introduction, the main difficulty with the Bayesian approach to parameter estimation is high–dimensional integration. To calculate the normalizing constant of the joint posterior PDF, for instance, requires $`d`$–dimensional integration. Having obtained the joint posterior PDF of $`𝜽`$, the posterior PDF of a single parameter $`\theta _i`$ of interest can be obtained by integrating out all the other components, i.e.

$$p(\theta _i|𝐳)=\mathrm{}p(𝜽|𝐳)𝑑\theta _1\mathrm{}𝑑\theta _{i1}𝑑\theta _{i+1}\mathrm{}𝑑\theta _d.$$
(4)
Calculation of the posterior mean of $`\theta _i`$ necessitates a further integration, e.g. $`E[\theta _i|𝐳]=\theta _ip(\theta _i|𝐳)𝑑\theta _i`$. This procedure is referred to as marginalization.

Only in the simplest situations can these integrals be solved analytically. The main approximate techniques are normal and Laplace approximations based on asymptotics, quadrature approximations, Monte Carlo integrations, and stochastic simulation. MCMC methods belong to the last category. For an overview see . The deterministic techniques rely upon approximate normality and asymptotic results in the sense of the sample size growing to infinity. These techniques were mostly developed before the immense technological advances that enabled computer–intensive methods to be applied. The complexity of these techniques increases substantially with the dimension of the parameter space. In very broad terms, experience suggests that deterministic techniques provide good results for low–dimensional models. Similar comments are valid for non–iterative simulation techniques, since finding a suitable auxiliary distribution in rejections or importance sampling, for instance, becomes an extremely difficult task for high dimensions. As the dimension of the model increases, only iterative simulation–based integration techniques such as sampling–importance resampling (SIR) or MCMC provide adequate solutions. One major advantage of using a sampling–based approach to posterior computation is that once a sample from the posterior PDF, say $`(𝜽_1,\mathrm{},𝜽_N)`$ is available, we can employ this to estimate the posterior mean of each parameter by the sample average of the corresponding component, the marginal PDF’s using kernel density estimates, and correlation between parameters using the sample correlations.

To avoid a time–consuming direct sampling of the joint posterior we propose a MCMC method . Instead of generating a sequence of independent samples from the joint posterior, in MCMC a Markov chain is constructed whose equilibrium distribution is just the joint posterior. Thus, after running the Markov chain for a certain burn–in period, one obtains (correlated) samples from the limiting distribution, provided that the Markov chain has reached convergence.

One method for constructing a Markov chain is via the Metropolis–Hastings algorithm. The MH algorithm was developed by Metropolis et al.  and generalized by Hastings . It is a MCMC method which means that it generates a Markov chain whose equilibrium distribution is just the target density, here the joint posterior PDF $`p(𝜽|𝐳)`$. This means that after sampling from this Markov chain for a sufficiently long time to allow the chain to reach equilibrium, the samples can be regarded as samples from the joint posterior PDF.

The MH algorithm shares the concept of a generating PDF with the well–known simulation technique of *rejection sampling*, where a candidate is generated from an auxiliary PDF and then accepted or rejected with some probability. However, the candidate generating PDF, $`q(𝜽|𝜽_n)`$ can now depend on the current state $`𝜽_n`$ of the Markov chain. A new candidate $`𝜽^{}`$ is accepted with a certain  *acceptance probability* $`\alpha (𝜽^{}|𝜽_n)`$ also depending on the current state $`𝜽_n`$ given by:

$$\alpha (𝜽^{}|𝜽_n)=\mathrm{min}\{\frac{p(𝜽^{})p(𝐳|𝜽^{})q(𝜽_n|𝜽^{})}{p(𝜽_n)p(𝐳|𝜽_n)q(𝜽^{}|𝜽_n)},1\}$$
(5)
if $`(p(𝜽_n)p(𝐳|𝜽_n)q(𝜽^{}|𝜽_n))>0`$ and $`\alpha (𝜽^{}|𝜽_n)=1`$ otherwise.

The steps of the MH algorithm are therefore:

| Step 0: | Start with an arbitrary value $`𝜽_0`$ |
| --- | --- |
| Step $`n+1`$: | Generate $`𝜽^{}`$ from $`q(𝜽|𝜽_n)`$ and $`u`$ from $`U(0,1)`$ |
|  | If $`u\alpha (𝜽^{}|𝜽_n)`$ set $`𝜽_{n+1}=𝜽^{}`$ (acceptance) |
|  | If $`u>\alpha (𝜽^{}|𝜽_n)`$ set $`𝜽_{n+1}=𝜽_n`$ (rejection) |

The MH algorithm does not require the normalization constant of the target density. Note that the outcomes from the MH algorithm can be regarded as a sample from the invariant density only after a certain burn–in period. Although the theory guarantees convergence for a wide variety of proposal PDFs, it does not say anything about the speed of convergence, i.e. how long the ”burn–in” period should be. Convergence rates of MCMC algorithms are important topics of ongoing statistical research with little practical findings so far. There is no formula for determining the minimum length of an MCMC run beforehand, nor a method to confirm that a given chain has reached convergence. The only tests available are based on an empirical time series analysis of the sampled values and can only detect non–convergence. Thus, by performing a whole sequence of so–called convergence diagnostics with negative results, one only gains more confidence in one’s hope that the chain reached equilibrium but never a guarantee. For issues concerning convergence diagnostics, the reader is referred to Cowles and Carlin .

Various convergence diagnostics have been developed and are implemented in CODA . CODA is a menu–driven collection of functions for analyzing the output of the Markov chain. Besides trace plots and the usual tests for convergence, CODA calculates statistical summaries of the posterior distributions and kernel density estimates.

## 4 Applying MCMC Methods to CMB Anisotropy Data

Probability distribution functions for cosmological parameters given CMB anisotropy data can be computed in a Bayesian fashion with the MCMC serving as a means of conducting a proper marginalization over parameters. The implementation of the MCMC method is relatively straightforward. Instead of calculating the likelihood at uniform locations in the parameter space , one lets the MCMC do its intelligent walk through the space. Uniform *a priori* distributions for the parameters seem reasonable, so the MCMC could sample the parameter space defined by them. Since the likelihood function can not be written explicitly in terms of the cosmological parameters, but instead in terms of the CMB anisotropy power spectra, $`C_l`$, it is necessary to implement a Metropolis–Hastings MCMC routine.

The Markov chain can commence at a randomly selected position in parameter space $`(\theta _1^{(0)},\mathrm{},\theta _d^{(0)})`$. With the parameter set one can then utilize a program like CMBfast  or CAMB  to generate an angular power spectrum and a likelihood can then be calculated. New values for the parameters $`(\theta _1^{(1)},\mathrm{},\theta _d^{(1)})`$ are selected via sampling from the *a priori* distributions. However, these values are not necessarily *accepted* as new values. First their likelihoods are evaluated. Then the new values are *accepted* or rejected according to the following test; a random number, $`u`$, would be generated between 0 and 1. If $`u\mathrm{min}[1,(p(\theta _1^{(1)},\mathrm{},\theta _d^{(1)})p(𝐳|\theta _1^{(1)},\mathrm{},\theta _d^{(1)}))/(p(\theta _1^{(0)},\mathrm{},\theta _d^{(0)})p(𝐳|\theta _1^{(0)},\mathrm{},\theta _d^{(0)}))]`$ where $`p(𝐳|\theta _1^{(1)},\mathrm{},\theta _d^{(1)})`$ is the likelihood in terms of data $`𝐳`$ and cosmological parameters $`(\theta _1,\mathrm{},\theta _d)`$ then the new parameters is *accepted* into the chain, if not the next chain element has values equal to that of the previous state. A new set of parameters is then randomly sampled from the a priori distributions and the procedure continues.

The generated chain of parameter values forms the set from which the statistical properties would be derived. After running the Markov chain for a certain burn–in period (in order for the Markov chain to reach convergence) one obtains (correlated) samples from the limiting distribution. This process continues for a sufficiently long time (as determined by convergence diagnostics ).

After the *burn–in* the frequency of appearance of parameters represents the actual posterior density of the parameter. From the posterior density one can then create confidence intervals. Summary statistics are produced from the distribution, such as posterior mean and standard deviation. A cross–correlation matrix is also easily produced—which could prove to be a useful part of quantifying the multi–dimensional constraints.

The above method constitutes the simplest implementation of the Metropolis–Hastings method, that of an independence chain. We have just used an *acceptance probability* $`\alpha (𝜽^{}|𝜽)`$, defined by

$$\alpha (𝜽^{}|𝜽)=\mathrm{min}\{\frac{p(𝜽^{})p(𝐳|𝜽^{})q(𝜽|𝜽^{})}{p(𝜽)p(𝐳|𝜽)q(𝜽^{}|𝜽)},1\}$$
(6)
where the generating density $`q(𝜽^{}|𝜽)`$ is the uniform density over the parameter space and thus, in particular, independent of the current state. By using uniform priors, the posterior PDFs in the acceptance probability calculation reduce to the likelihoods. However, the efficiency of a Metropolis–Hastings algorithm depends crucially on the form of the generating density $`q(𝜽^{}|𝜽)`$. Just using a uniform distribution that does not even depend on the current state $`𝜽`$ is the simplest but probably least efficient way to accomplish the task. Even with a uniform distribution, the algorithm will be irreducible/aperiodic/reversible and thus the Markov chain will converge towards its stationary distribution.

A slightly better way might be to use a uniform distribution in a neighborhood of the current $`𝜽`$. Any prior information could be useful, such as correlations that one could use to specify a multivariate normal centered around the current $`𝜽`$ with a covariance matrix that takes said correlations into account. The optimization of the Metropolis–Hastings MCMC strategy will inevitably require experimentation with the generating density $`q(𝜽^{}|𝜽)`$. While this may require some detailed study, the benefit will be the ability to generate posterior distributions for a large number of cosmological parameters.

An approximate formula for the Fisher matrix of the cosmological parameters which has often been used for forecasting parameter errors for CMB experiments may provide us with a useful generating density. The Fisher matrix is the expectation value of the second–derivative of the log of the likelihood and is given by:

$$F_{pp^{}}\left(𝜽\right)=\underset{l}{}\frac{C_l}{\theta _p}\frac{C_l}{\theta _p^{}}\frac{1}{\sigma _l^2}$$
(7)
where the variance in each $`C_l`$ determination can be approximated by

$$\sigma _l^2=\frac{2}{(2l+1)f_{\mathrm{sky}}}\left(C_l+w^1B_l^2\right)^2$$
(8)
where $`f_{\mathrm{sky}}`$ is the fraction of sky observed, $`w^1=\sigma _{\mathrm{pix}}^2\mathrm{\Omega }_{\mathrm{pix}}`$, $`\sigma _{\mathrm{pix}}`$ is the standard error in each map pixel, $`\mathrm{\Omega }_{\mathrm{pix}}`$ is the size of the pixel and $`B_l`$ is the Legendre transform of the beam profile. The generating density would then be

$$q(𝜽^{}|𝜽)\mathrm{exp}\left(\underset{pp^{}}{}(\theta _p\theta _p^{})F_{pp^{}}\left(𝜽^{\mathbf{}}\right)(\theta _p^{}\theta _p^{}^{})/2\right)$$
(9)
Note that the Fisher matrix is always evaluated at the same $`𝜽^{}`$ to avoid having to calculate the derivatives in Eq. 7 at each iteration.

One may wish to know the posterior distribution of the parameters for a range of choices of priors. For example, in  the analysis was performed with 13 different priors each of which corresponded to a different choice about our knowledge of, e.g., the Hubble constant and the baryon density. For some parameter–determination methods adoption of a variety of priors is very costly since for each choice of prior the entire calculation must be redone. This is true, for example, of the Levenberg–Marquardt–type search algorithm used in . In contrast, as emphasized in , evaluation of the likelihood on a grid allows for rapid calculation of various posteriors with different prior assumptions.

Fortunately, the sampled values from a Markov chain constructed assuming one prior can be used for other priors as well using importance sampling . Suppose we have devised a MCMC sampler with stationary distribution $`p^{}(\theta |z)p_1(\theta )p(z|\theta )`$ with a prior $`p_1`$. Then we can estimate the expectation of an arbitrary function $`g(\theta )`$ of interest under the modified posterior $`p(\theta |z)p_2(\theta )p(z|\theta )`$, i.e. using a new prior $`p_2`$, by importance reweighting the output $`\theta _1,\mathrm{},\theta _N`$ from the chain with stationary distribution $`p^{}`$. Thus,

$$E_p[g(\theta )|z]\underset{i=1}{\overset{N}{}}\frac{w_ig(\theta _i)}{_{i=1}^Nw_i}$$

where the importance weight is $`w_i=p_2(\theta _i)/p_1(\theta _i)`$.

## 5 Test of Metropolis–Hastings Algorithm

We have successfully tested our MH method on simulated data. The ‘data’ we simulated is an angular power spectrum from $`l=100`$ to $`l=800`$ with normally distributed errors. The underlying model from which we made the realization is a cold dark matter model with baryonic matter density $`\omega _b=0.019`$, dark matter density $`\omega _d=0.154`$, $`A=100.0`$ (in some units) and $`n=1.0`$ . The errors were realized as uncorrelated Gaussian random variables with zero mean and variance, $`\sigma _l^2`$ given by Eq. 8 assuming a noise–free map covering 10% of the sky.

To evaluate the likelihood of the data at any given point in the four–dimensional parameter space we evaluated:

$$2\mathrm{ln}=\underset{l}{}\left(C_l(\omega _b,\omega _d,A,n)C_l^{\mathrm{data}}\right)^2/\sigma _l^2.$$
(10)
The most time–consuming step in the likelihood calculation is calculation of $`C_l(\omega _b,\omega _d,A,n)`$. We sped this up using a fast $`C_l`$ calculator described very briefly here and in more detail elsewhere . We use CMBfast to pre–compute $`\mathrm{\Delta }_l^2(k)`$ on an 8 by 8 grid of values of $`\omega _b`$ and $`\omega _d`$. The quantity $`\mathrm{\Delta }_l^2(k)`$ gives the contribution from each wave number to the angular power spectrum: $`C_l=(4\pi )^2𝑑kk^2\mathrm{\Delta }_l^2(k)P(k)`$ where $`P(k)`$ is the primordial spectrum of density perturbations which we parametrize as a power–law: $`P(k)=Ak^n`$. Thus, for a given $`\omega _b,\omega _d,A`$ and $`n`$ we perform a multi–linear interpolation over the $`\omega _b`$, $`\omega _d`$ grid of $`\mathrm{\Delta }_l^2(k)`$ values and then do the integral over $`k`$.

A simple Fortran code was used for the MH routine, and it can be obtained at http://physics.carleton.edu/Faculty/nelsonhome.html. This routine calls the fast likelihood calculator for the likelihood values. The  a priori distributions were all uniform, with ranges of $`0.015<\omega _b<0.025`$, $`0.1<\omega _d<0.2`$., $`0<n<2`$ and $`0<A<200`$. The candidate generating densities were uniform and centered at the current chain parameter values and extending with ranges of $`5\times 10^4`$ for $`\omega _b`$, $`5\times 10^3`$ for $`\omega _d`$, $`0.25`$ for $`A`$, and $`0.01`$ for $`n`$. The results presented here were from a run of 200000 iterations, with an acceptance ratio of 41%. This calculation took ~24 hours on a Sun Ultra 10 (440 MHz) workstation.

For our analysis we thinned the 200,000 cycle chain by accepting every 20th observation in order to avoid highly correlated values. Of the remaining 10,000 samples we used a burn–in of 1,000 which yields a final chain length of 9,000. Extensive convergence diagnostics were calculated for the four parameters using the CODA software . All chains passed the Heidelberger–Welch stationarity test. The Raferty–Lewis convergence diagnostics confirmed that the thinning and burn–in period were sufficient. Lags and autocorrelations within each chain were reasonably low for all parameters. These convergence diagnostics are summarized in , and references therein.

The trace plots and resulting kernel densities for the four parameters are shown in figure 1. Summary statistics including posterior mean, standard deviation, the time series standard error (the square root of the spectral density estimate divided by the sample size), and the $`2.5\%`$, $`50\%`$, and $`97.5\%`$ credible regions are listed in table 1. The cross–correlation matrix is presented in table 2.

## 6 Discussion

The computational demands of Bayesian inference with large numbers of parameters are best met with MCMC methods. These methods have demonstrated their importance in numerous applications. As cosmological models grow in complexity it will become necessary to use techniques such as those discussed here in order to handle marginalization of parameters.

The demonstration of the Metropolis–Hastings routine with the toy model was successful. We plan to apply this technique to real CMB anisotropy power spectrum measurements. With MCMC techniques one can easily extend this method to ~10 parameters. We are presently working on optimizing the speed of the likelihood calculator for the larger parameter number. Markov chain Monte Carlo methods provide a means of handling large parameter numbers, and maintain a rigorous approach to Bayesian inference. They may prove to be essential for certain CMB parameter–determination efforts in the near future.

## 7 Acknowledgements

This work was supported by the National Science Foundation, Royal Society of New Zealand Marsden Fund, the University of Auckland Research Committee, and Carleton College. LK thanks the Institut d’Astrophysique de Paris for their hospitality while some of this work was completed.

## References

Figure and Table Captions

Figure 1. Trace and kernel density plots of the marginal posterior distributions for the parameters $`A`$ (amp), $`\omega _b`$ (b), $`\omega _d`$ (d), and $`n`$ (tilt). The true parameters are $`A=100.0`$, $`\omega _b=0.019`$, $`\omega _d=0.154`$ and $`n=1.0`$.

Table 1. The posterior mean, standard deviation, time series standard error (SE), and the $`2.5\%`$, $`50\%`$ (Median) and $`97.5\%`$ credible regions of the parameters $`A,`$ $`\omega _b`$, $`\omega _d`$ and $`n`$ (cf. figure 1).

Table 2. The cross–correlation matrix of the parameters $`A,`$ $`\omega _b`$, $`\omega _d`$ and $`n`$ (cf. figure 1).

TABLE 1

| Parameter | Mean | SD | SE | $`2.5\%`$ | Median | $`97.5\%`$ |
| --- | --- | --- | --- | --- | --- | --- |
| $`A`$ | $`100.0`$ | $`0.729`$ | $`7.69x10^3`$ | $`99.1`$ | $`100.0`$ | $`102.0`$ |
| $`\omega _b`$ | $`0.0187`$ | $`5.28x10^4`$ | $`5.56x10^6`$ | $`0.0177`$ | $`0.0187`$ | $`0.0198`$ |
| $`\omega _d`$ | $`0.156`$ | $`2.27x10^3`$ | $`2.39x10^5`$ | $`0.151`$ | $`0.156`$ | $`0.160`$ |
| $`n`$ | $`0.971`$ | $`1.81x10^2`$ | $`1.90x10^4`$ | $`0.936`$ | $`0.971`$ | $`1.010`$ |

TABLE 2

| Parameter | $`A`$ | $`\omega _b`$ | $`\omega _d`$ | $`n`$ |
| --- | --- | --- | --- | --- |
| $`A`$ | 1.0 |  |  |  |
| $`\omega _b`$ | 0.220 | 1.0 |  |  |
| $`\omega _d`$ | 0.614 | 0.014 | 1.0 |  |
| $`n`$ | -0.189 | 0.420 | -0.440 | 1.0 |