text:
# Brane World Cosmology Without the 𝑍₂ Symmetry

## 1 Introduction

Recently there has been considerable interest in the novel suggestion that we live in a Universe that possesses more than four dimensions. The standard model fields are assumed to be confined to a hypersurface (or 3-brane) embedded in this higher dimensional space, in contrast the gravitational fields propagate through the whole of spacetime . In order for this to be a phenomenologically relevant model of our universe, standard 4D gravity must be recoverd on our brane. There are various ways to do this, the most obvious being to assume that the extra dimensions tranverse to our brane are compact. In this case gravity can be recovered on scales larger than the size of the extra dimensions . This is different from earlier proposals since the restrictions on the size of the extra dimensions from particle physics experiments no longer apply, as the standard model fields are confined to the brane. The extra dimensions only have to be smaller than the scale on which gravity experiments have probed, currently of order 1mm. Another way to recover 4D gravity at large distances is to embed a positive tension 3-brane into an $`AdS_5`$ bulk . In this scenario 4D gravity is obtained at scales larger than the $`AdS`$ radius. Randall and Sundrum showed that this could produce sensible gravity even if the extra dimension was not compact.

The cosmology of these extra dimension scenarios has been investigated and the Friedmann equation derived and shown to contain important deviations from the usual 4-dimensional case . Some inflationary models have been investigated , as have brane world phase transitions, topological defects and baryogenesis . Most brane world scenarios that involve one extra dimension assume a $`Z_2`$ symmetry about our brane, motivated by a model derived from M-Theory proposed by Horava and Witten . However, many recent papers examine models that are not directly derived from M-Theory: for example there has been much interest in the one infinite extra dimension proposal. The motivation for maintaining the $`Z_2`$ symmetry is seemingly no longer adequate and it is therefore interesting to analyse a brane world model without this symmetry and to assess its phenomenological implications. There have been some multi-brane scenarios suggested, which involve branes that although lying between two bulk space times with the same cosmological constant, do not possess a $`Z_2`$ symmetry of the metric itself . This approach generates an altered Friedmann equation as well as giving different bulk solutions. The cosmological solutions have not been analysed and it is therefore interesting to entertain the possibility that we live on such a brane and to determine whether it possesses a sensible cosmology. A different approach has been taken by  where they have looked at the effect of having different cosmological constants either side of the brane. Our approach is not equivalent to theirs, as will be discussed in more detail below.

In Section 2 of this paper we utilise the 5D Einstein’s equations to derive the Friedmann equation for such a non-$`Z_2`$ symmetric 3-brane and we then investigate the associated cosmological consequences of such a model. We show that the Friedmann equation acquires an extra term when there is no $`Z_2`$ symmetry, which can give rise to a period of expansion on the brane, and consider constraints on such a term from nucleosynthesis. These constraints apply to any scenario in which we reside on a 3-brane without the $`Z_2`$ symmetry, and demonstrate that scenarios such as  are physically unrealistic. In Section 3 we restrict ourselves to the infinite fifth dimensional case and solve Einstein’s equations in the bulk in order to generate the global solutions, examining the difference between them and the symmetric solutions and also demonstrating that they remain well-defined. Our conclusions are summarised in Section 4.

## 2 The Friedmann Equation with no $`Z_2`$ Symmetry

In this section we examine a cosmologically realistic positive tension 3-brane in 5 dimensions. We assume that the 5D cosmological constants either side of the brane are identical and then derive the general solution to Einstein’s equations without assuming that the metric across the brane is $`Z_2`$ symmetric. For the single brane case, this amounts to having non-$`Z_2`$ symmetric initial conditions with a symmetric bulk, however many of the results also apply to multi-brane scenarios such as .

Since we are interested in cosmological solutions, we take a metric of the form:

$$ds^2=n^2(t,y)dt^2+a^2(t,y)\gamma _{ij}dx^idx^j+b^2(t,y)dy^2,$$
(1)
where $`y`$ is the coordinate of the fifth dimension and we adopt a brane-based approach where the brane is the hypersurface defined by $`y=0`$. Here $`\gamma _{ij}`$ is a maximally symmetric 3-dimensional metric with $`k=1,0,1`$ parameterising the spatial curvature. We assume immediately that $`b^2(t,y)`$ is not a function of time and therefore $`y`$ can be scaled so that $`b(y)=1`$. The metric is found by solving the 5D Einstein’s equations,

$$G_{AB}=\kappa ^2T_{AB},$$
(2)
where we define $`\kappa ^2=1/\stackrel{~}{M}_5^3`$ where $`\stackrel{~}{M}_5`$ is the fundamental (reduced) 5D Planck Mass. The stress-energy-momentum tensor can be written as,

$$T^A{}_{B}{}^{}=T^A{}_{B}{}^{}|_{\mathrm{𝑏𝑟𝑎𝑛𝑒}}^{}+T^A{}_{B}{}^{}|_{\mathrm{𝑏𝑢𝑙𝑘}}^{}.$$
(3)
Again with cosmology in mind, we assume a homogeneous and isotropic geometry in the brane and this makes it possible to write the first term as:

$$T^A{}_{B}{}^{}|_{\mathrm{𝑏𝑟𝑎𝑛𝑒}}^{}=\delta (y)\mathrm{diag}(\rho _b,p_b,p_b,p_b,0).$$
(4)
The second term, which describes the same negative bulk cosmological constant $`\rho _B`$ either side of the brane, is of the form,

$$T^A{}_{B}{}^{}|_{\mathrm{𝑏𝑢𝑙𝑘}}^{}=\mathrm{diag}(\rho _B,\rho _B,\rho _B,\rho _B,\rho _B).$$
(5)
Note that here we have taken an alternative approach from  who assume different cosmological constants either side of the brane. Our lack of $`Z_2`$ symmetry enters when we consider the solution of the metric itself. As is shown in , for the setup defined above, any set of functions $`a(t,y)`$ and $`n(t,y)`$ satisfying,

$$\left(\frac{\dot{a}}{na}\right)^2=\frac{1}{6}\kappa ^2\rho _B+\left(\frac{a^{}}{a}\right)^2\frac{k}{a^2}+\frac{𝒞}{a^4},$$
(6)
together with $`G_{05}=0`$, will be solutions to all the Einstein’s equations locally in the bulk. The last term in (6) comes from the electric part of the 5D Weyl tensor (where $`𝒞`$ is some constant) and its existence and possible problematic effects will be discussed below. Now $`G_{05}=0`$ is satisfied if

$$n(t,y)=\frac{\dot{a}(t,y)}{\dot{a}(t,0)}.$$
(7)
Here, we have normalised $`n(t,y)`$ so that $`n(t,0)=1`$. In order to obtain the Friedmann Equation on the brane, we evaluate (6) at $`y=0`$. This is easily done except for the $`(a^{}/a)^2`$ term. To evaluate this we need the junction conditions,

$`{\displaystyle \frac{[a^{}]}{a_0}}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{3}}\rho _b,`$ (8)
$`{\displaystyle \frac{[n^{}]}{n_0}}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{3}}(3p_b+2\rho _b),`$ (9)
where $`[Q]=Q(0^+)Q(0^{})`$ and $`Q(0)=Q_0`$. Now instead of assuming the $`Z_2`$ symmetry $`yy`$ on the metric itself which would give $`a^{}(0^+)=a^{}(0^{})`$ and therefore $`[a^{}]=2a^{}(0^+)`$, we write,

$$a^{}(0^+)=a^{}(0^{})+d(t).$$
(10)
where $`d(t)`$ is some function of time only and has yet to be determined. $`d(t)`$ represents the asymmetry of the metric across the brane and the fact that it is not identically zero for all time is what is referred to as the brane being ‘non-$`Z_2`$ symmetric’ in the rest of this paper. Now in order to determine $`d(t)`$ we use (8) which gives,

$`a^{}(0^+)`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{6}}\rho _ba_0+{\displaystyle \frac{d(t)}{2}},`$ (11)
$`a^{}(0^{})`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{6}}\rho _ba_0+{\displaystyle \frac{d(t)}{2}}.`$ (12)
Assuming that $`a^2(0)=(a^2(0^+)+a^2(0^{}))/2`$ results in,

$$\frac{a_0^2}{a_0^2}=\frac{\kappa ^4}{36}\rho _b^2+\frac{d^2(t)}{4a_0^2}.$$
(13)
To find an expression for $`d(t)`$ we take the jump of the (5,5) component of Einstein’s equations as is done in ,

$$\frac{\overline{a}^{}}{a_0}p_b=\frac{1}{3}\rho _b\frac{\overline{n}^{}}{n_0},$$
(14)
where $`\overline{Q}=(Q(0^+)+Q(0^{}))/2`$. Replacing $`a`$ and $`n`$ using equations (10) and (7), shows that,

$$\dot{d(t)}=\frac{3p_b\dot{a_0}}{\rho _ba_0}d(t).$$
(15)
The energy conservation equation is derived directly from the junction conditions (8) and (9) as shown by ,

$$\dot{\rho _b}=3(\rho _b+p_b)\frac{\dot{a}_0}{a_0}.$$
(16)
By using this to solve the differential equation (15), it gives us the desired expression for $`d(t)`$,

$$d(t)=\frac{2F}{\rho _ba_0^3},$$
(17)
where $`F`$ is an integration constant which, when non-zero dictates to what extent the $`Z_2`$ symmetry is broken. Combining this with equation (6) and doing the usual replacements (first found by ) to obtain the standard $`H^2\rho `$ relation at late times: $`\rho _b=\rho +\rho _\lambda `$ and $`\kappa ^2\rho _B/6+\kappa ^4\rho _\lambda ^2/36=0`$, (where $`\rho _\lambda `$ is the brane tension and $`\rho `$ is the physical brane energy density) results in the Friedmann equation for a brane without the $`Z_2`$ symmetry,

$$\left(\frac{\dot{a}_0}{a_0}\right)^2=\frac{\kappa ^4\rho _\lambda }{18}\rho +\frac{\kappa ^4}{36}\rho ^2\frac{k}{a_0^2}+\frac{𝒞}{a_0^4}+\frac{F^2}{(\rho +\rho _\lambda )^2a_0^8}.$$
(18)
So the absence of the $`Z_2`$ symmetry gives rise to an extra term in the Friedmann equation. For a radiation dominated Universe where $`\rho =\gamma /a_0^4`$, the extra term behaves as $`F^2/\gamma ^2`$ as $`\rho \mathrm{}`$ and $`(F\rho /\gamma \rho _\lambda )^2`$ as $`\rho 0`$. Note that equation (18) is the same as that found by , but differs from the Friedmann equations found by . Therefore all the following analysis and results of this section also apply to the setup described in . In order to obtain standard cosmology at late times we need to make the identification,

$$\frac{\kappa ^4\rho _\lambda }{18}=\frac{8\pi G_4}{3}=\frac{1}{3\stackrel{~}{M}_4^2},$$
(19)
where we have used the reduced 4D Planck mass defined by $`\stackrel{~}{M}_4^2=M_4^2/8\pi `$ and will use the 5D reduced Planck mass defined by $`\stackrel{~}{M}_5^3=M_5^3/8\pi `$. This implies that the brane tension, $`\rho _\lambda =6\stackrel{~}{M}_5^6/\stackrel{~}{M}_4^2`$. Using this, the fact that $`\kappa ^2=1/\stackrel{~}{M}_5^3`$ and also defining the dimensionless constants $`f=F\stackrel{~}{M}_4^2/\gamma \stackrel{~}{M}_5^3`$, $`c=3\stackrel{~}{M}_4^2𝒞/\gamma `$ and assuming $`k=0`$, allows us to write (18) in terms of $`f`$, $`c`$, $`\stackrel{~}{M}_4`$, $`\stackrel{~}{M}_5`$ (and $`\rho _\lambda `$),

$$\left(\frac{\dot{a}_0}{a_0}\right)^2=\frac{1+c}{3\stackrel{~}{M}_4^2}\rho +\frac{1}{(6\stackrel{~}{M}_5^3)^2}\left[\rho ^2+\frac{f^2\rho _\lambda ^2\rho ^2}{\left(\rho +\rho _\lambda \right)^2}\right].$$
(20)
This shows that the expansion of the universe is initially dominated by the $`\rho ^2`$ term, while at late times the standard cosmology phase with the usual $`H^2\rho `$ behaviour is obtained. If $`f<\sqrt{8(1+c)}`$ the third, ‘$`f`$’-term of (20) is always less significant than the other terms, and the resulting cosmology is similar to a brane cosmology with a $`Z_2`$ symmetry.

If $`f>3+2c`$ there will be a period between the $`\rho ^2`$ and $`\rho `$ driven phases when the ‘$`f`$’-term is dominant. It is reasonable to suggest that $`f`$ is significantly greater than $`c`$, since as will be demonstrated later in this section, nucleosynthesis constrains $`c`$ to be $`0.2`$. Now in order to get an explicit solution for $`a_0(t)`$ so that we can ascertain the effect of the new $`f`$ term, we are forced to make several approximations. First we use,

$$\frac{\rho ^2\rho _\lambda ^2}{(\rho +\rho _\lambda )^2}=\{\begin{array}{cc}\rho ^2& \rho <\frac{1}{4}\rho _\lambda \\ \frac{1}{4}\rho _\lambda \rho & \frac{1}{4}\rho _\lambda <\rho <4\rho _\lambda \\ \rho _\lambda ^2& \rho >4\rho _\lambda \end{array},$$
(21)
to approximate the ‘$`f`$’-term. Now it is possible to solve for $`a_0(t)`$ in equation (20) which demonstrates that the time dependance of $`a_0(t)`$ in each of the three phases takes the form:

$$a_0(t)\{\begin{array}{cc}\gamma \left(\frac{4(1+c)}{3\stackrel{~}{M}_4^2}t^2+\frac{2\sqrt{1+f^2}}{3\stackrel{~}{M}_5^3}t\right)^{\frac{1}{4}}& \rho <\frac{1}{4}\rho _\lambda \\ \gamma \left(\frac{16(1+c)+f^2}{12\stackrel{~}{M}_4^2}t^2+\frac{2}{3\stackrel{~}{M}_5^3}t\right)^{\frac{1}{4}}& \frac{1}{4}\rho _\lambda <\rho <4\rho _\lambda \\ \gamma \left(\frac{(1+c)}{\rho _\lambda }\left[\mathrm{cosh}\frac{2\rho _\lambda f}{3\stackrel{~}{M}_5^3}t1\right]+\frac{1}{\rho _\lambda }\mathrm{sinh}\frac{2\rho _\lambda f}{3\stackrel{~}{M}_5^3}t\right)^{\frac{1}{4}}& \rho >4\rho _\lambda \end{array},$$
(22)
Evaluating the exact constants that replace the proportional signs in (22) is rather cumbersome, so in order to obtain an order of magnitude approximation of the effects of the ‘f’-term we instead ignore all subdominant terms in each phase of the universe, assume that $`f1`$ and then solve for $`a_0(t)`$ again using (20). The resulting evolution of the universe then divides into 5 phases as described below. In the following it is useful to define the time scale $`t_\mathrm{F}=\stackrel{~}{M}_4^2/(4\stackrel{~}{M}_5^3f)`$,

* PHASE 1 ($`0<t<t_\mathrm{F}`$): Initially (20) is dominated by the $`\rho ^2`$ term. This continues until $`t=t_\mathrm{F}`$ when $`\rho =f\rho _\lambda `$, after which the ‘$`f`$’-term becomes dominant. Until that happens,

$$a_0(t)=\left[\frac{2\gamma }{3\stackrel{~}{M}_5^3}\right]^{1/4}t^{1/4}$$
(23)
* PHASE 2a ($`t_\mathrm{F}<t<t_\mathrm{F}+t_\mathrm{I}`$), with $`t_\mathrm{I}=t_\mathrm{F}\mathrm{ln}(f/4)`$: For large $`\rho `$ the ‘$`f`$’ term can be approximated by a constant. During this period the universe expands exponentially, as in inflation,

$$a_0(t)=\left[\frac{\stackrel{~}{M}_4^2\gamma }{6\stackrel{~}{M}_5^6fe}\right]^{1/4}e^{Ht},H=\frac{\stackrel{~}{M}_5^3}{\stackrel{~}{M}_4^2}f.$$
(24)
* PHASE 2b $`(t_\mathrm{I}+t_\mathrm{F}<t<t_\mathrm{I}+7t_\mathrm{F}`$): For $`\rho \rho _\lambda `$ the ‘f’-term starts to decrease, and is approximately proportional to $`f^2\rho _\lambda \rho `$. During this phase,

$$a_0(t)=\left[\frac{\gamma f^2}{6\stackrel{~}{M}_4^2}\right]^{1/4}(tt_\mathrm{I}+t_\mathrm{F})^{1/2}.$$
(25)
* PHASE 2c $`(t_\mathrm{I}+7t_\mathrm{F}<t<t_\mathrm{B}=t_\mathrm{I}+(3+f^2/2)t_\mathrm{F}`$): For small $`\rho `$ the ‘f’-term is approximately $`f^2\rho ^2`$. It ceases to be the dominant term in (20) when $`\rho =2\rho _\lambda /f^2`$. Until then,

$$a_0(t)=\left[\frac{2\gamma f}{3\stackrel{~}{M}_5^3}\right]^{1/4}(tt_\mathrm{I}3t_\mathrm{F})^{1/4}.$$
(26)
* PHASE 3 $`(t>t_\mathrm{B})`$: Finally, at late times the $`\rho `$ term dominates (20) to give the standard cosmology,

$$a_0(t)=\left[\frac{4\gamma (1+c)}{3\stackrel{~}{M}_4^2}\right]^{1/4}(tt_\mathrm{B}+f^2t_\mathrm{F})^{1/2}.$$
(27)

The above solution for $`a_0(t)`$ presents a very different picture of the evolution of our universe: it has the unconventional early $`\rho ^2`$ behaviour as seen in most brane world models, but now this is broken up by a period of exponential expansion. During this extra phase the scale factor increases by a factor of $`f^{1/4}`$. Like inflation this could help to solve the flatness problem, however unlike inflation it is not followed by reheating and so cannot help with the horizon and monopole problems. Eventually, as expected, the standard cosmology is obtained. The above approximate solution (equations (23–27)), suggests that there will be no exponential expansion unless $`f4+2c`$. If $`f`$ is very large ($`f\stackrel{~}{M}_4^2/\stackrel{~}{M}_5^2`$) then $`f\rho _\lambda >\stackrel{~}{M}_5^4`$ and Phase 1, and some of the succeeding phases will be above the 5D Planck scale, and may not actually occur.

For $`\sqrt{8(1+c)}<f<3+2c`$ there will be a short period of ‘$`f`$’-term domination. This occurs after the time when the $`\rho `$-term of (20) starts to dominate the $`\rho ^2`$-term. At this stage the ‘$`f`$’-term is no longer approximately constant, and so there is no exponential expansion.

Nucleosynthesis provides constraints on both $`c`$ and $`f`$. If $`c`$ is nonzero then it gives rise to a term of the form $`𝒞/a_0^4`$ in the Friedman equation (18) which at the time of nucleosynthesis behaves like additional relativistic degrees of freedom. The energy density at this time is given as $`\rho (t_N)=g_{}\frac{\pi ^2}{30}T_N^4`$, where $`t_N`$ and $`T_N`$ are the time and temperature of Nucleosynthesis, and $`g_{}`$, the number of effective relativistic degrees of freedom is strongly constrained by the observed abundances of light elements which show that any deviation from $`g_{}`$ given by $`\mathrm{\Delta }g_{}`$ satisfies: $`\mathrm{\Delta }g_{}<2`$. This implies the following constraint for $`𝒞`$, as pointed out by ,

$$\frac{𝒞}{a_0^4(t_N)}\frac{\pi ^2}{30}\frac{\mathrm{\Delta }g_{}T_N^4}{3\stackrel{~}{M}_4^2}.$$
(28)
Using the fact that $`\rho (t_N)=\gamma /a_0^4(t_N)=g_{}\frac{\pi ^2}{30}T_N^4`$, leads to,

$$\frac{𝒞}{\gamma }\frac{\mathrm{\Delta }g_{}}{3\stackrel{~}{M}_4^2g_{}},$$
(29)
which implies that

$$c\frac{\mathrm{\Delta }g_{}}{g_{}}0.2,$$
(30)
where we have taken the standard values $`g_{}=10.75`$ and $`\mathrm{\Delta }g_{}<2`$. Hence it is reasonable to suggest that $`f`$ could play a cosmologically significant role but that $`c`$ cannot. Note that this constraint on $`c`$ is independant of the approximations used to derive equations (22) and (23–27) and comes purely from our knowledge of the restrictions on extra relativistic degrees of freedom at the time of nucleosynthesis.

We can also obtain rough restrictions on $`f`$ in terms of $`\stackrel{~}{M}_5`$ by demanding that standard cosmology is in place by the time of nucleosynthesis. This is equivalent to requiring that phase 2 is over well before nucleosynthesis begins, which implies that,

$$\frac{1+c}{3\stackrel{~}{M}_4^2}\rho (t_N)\frac{1}{(6\stackrel{~}{M}_5^3)^2}f^2\rho (t_N)^2.$$
(31)
At this time the universe is radiation dominated and $`\rho (t_N)`$ can be written in terms of the temperature at nucleosynthesis, $`T_N1`$MeV, and the number of effective relativistic degrees of freedom, $`g_{}`$,

$$\rho (t_N)=g_{}\frac{\pi ^2}{30}T_N^4.$$
(32)
Substituting this into (31), assuming $`c`$ to be negligible due to (30) and switching to the standard 4D and 5D Planck masses, gives the following relation between the 5D Planck mass $`M_5`$ and the $`Z_2`$ symmetry breaking parameter $`f`$,

$$M_5\left(\frac{g_{}\pi ^3}{45}T_N^4M_4^2f^2\right)^{1/6}30f^{1/3}\text{TeV}.$$
(33)
For the case where $`f=0`$, only phase 1 and 3 exist and demanding just that phase 3 is over well before $`t_N`$ gives, from (20),

$$\frac{1+c}{3\stackrel{~}{M}_4^2}\rho (t_N)\frac{1}{(6\stackrel{~}{M}_5^3)^2}\rho (t_N)^2.$$
(34)
which leads as above to:

$$M_530\text{TeV}.$$
(35)
This is a fairly weak bound on $`M_5`$. For the infinite extra dimension scenario, experiments testing gravity at small distances have already demonstrated using the corrections to Newton’s gravity law calculated by Randall and Sundrum  that $`M_5>10^5`$TeV. This experimental constraint however, is not applicable to compactified scenarios. Supposing that we do live in an infinite 5th dimension and that $`M_5`$ has a value that is just outside our experimental reach, we can then use (33) to constrain $`f`$; $`f10^{11}`$. In this case, the period of exponential expansion (Phase 2a) would only last for $`t_\mathrm{I}=M_4^2/(4M_5^3f)\mathrm{ln}f10^7\text{TeV}^1`$ and the scale factor would increase during this time by a factor of around $`500`$. Increasing $`M_5`$ relaxes the bound on $`f`$ and would appear to lead to more inflation. However, if $`M_5`$ is too large then some, or all, of the resulting expansion occurs at energies higher than the 5D Planck scale. Consequently the maximum inflation occurs for $`M_5=5\times 10^6`$ TeV and $`f=10^{17}`$, which leads to an expansion of only $`10^4`$. Unfortunately, this is not cosmologically relevant. These order of magnitude expansion rates due to a non-zero $`f`$ are dependant on the approximations made to generate equations (23-27), however an exact treatment would likely lead to the same conclusion: that the cosmological effect of $`f`$ is negligible.

We have therefore derived constraints on $`c`$ and $`f`$, showing both that $`c`$ must be negligible in size: $`c0.2`$ and that although $`f`$ can be large, its cosmological effect is insignificant in that it causes inflation that results in an increase in $`a_0(t)`$ only of order $`10^4`$.

Note that, from equations (11), (12), (17) and (31) the effect of the $`Z_2`$ breaking term decreases with increasing time such that the universe reverts to standard cosmology. This suggests that brane world scenarios where the physical Universe is on a brane without this symmetry are not viable after nucleosynthesis. The proposal made in , which described a setup that would solve the hierarchy problem despite our physical universe existing on a positive tension brane, is unfortunately found to be unrealistic due to this reason.

## 3 Global Solutions

After having examined the phenomenological effects resulting from the Friedmann equation of a non-$`Z_2`$ symmetric brane world, we will now solve the 5D Einstein equations in the bulk and hence derive the corresponding global solutions for $`a(t,y)`$ and $`n(t,y)`$ for the non-$`Z_2`$ symmetric brane in an infinite extra dimension. We derive this solution for the general case first and then for a specific cosmologically realistic brane. To do this we can adapt the previously known general global solution for a brane with tension $`\rho _b`$ and negative 5D bulk cosmological constant $`\rho _B`$ to the non-$`Z_2`$ symmetric case.

We know from the (0,0) component of Einstein’s equations that the new non-symmetric solution will have a form similar to the symmetric case as derived in ;

$$a^2(t,y)=a_0^2\left(A(t)\mathrm{cosh}\mu y+B(t)\mathrm{sinh}\mu |y|+C(t)\right),$$
(36)
where $`\mu =\sqrt{2\rho _B/(3\stackrel{~}{M}_5^3)}`$. The requirement $`a^2(t,0)=a_0^2`$ trivially implies that $`A(t)+C(t)=1`$ for all $`t`$. Now using equation (11) it can be seen that,

$`a^{}(t,0^\pm )`$ $`=`$ $`\left\{a_{sym}^{}(t,0^\pm )\right\}+{\displaystyle \frac{d(t)}{2}},`$ (37)
$`B(t)`$ $`=`$ $`\left\{{\displaystyle \frac{\rho _b}{\sqrt{6\stackrel{~}{M}_5^3\rho _B}}}\right\}\pm {\displaystyle \frac{d(t)}{a_0\mu }}.`$ (38)
Here the $`\pm `$ in the expression for $`B(t)`$ corresponds to the solution on either side of the brane and we use $`\{\mathrm{}\}`$ to denote the solution found in the $`Z_2`$ symmetric case. $`C(t)`$ is found from the differential equation for $`a^2(t,y)`$ which is derived from the Einstein equations, (see ),

$$C(t)=\frac{3\stackrel{~}{M}_5^3(\dot{a}_0^2+k)}{\rho _Ba_0^2}.$$
(39)
Rewriting this using the ‘new’ Friedmann equation (18) gives,

$$C(t)=\left\{\frac{1}{2}\left(1+\frac{\rho _b^2}{6\stackrel{~}{M}_5^3\rho _B}\right)+\frac{3\stackrel{~}{M}_5^3𝒞}{\rho _Ba_0^4}\right\}+\frac{3\stackrel{~}{M}_5^3F^2}{\rho _B\rho _b^2a_0^8},$$
(40)
and therefore $`A(t)`$ is trivially given by,

$$A(t)=\left\{\frac{1}{2}\left(1\frac{\rho _b^2}{6\stackrel{~}{M}_5^3\rho _B}\right)\frac{3\stackrel{~}{M}_5^3𝒞}{\rho _Ba_0^4}\right\}\frac{3\stackrel{~}{M}_5^3F^2}{\rho _B\rho _b^2a_0^8}.$$
(41)
Using (37) and (17) leads to,

$$B(t)=\left\{\frac{\rho _b}{\sqrt{6\stackrel{~}{M}_5^3\rho _B}}\right\}\pm \sqrt{\frac{6\stackrel{~}{M}_5^3}{\rho _B}}\frac{F}{\rho _ba_0^4}.$$
(42)
Again the $`\pm `$ signs in the expression for $`B(t)`$ give the two different solutions on either side of the brane. The solution for $`n(t,y)`$ in the non-symmetric case is found from the above solution for $`a(t,y)`$ by using equation (7) as before. It is easily seen that setting $`F`$ to zero in the above solutions recovers the $`Z_2`$ symmetric situation.

We are interested in these solutions for a cosmologically realistic brane, so we make the same substitutions as were made to generate (18) and also assume a radiation dominated Universe. Setting $`\rho =\gamma /a_0^4`$ where $`\gamma `$ is a constant, we obtain expressions for $`A(t)`$, $`B(t)`$ and $`C(t)`$ corresponding to a brane with a viable cosmology,

$`A(t)`$ $`=`$ $`1+\chi +{\displaystyle \frac{1}{2}}\chi ^2+c\chi +{\displaystyle \frac{f^2\chi ^2}{2\left(1+\chi \right)^2}},`$ (43)
$`B(t)`$ $`=`$ $`(1+\chi )\pm {\displaystyle \frac{f\chi }{(1+\chi )}},`$ (44)
$`C(t)`$ $`=`$ $`\chi {\displaystyle \frac{1}{2}}\chi ^2c\chi {\displaystyle \frac{f^2\chi ^2}{2\left(1+\chi \right)^2}}.`$ (45)
Where we have defined the dimensionless variable $`\chi =\stackrel{~}{M}_4^2\rho /6\stackrel{~}{M}_5^6`$. Equation (43) shows that the metric off the brane is greatly altered by a large value of $`f`$, in fact the ‘f’-terms dominate during the period when $`f\chi 1/f`$. During this time, $`a(t,y)`$ vanishes on the side of the brane corresponding to the negative sign in equation (38), at position $`y_0`$ which is given by:

$$y_0\frac{1}{\mu }\mathrm{ln}\left[\frac{(f+2)\chi +2}{(f2)\chi 2}\right],$$
(46)
and $`n(t,y)`$ will only become small on the opposite side at position $`y_0`$ given by:

$$y_0\frac{1}{\mu }\left|\frac{4(\chi ^21)}{f\chi (\chi 3)}\right|,$$
(47)
although the exact position of the zeros will depend upon the negleglected terms in (43). After this time, $`a(t,y)`$ no longer vanishes, however $`n(t,y)`$ vanishes roughly as for the symmetric case where $`y_0`$ is now given by:

$$y_0\frac{1}{\mu }\mathrm{ln}\left[1+\frac{2\chi +\sqrt{2\chi (2\chi +c)}}{\chi (3\chi +2c)}\right].$$
(48)
Calculation of the Ricci tensor and scalar show that they are both finite at the points where $`a(t,y)`$ and $`n(t,y)`$ vanish, which implies that these points just correspond to coordinate singularities. A similar result has been obtained in  for the symmetric case. Since there are no obvious problems with the global solutions, the only restrictions on $`f`$ and $`c`$ are those obtained in the previous section, which are applicable to all scenarios, not just the infinite extra dimensional one considered here.

We have seen how in a cosmological context, the various new features that are presented by a non-$`Z_2`$ symmetric brane world as opposed to an ordinary 4D Universe, are not necessarily dangerous and providing that the basic parameters $`M_5`$, $`c`$ and $`f`$ satisfy some fairly relaxed constraints, then the phenomenological predictions of the scenario proposed in this paper are still in agreement with current experimental observation. This has been shown to be true, however only for the cosmological case. While it is comforting to know that the Friedmann equation for a homogeneous, isotropic brane world gives sensible late time cosmology, the extra features such as the electric part of the 5D Weyl tensor and the extra degrees of freedom in the metric due to dropping the $`Z_2`$ symmetry condition, could cause major problems when pertubations to the background metric are considered. Although several papers have been produced addressing the issue (see ), the complete analysis remains to be done. Only a fully rigorous treatment would be able to show whether or not the Weyl tensor or the lack of a $`Z_2`$ symmetry would lead to pertubations growing in a disasterous manner and whether a possible choice of initial conditions could help prevent this.

## 4 Conclusions

We have derived the Friedmann equation and bulk solutions in a brane world scenario where there is no $`Z_2`$ symmetry of the metric across the brane. The bulk solutions were shown to be well behaved. Relaxing the $`Z_2`$ symmetry introduces an extra term that behaves as an effective cosmological constant at early times. Approximate solutions of the Friedmann equation on the brane were found with this extra term. If the asymmetry is sufficient this will introduce a period of exponential expansion during the early brany evolution.

Demanding that standard cosmology is in place by the time of nucleosynthesis constrains the $`Z_2`$ symmetry breaking parameter. This limits the amount of exponential expansion of the scale factor at early times to be of order $`10^4`$, which is not cosmologically significant. Thus this would only be a small mitigating factor in the flatness problem.

We also note that the effects of the $`Z_2`$ breaking term decrease with time. This is essential to ensure that the standard cosmology is recovered at late times. This suggests that the scenarios without this symmetry at late times, such as  are not viable.

## Acknowledgements

We would like to thank Pierre Binetruy, Tim Hollowood, Nathalie Deruelle and David Wands for discussions and PPARC for financial support.


format: text

text:
# 1 Introduction

## 1 Introduction

A spectral set $`\mathrm{\Omega }`$ in $`^n`$ is a closed set of finite Lebesgue measure such that $`L^2(\mathrm{\Omega })`$ has an orthogonal basis of exponentials $`\{e^{2\pi i\lambda ,x}:\lambda \mathrm{\Lambda }\}`$ restricted to $`\mathrm{\Omega }`$. Any discrete set $`\mathrm{\Lambda }`$ in $`^n`$ with this property is called a spectrum for $`\mathrm{\Omega }`$ and $`(\mathrm{\Omega },\mathrm{\Lambda })`$ is called a spectral pair.

In 1974 B. Fuglede  studied the problem of finding self-adjoint commuting extensions of the operators $`i\frac{}{x_1},\mathrm{},i\frac{}{x_n}`$ inside $`L^2(\mathrm{\Omega })`$, and related the existence of such an extension to $`\mathrm{\Omega }`$ being a spectral set. He formulated the following conjecture.

### Spectral Set Conjecture.

Let $`\mathrm{\Omega }`$ be a measurable set of $`^n`$ with finite Lebesgue measure. Then $`\mathrm{\Omega }`$ is a spectral set if and only if $`\mathrm{\Omega }`$ tiles $`^n`$ by translations.

Here $`\mathrm{\Omega }`$ tiles $`^n`$ with a set $`𝒯`$ of translations if $`\mathrm{\Omega }+𝒯=^n`$ and, for $`t,t^{}𝒯`$,

$$meas((\mathrm{\Omega }+t)(\mathrm{\Omega }+t^{}))=0\text{if}tt^{}.$$

Despite extensive study, this conjecture remains open in all dimensions, in either direction. For work on this problem see , , , , , , , , .

Lagarias and Wang  approached the spectral set conjecture in terms of tiling sets. They proved that a large class of tiling sets $`𝒯`$ had a universal spectrum $`\mathrm{\Lambda }`$ in the sense that every set $`\mathrm{\Omega }`$ that tiles $`^n`$ with the tiling set $`𝒯`$ was a spectral set with the spectrum $`\mathrm{\Lambda }`$. They considered periodic tiling sets of the form111General periodic tilings can always be reduced to this form by a linear transformation of $`^n`$. $`𝒯=N_1\times \mathrm{}\times N_n+𝒜`$ with $`𝒜^n`$, for some $`N_1,\mathrm{},N_n`$. Their results were obtained by reduction to questions about factorizations 222Factorizations of groups are defined in §2. of abelian groups. (A connection to such factorizations was originally noted in Fuglede .) In the one-dimensional case, they formulated the following conjecture.

### Universal Spectrum Conjecture.

Let $`𝒯=+\frac{1}{m}A`$ where $`A`$ reduced (mod $`m`$) admits some factorization $`AB=/m`$. Then $`𝒯`$ has a universal spectrum of the form $`m+\mathrm{\Gamma }`$, where $`\mathrm{\Gamma }`$.

In support of this conjecture, they proved \[13, Theorem 1.2\] that if the cyclic group $`_m`$ has a property called the strong Tijdeman property, then any tile set $`𝒯=+\frac{1}{m}𝒜`$ has a universal spectrum of the form $`\mathrm{\Lambda }=m+\mathrm{\Gamma }`$ for some $`\mathrm{\Gamma }`$. They showed that the strong Tijdeman property holds for many cyclic groups, and formulated the conjecture that all cyclic groups have the strong Tijdeman property. We do not define the strong Tijdeman property here (see ) but observe that a necessary requirement for its truth for $`G=/m`$ is the truth of the Tijdeman conjecture ($`modm`$) stated in §2. In addition to these results, Petersen and Wang \[15, Theorem 4.5\] also found other tiling sets $`𝒯`$ which have universal spectra.

Recently Coven and Meyerowitz  observed that a conjecture equivalent to Tijdeman’s Conjecture $`(modm)`$ had been made earlier, by Sands  in 1977. Sands proved this conjecture holds for all $`m`$ divisible by at most two distinct primes. Sands’ conjecture $`(mod)`$ was disproved by Szabó  in 1985, by a direct construction which applies to certain integers $`m`$ divisible by three or more distinct primes. The smallest counterexample covered by that construction is $`m=2^33^35^2=5400.`$ In §2 we present a counterexample to this conjecture for $`m=2^23^25^2=900,`$ found using similar ideas to the construction in Szabó . The counterexamples disprove the Tijdeman Conjecture ($`mod`$ $`m`$) for such $`m`$, and this shows that the method for proving the existence of universal spectra given in  does not work in general. In consequence, new methods are needed to resolve the universal spectrum conjecture.

In §3 we present a new sufficient condition for the existence of a universal spectrum for a periodic tiling set in $`^n`$ (Theorem 3.1). This criterion is easier to check computationally than a necessary and sufficient condition given in . It seems conceivable that the condition of Theorem 3.1 is actually necessary and sufficient; this remains an unresolved question. If so it might prove useful in resolving the spectral set conjecture for periodic tiling sets.

In §4 we apply the sufficient condition of §3 to show that the tiling sets $`𝒯_A=+\frac{1}{900}A`$ and $`𝒯_B=+\frac{1}{900}B`$ associated to the counterexample in §2 do have universal spectra. These sets give examples supporting the universal spectrum conjecture which are not covered by the methods of  and .

## 2 Counterexample to Tijdeman’s Conjecture

A factorization $`(A,B)`$ of the finite cyclic group $`G=/m`$, written

$$\overline{A}\overline{B}=G,$$
(2.1)
is one in which every element $`gG`$ has a unique representation

$$g=\overline{a}+\overline{b},\overline{a}\overline{A}\text{and}\overline{b}\overline{B}.$$

We write $`\overline{A}=\{\overline{a}_0,\mathrm{},\overline{a}_{n1}\}`$ and $`\overline{B}=\{\overline{b}_0,\overline{b}_1,\mathrm{},\overline{b}_{l1}\}`$ with $`ln=m`$. Let $`A=\{a_0,\mathrm{},a_{n1}\}`$ be a lifting of $`\overline{A}`$ to $``$, with $`\overline{a}_i=a_i`$ $`(modm)`$, and similarly let $`B`$ be a lifting of $`\overline{B}`$. Any factorization $`(A,B)`$ of $`/m`$ yields the direct sum decomposition of $``$, as

$$A𝒯_B=,$$
(2.2)
in which

$$𝒯_B=m+B$$
(2.3)
is a periodic set with period $`m`$. There has been extensive study of the structure of factorization of finite cyclic groups (more generally finite abelian groups) and of direct sum decomposition (1.2) of the integers $``$, the history of which can be found in Tijdeman , see also Coven and Meyerowitz . The original conjecture of Tijdeman \[18, p. 266\] is as follows.

### Tijdeman’s Conjecture.

If $`A𝒯=`$ with $`0A𝒯`$ and $`A`$ is a finite set with $`n`$ elements, with g.c.d. $`\{a:aA\}=1`$, then there exists some prime factor $`p`$ of $`n`$ such that all elements of $`𝒯`$ are divisible by $`p`$.

Hajós  and de Bruijn  showed for any direct sum decomposition $`A+𝒯=`$ where $`|A|`$ is finite, the infinite set $`𝒯`$ is periodic and necessarily has the form (2.3) and thus corresponds to some factorization (2.1) of a cyclic group $`/m`$ in which $`|A|`$ divides $`m`$. This permits Tijdeman’s conjecture to be reformulated as a series of conjectures for each finite cyclic group $`/m`$, as follows.

### Tijdeman’s Conjecture ($`modm`$).

If $`\overline{A}\overline{B}=/m`$ with $`\overline{0}\overline{A}\overline{B}`$ and if $`\overline{A}`$ lifts to set $`A`$ with $`0A`$ and g.c.d. $`\{a:aA\}=1`$, then any lifting $`B`$ with $`0B`$ has
g.c.d. $`\{b:bB\}1`$.

Tijdeman showed that this conjecture holds for any $`m`$ for which $`/m`$ is a “good” group in the sense of Hajós  and de Bruijn . The complete list of cyclic “good” groups are known to be those of order $`p^n`$ $`(n1)`$, $`pq`$, $`pqr`$, $`p^nq`$ $`(n>1)`$, $`p^2q^2`$, $`p^2qr`$ and $`pqrs`$, where $`p`$, $`q`$, $`r`$ and $`s`$ are distinct primes, c.f.  for references. As indicated in the introduction, this conjecture was actually made earlier by Sands, who proved it also holds in the cases $`m=p^nq^k`$ $`(n,k1)`$.

In 1985 Szabó  gave a construction which gave counterexamples to Sands’ conjecture for certain $`m`$ divisible by three or more distinct primes. Szabó actually constructs sets $`A`$ with $`0A`$ and $`\mathrm{gcd}(A)=1`$, which tile the integers, whose members are not uniformly distributed $`modk`$ for any $`k2`$. Coven and Meyerowitz\[3, Lemma 2.5\] observe that it follows that any tiling set $`C`$ for $`A`$ with $`0C`$ must have $`\mathrm{gcd}(C)=1,`$ hence cannot be contained in any subgroup of $`/m`$, where $`m`$ is the minimal period of $`C`$. The Szabó construction applies to $`m=m_1m_2\mathrm{}m_r`$ with $`r3`$ in which each $`m_i=u_iv_i`$ with the $`m_i`$ pairwise relatively prime, $`u_i`$ is the smallest prime dividing $`m_i`$ and each $`v_i4.`$ The smallest $`m`$ satisfying these conditions is $`m=2^33^35^2=5400.`$ Here we present a counterexample for $`m=2^23^25^2=900`$ which was found using similar ideas.

###### Theorem 2.1

Tijdeman’s Conjecture $`(mod900)`$ is false.

### Proof.

We take sets $`A`$ and $`B`$ with $`|A|=|B|=30`$. The set

$$A=\{0,36,72,108,144\}\{0,100,200\}\{0,225\}.$$
(2.4)
It has $`g.c.d.\{a:aA\}=1,`$ since g.c.d.$`\{2^23^2,2^25^2,3^25^2\}=1`$. We choose

$$B=\left\{\begin{array}{cccccccccc}\hfill 0& \hfill 30& \hfill 60& \hfill 126& \hfill 180& \hfill 210& \hfill 220& \hfill 240& \hfill 300& \hfill 306\\ \hfill 330& \hfill 360& \hfill 375& \hfill 390& \hfill 480& \hfill 486& \hfill 510& \hfill 520& \hfill 540& \hfill 570\\ \hfill 660& \hfill 666& \hfill 690& \hfill 750& \hfill 780& \hfill 820& \hfill 825& \hfill 840& \hfill 846& \hfill 870\end{array}\right\}.$$
(2.5)
We claim that $`\overline{A}\overline{B}=/900`$ is a direct sum. This can be verified by a calculation333More generally one can consider the construction in Szabó. (a short computer program). However

$$g.c.d.\{b:bB_0\}=1,$$

by considering the values 126, 220 and 375. Thus Tijdeman’s conjecture ($`mod900`$) is false.

### Remark.

Hajós  advanced a weaker conjecture concerning direct sum decompositions of a cyclic group $`G`$, which is that every factorization $`\overline{A}\overline{B}=G`$ is quasiperiodic. We say that a factorization is quasiperiodic if one of $`\overline{A}`$ or $`\overline{B}`$, say $`\overline{B}`$, can be partitioned into disjoint subsets $`\{\overline{B}_1,\mathrm{},\overline{B}_m\}`$ such that there is a subgroup $`H=\{h_1,\mathrm{},h_m\}`$ of $`G`$ with

$$\overline{A}+\overline{B}_i=\overline{A}+\overline{B}_1+h_i,1im.$$

The example $`(\overline{A},\overline{B})`$ above for $`G=/900`$ is quasiperiodic. The choices are $`H=\{0,300,600\}`$ and

$`\overline{B}_1`$ $`=`$ $`\{\begin{array}{cccccccccc}\hfill 0& \hfill 126& \hfill 180& \hfill 306& \hfill 360& \hfill 486& \hfill 540& \hfill 666& \hfill 820& \hfill 846\end{array}\}`$
$`\overline{B}_2`$ $`=`$ $`\{\begin{array}{cccccccccc}\hfill 30& \hfill 210& \hfill 220& \hfill 300& \hfill 390& \hfill 480& \hfill 570& \hfill 660& \hfill 750& \hfill 840\end{array}\}`$
$`\overline{B}_3`$ $`=`$ $`\{\begin{array}{cccccccccc}\hfill 60& \hfill 240& \hfill 330& \hfill 375& \hfill 510& \hfill 520& \hfill 690& \hfill 780& \hfill 825& \hfill 870\end{array}\}.`$
The quasiperiodicity conjecture remains open.

## 3 Criterion for Universal Spectrum

We formulate a sufficient condition for a universal spectrum for a periodic tiling set in $`^n`$, which is simpler to check than the necessary and sufficient condition given in \[13, Theorem 1.1\].

Given a finite set $`^n`$ let $`f_{}(\lambda )`$ denote the exponential polynomial

$$f_{}(\lambda ):=\underset{b}{}e^{2\pi i<\lambda ,b>}$$
(3.9)
and let $`Z(f_{})`$ denote its set of real zeros, i.e.

$$Z(f_{})=\{\lambda ^n:f_{}(\lambda )=0\}.$$
(3.10)
We recall the following criterion for a set $`\mathrm{\Lambda }`$ to be a spectrum, taken from .

###### Proposition 3.1

Let $`\mathrm{\Omega }=[0,\frac{1}{N_1}]\times \mathrm{}\times [0,\frac{1}{N_n}]+`$ where $`\frac{1}{N_1}\times \mathrm{}\times \frac{1}{N_n}`$ is a finite set. Suppose that $`\mathrm{\Gamma }^n`$ is a set of distinct residue classes $`(modN_1\times \mathrm{}\times N_n)`$, i.e. $`(\mathrm{\Gamma }\mathrm{\Gamma })(N_1\times \mathrm{}\times N_n)=\{\mathrm{𝟎}\}`$. Then $`\mathrm{\Lambda }=(N_1\times \mathrm{}\times N_n)+\mathrm{\Gamma }`$ is a spectrum for $`\mathrm{\Omega }`$ if and only if $`|\mathrm{\Gamma }|=||`$ and

$$\mathrm{\Gamma }\mathrm{\Gamma }Z(f_{})\{\mathrm{𝟎}\}.$$
(3.11)
### Proof.

This is \[13, Theorem 2.3\], after a linear rescaling of Euclidean space $`^n`$ of $`\mathrm{\Omega }`$ by a factor $`(\frac{1}{N_1},\mathrm{},\frac{1}{N_n})`$ and a corresponding dilation of Fourier space by a factor $`(N_1,\mathrm{},N_n)`$.

###### Theorem 3.1

$`(`$Universal Spectrum Criterion$`)`$. Let $`𝒯=^n+𝒜`$ where $`𝒜\frac{1}{N_1}\times \mathrm{}\times \frac{1}{N_n}`$ is a finite set, and suppose there exists some set $`\mathrm{\Omega }`$ such that $`\mathrm{\Omega }`$ tiles $`^n`$ by translations using the tiling set $`𝒯`$. Consider a set $`\mathrm{\Lambda }=(N_1\times \mathrm{}\times N_n)+\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }^n`$ such that the residue classes $`\mathrm{\Gamma }(modN_1\times \mathrm{}\times N_n)`$ are all distinct, i. e.

$$(\mathrm{\Gamma }\mathrm{\Gamma })(N_1\times \mathrm{}\times N_n)=\{0\}$$
(3.12)
Then $`\mathrm{\Lambda }`$ is a universal spectrum for $`𝒯`$ provided that $`|\mathrm{\Gamma }|=\frac{N_1N_2\mathrm{}N_n}{|𝒜|}`$, and

$$(\mathrm{\Gamma }\mathrm{\Gamma })Z(f_𝒜)=\mathrm{}.$$
(3.13)
### Proof.

By Theorem 3.1 of  a set $`\mathrm{\Omega }`$ tiles $`^n`$ by translations with a periodic tiling set $`𝒯`$ if and only if there exists some finite set $`\frac{1}{N_1}\times \mathrm{}\times \frac{1}{N_n}`$ giving a factorization

$$𝒜=_{N_1}\times \mathrm{}\times _{N_n}=(\frac{1}{N_1}/)\times \mathrm{}\times (\frac{1}{N_n}/),$$

in which case $`\stackrel{~}{\mathrm{\Omega }}=[0,1]^n+`$ also has $`𝒯`$ as a tiling set.

By Theorem 1.1 of  it suffices to verify that $`\mathrm{\Lambda }`$ is a spectrum for each set

$$\mathrm{\Omega }_{}=[0,\frac{1}{N_1}]\times \mathrm{}\times [0,\frac{1}{N_n}]+$$

where

$$𝒜=(\frac{1}{N_1}/)\times \mathrm{}\times (\frac{1}{N_n}/)$$

is a direct sum decomposition. This property shows that

$$\mathrm{\Omega }_{}+𝒜=[0,\frac{1}{N_1}]\times \mathrm{}\times [0,\frac{1}{N_n}]++𝒜$$

is a fundamental domain for the $`n`$torus $`^n/^n.`$ It follows that the Fourier transform

$$\widehat{\chi }_{\mathrm{\Omega }+𝒜}(\lambda )=_{\mathrm{\Omega }+𝒜}e^{2\pi i<\lambda ,x>}𝑑x$$

for $`\lambda =k^n`$ satisfies

$$\widehat{\chi }_{\mathrm{\Omega }+𝒜}(k)=\{\begin{array}{ccc}1\hfill & \text{for}\hfill & k=\mathrm{𝟎},\hfill \\ 0\hfill & \text{for}\hfill & k^n\{0\}.\hfill \end{array}$$
(3.14)
Now

$`\widehat{\chi }_\mathrm{\Omega }(\lambda )`$ $`=`$ $`{\displaystyle _\mathrm{\Omega }}e^{2\pi i<\lambda ,x>}𝑑x={\displaystyle \underset{b}{}}e^{2\pi i<\lambda ,b>}{\displaystyle \underset{j=1}{\overset{n}{}}}{\displaystyle _0^{\frac{1}{N_j}}}e^{2\pi i\lambda _jx_j}𝑑x_j`$
$`=`$ $`f_{}(\lambda ){\displaystyle \underset{j=1}{\overset{n}{}}}{\displaystyle \frac{e^{\frac{2\pi i\lambda _j}{N_j}}1}{2\pi i\lambda _j}},`$
which gives

$$\widehat{\chi }_{\mathrm{\Omega }+𝒜}(\lambda )=f_𝒜(\lambda )f_{}(\lambda )\underset{j=1}{\overset{n}{}}\left(\frac{\mathrm{sin}\frac{\pi \lambda _j}{N_j}}{\pi \lambda _j}\right)e^{\pi i\frac{\lambda _j}{N_j}}.$$
(3.15)
Comparing (3.14) and (3.15) yields

$$f_𝒜(k)f_{}(k)=0\text{if}k^n(N_1\times \mathrm{}\times N_n).$$
(3.16)
Thus we obtain

$$f_{}(k)=0\text{if}kZ(f_𝒜)(^n(N_1\times \mathrm{}\times N_n)).$$
(3.17)
By hypothesis $`\mathrm{\Lambda }=(N_1\times \mathrm{}\times N_n)+\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }^n`$, $`|\mathrm{\Gamma }|=\frac{N_1\mathrm{}N_n}{|𝒜|}`$,

$$(\mathrm{\Gamma }\mathrm{\Gamma })(N_1\times \mathrm{}\times N_n)=\{\mathrm{𝟎}\},$$
(3.18)
and $`(\mathrm{\Gamma }\mathrm{\Gamma })Z(f_𝒜)=\mathrm{}`$. Thus $`|\mathrm{\Gamma }|=||`$. We claim that

$$\mathrm{\Lambda }\mathrm{\Lambda }Z(f_{})\{\mathrm{𝟎}\}.$$
(3.19)
This claim holds since $`\mathrm{\Lambda }\mathrm{\Lambda }^n`$, then noting that $`Z(f_{})`$ contains all points of $`^n(N_1\times \mathrm{}\times N_n)`$ not in $`Z(f_𝒜)`$ by (3.17), while (3.18) takes care of points in $`N_1\times \mathrm{}\times N_n`$. Now Proposition 3.1 shows that $`\mathrm{\Omega }`$ has $`\mathrm{\Lambda }`$ as a spectrum, and the theorem follows.

### Remarks.

(1). The main hypothesis in Theorem 3.1 is (3.13), which requires determining the finite set $`\mathrm{\Gamma }\mathrm{\Gamma }^n`$ and evaluating $`f_𝒜`$ at these points. This condition is computationally simpler to check than the criterion of \[13, Theorem 1.1\], which requires determining all the complementing sets $``$ to $`𝒜`$.

(2). It seems conceivable that the sufficient conditon of Theorem 3.1 might also be a necessary condition for a universal spectrum of the given form $`\mathrm{\Lambda }=m+\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }`$. To show this one would have to show that for each integer vector $`𝐤Z(f_𝒜)(^nN_1\times \mathrm{}\times N_n)`$ there exists some set $``$ with $`𝒜=(\frac{1}{N_1}/)\times \mathrm{}\times (\frac{1}{N_n}/)`$ such that $`𝐤Z(f_{})`$.

## 4 Universal Spectra

We apply the criterion of Theorem 3.1 to show that the the tiling sets $`+\frac{1}{900}A`$ and $`+\frac{1}{900}B`$ associated to this counterexample in §2 both have universal spectra.

###### Theorem 4.1

Let $`𝒯_A=+\frac{1}{900}A`$ with $`A=\{0,36,72,108,144\}\{0,100,200\}\{0,225\}`$. Then $`𝒯_A`$ has the universal spectrum $`\mathrm{\Lambda }_A=900+A`$.

### Proof.

We apply Theorem 3.1 with $`n=1`$ and $`N_1=900`$. Then $`𝒜=\frac{1}{900}A`$ and $`=\frac{1}{900}B`$ with $`B`$ given in (2.5) gives a factorization $`𝒜=\frac{1}{900}/`$. A calculation gives

$$f_𝒜(\lambda )=\left(\frac{1e^{\frac{2\pi i\lambda }{5}}}{1e^{\frac{2\pi i\lambda }{25}}}\right)\left(\frac{1e^{\frac{2\pi i\lambda }{3}}}{1e^{\frac{2\pi i\lambda }{9}}}\right)\left(\frac{1e^{\frac{2\pi i\lambda }{2}}}{1e^{\frac{2\pi i\lambda }{4}}}\right).$$
(4.20)
It follows that the set $`Z(f_𝒜)`$ consists of all integers $`k`$ such that one or more of the following three conditions hold:

* 5 divides $`k`$ and 25 doesn’t divide $`k`$.
* 3 divides $`k`$ and 9 doesn’t divide $`k`$.
* 2 divides $`k`$ and 4 doesn’t divide $`k`$.

The set $`\mathrm{\Lambda }_A=900+A`$ has $`\mathrm{\Lambda }`$ and $`(AA)900=\{0\}`$, since all $`aA`$ have $`0a900`$ and are distinct. Also $`|A|=130=\frac{900}{|𝒜|}`$. To apply Theorem 3.1 it remains to verify

$$(AA)Z(f_𝒜)=\mathrm{}.$$
(4.21)
While $`a_iA`$ for $`i=1,2`$ as

$$a_i=36k_i+100l_i+225m_i$$

with $`0k_i4`$, $`0l_i2`$ and $`0m_i1`$. Then

$$a_1a_2=36(k_1k_2)+100(l_1l_2)+225(m_1m_2)$$

with $`|k_1k_2|<5`$, $`|l_1l_2|<3`$ and $`|m_1m_2|<2`$. Thus if 5 divides $`a_1a_2`$, then 5 divides $`k_1k_2`$ so $`k_1=k_2`$, and we conclude 25 divides $`a_1a_2`$. By similar arguments if 3 divides $`a_1a_2`$ then 9 divides $`a_1a_2`$, while if 2 divides $`a_1a_2`$ then 4 divides $`a_1a_2`$. Thus none of (i), (ii), (iii) hold, and (4.21) follows.

### Remark.

The proof of Theorem 4.1 easily generalizes to the sets $`A`$ as appearing in the general construction of Szabó : All such tiling sets $`+\frac{1}{m}A`$ have a universal spectrum.

###### Theorem 4.2

Let $`𝒯_B=+\frac{1}{900}B`$ with $`B`$ given by (2.5). Then $`𝒯_B`$ has the universal spectrum $`\mathrm{\Lambda }_B=900+B`$.

### Proof.

We do not have a conceptual proof of this result; however the conditions of Theorem 3.1 can be verified by a direct calculation (on the computer.) A key fact is that $`Z(f_{})900`$ and that the complement of the set $`(f_{})`$ in $``$ is exactly444To verify this, by (3.18) it suffices to check that $`f_{}(k)0`$ whenever $`kZ(f_𝒜)`$. $`Z(f_𝒜)900`$. Because this fact holds, one can prove $`\mathrm{\Lambda }_B`$ is a universal spectrum for $`𝒯_B`$ by checking that

$$BBZ(f_𝒜)\{0\}.$$

Since $`Z(f_𝒜)`$ is given by conditions (i)–(iii) in the proof of Theorem 4.1, it suffices to verify that every nonzero element of $`BB`$ satisfies one of (i)–(iii). This can be done by hand.

### Remarks.

In Theorem 4.1 and Theorem 4.2 the universal spectrum $`\mathrm{\Lambda }`$ exhibited is a scaled version of the tiling set $`𝒯`$. This fact is special to these examples, and cannot hold in general. Any tiling set $`𝒯=+\frac{1}{m}𝒞`$ with the property that $`\mathrm{\Lambda }=m+𝒞`$ is a universal spectrum must have $`m=|𝒞|^2`$, by Proposition 3.1.

(2). The proof of Theorem 4.1 easily generalizes to apply to the sets $`A`$ appearing in the general construction of Szabó : All such tiling sets $`+\frac{1}{m}A`$ have a universal spectrum.

### Acknowledgment.

The first author is indebted to J. A. Reeds for helpful computations, and especially to M. Szegedy for computations and references.

AT&T Labs–Research, Florham Park, NJ 07932-0971, USA
E-mail address: jcl@research.att.com

Dept. of Mathematics, Univ. of Bahrain, P.O. Box 30238 , Isa Town, BAHRAIN
E-mail address: sszabo7@hotmail.com


text:
# Noisy Monte Carlo revisited

## I Introduction

When sampling by Monte Carlo the partition function $`Z=dUe^{H(\{U\})}`$, the most common algorithm is that of Metropolis : at each step, starting from the current configuration $`\{U\}`$, a candidate configuration $`\{U^{}\}`$ is proposed, and it is accepted with probability

$$P_{acc}=\mathrm{min}(1,\mathrm{e}^{(\mathrm{H}(\{\mathrm{U}^{}\})\mathrm{H}(\{\mathrm{U}\}))}).$$
(1)
This acceptance test is realized by comparing the right-hand side of (1) to a random number uniformly distributed in $`[0,1]`$. This seems like a waste of information: why compute $`H(\{U^{}\})`$ exactly, then compare it with a random number? It should be sufficient to estimate it. Indeed, this logical proposition has been studied several times . Two difficulties have been identified, both caused by the non-linear relationship between the energy $`H`$ and the probability $`e^H`$: $`(i)`$ what is needed is an unbiased estimate of $`e^H`$, which must be obtained from unbiased estimate(s) of $`H`$; $`(ii)`$ to be interpreted as a probability, the noisy estimator of $`P_{acc}`$ must be bounded, and in particular stay positive. Difficulty $`(i)`$ was overcome in , which however showed that violations of $`(ii)`$ caused intolerable systematic errors unless the amount of noise in the estimate of $`H`$ was minuscule. Difficulty $`(ii)`$ was overcome in , which showed that exact results could be obtained even in the presence of a large amount of noise in the estimate of $`H`$. Ref., however, introduces an infinite number of auxiliary variables, so that it may take an infinite amount of work to bring these auxiliary variables to equilibrium. And tests of the method are performed on a toy model with 5 degrees of freedom only, whose relevance may be questioned. Here, we simplify the method of , by introducing only 1 auxiliary variable per term in $`H`$. Moreover, we separate $`H`$ into a leading part to be calculated exactly, and a sum of small correction terms, which we treat stochastically. This separation is essential: because stochastic estimates are used for correction terms only, large amounts of noise can be tolerated. As a consequence, our algorithm is a very efficient approach to the simulation of complicated Hamiltonians.

Consider a generic Hamiltonian of the type

$$H=\underset{k=0}{\overset{m}{}}c_kW_k$$
(2)
where as $`k`$ increases, $`|c_k|`$ decreases and the successive terms $`W_k`$ typically become less and less local. For instance, in a spin model $`\{\stackrel{}{\sigma _i}\}`$, $`W_0`$ would be the nearest-neighbour interaction $`_{<ij>}\stackrel{}{\sigma _i}.\stackrel{}{\sigma _j}`$, $`W_1`$ would represent next-nearest-neighbour interactions, etc… Here, we will illustrate our method for lattice gauge theory. In that context, $`W_k`$ are the traces of Wilson loops of increasing size: $`W_0=_{x,\mu ,\nu }Tr_4U`$ around elementary plaquettes, $`W_{1,2,3}`$ correspond to different geometries of 6-link loops, etc… It is often the case that one would like to study a Hamiltonian of type (2) resulting from an expansion, be it perturbative , non-perturbative , or based on the fixed point of a renormalization group transformation . In all these situations, the expansion is truncated to a maximal order $`m`$ dictated by technical reasons. As $`k`$ increases in (2), the number of geometrically equivalent terms grouped into $`W_k`$ increases exponentially: in a spin model on a hypercubic lattice in $`d`$ dimensions, each spin has $`2d`$ nearest-neighbours (these interactions are grouped into $`W_0`$), $`d(d1)/2`$ next-nearest neighbours (grouped into $`W_1`$), $`d(d1)(d2)/6`$ 3rd-neighbours, etc… This combinatoric explosion normally makes the simulation of extended Hamiltonians prohibitively expensive. This is the reason for a truncation to very low $`m`$, often taken to be 1 or 2. However, in most cases the couplings $`c_k`$ in (2) decrease exponentially with $`k`$, so that the overall Hamiltonian is dominated by $`W_0`$, with small corrections. In lattice field theory, this is actually required if the Hamiltonian is to make sense and tend to a local operator as the continuum limit is approached. By making use of stochastic methods to estimate the correction terms $`W_k,k1`$, we aim at postponing the combinatoric explosion of the simulation costs incurred when including higher terms $`W_k`$. This opens the possibility of studying numerically much more complicated Hamiltonians including higher-order correction terms. In lattice field theory, these correction terms are crucial to suppress discretization errors, and form the building blocks of so-called “improvement” strategies. Also the inclusion of higher-order terms can be very useful in the approaches to the fermion determinant simulations involving the loop expansion .

We present our method in Section II, and illustrate it in Section III with simulations of a 5-loop perturbatively improved action for $`SU(2)`$ lattice theory. We conclude with prospective applications of our method, in particular for dynamical fermion QCD simulations.

## II Noisy Monte Carlo: the method

Given the Hamiltonian (2) let us suppose that the terms $`c_kW_k`$ are nonpositive starting from $`k=1`$:

$$k1:c_kW_k(U)0U$$
(3)
This can be easily arranged by adding to each term of the Hamiltonian a nonessential constant. Here $`U`$ are the fields of the model under consideration. The key idea of the method is to estimate the contribution of the terms $`W_k(U),k1`$ stochastically by introducing auxiliary fields. This will lead to a significant reduction of computational effort if the coefficients $`c_k,k1`$ are small enough. In all cases, the algorithm remains exact.

We introduce auxiliary fields $`\sigma _k,k1`$ (associated with the terms $`W_k`$), which can take two values: 0 and 1. Using the identity

$$a+b=\underset{\sigma =0,1}{}[a\delta _{\sigma ,0}+b\delta _{\sigma ,1}]$$
(4)
we represent the probability $`e^H`$ in the form

$$e^H=P_0[U]P_1[U,\sigma ]$$
(5)
where

$$P_0[U]=e^{c_0W_0(U)};P_1[U,\sigma ]=\underset{k=1}{\overset{m}{}}\underset{\sigma _k=0,1}{}[\delta _{\sigma _k,0}+\delta _{\sigma _k,1}(e^{c_kW_k(U)}1)]$$
(6)
The r.h.s. of (5) can be interpreted as the joint probability distribution for the original fields of the model and the new $`\sigma `$ fields. Due to the inequalities (3) this distribution is well defined: $`P_1[U,\sigma ]0\{U,\sigma \}`$, and the probabilities for $`\sigma _k`$ to take value 0 or 1 when the $`U`$ fields are fixed lie in the interval :

$$p_{\sigma _k=0}=e^{c_kW_k(U)};p_{\sigma _k=1}=1e^{c_kW_k(U)}$$
(7)
This means that our algorithm has no probability bound violations, which plagued previous attempts to construct an efficient stochastic algorithm .

One can easily see why the introduction of auxiliary $`\sigma `$ fields can be useful. Starting from the current $`\{U_1,\sigma \}`$ configuration, a candidate configuration $`\{U_2\}`$ distributed with the weight $`P_0[U_2]`$ is proposed, and accepted with probability

$$P_{acc}=\mathrm{min}(1,\frac{\mathrm{P}_1[\mathrm{U}_2,\sigma ]}{\mathrm{P}_1[\mathrm{U}_1,\sigma ]})=\mathrm{min}(1,\underset{\mathrm{k}:\sigma _\mathrm{k}=1}{}\frac{\mathrm{e}^{\mathrm{c}_\mathrm{k}\mathrm{W}_\mathrm{k}(\mathrm{U}_2)}1}{\mathrm{e}^{\mathrm{c}_\mathrm{k}\mathrm{W}_\mathrm{k}(\mathrm{U}_1)}1})$$
(8)
Since the terms $`c_kW_k(U)`$ contribute in $`P_{acc}`$ only if $`\sigma _k=1`$, the amount of computational work is greatly reduced if the configurations with $`\sigma _k=0`$ are dominating. That is certainly the case when the absolute values of the coupling coefficients $`|c_k|`$ are small: the probabilities for $`\sigma _k`$ to be unity, averaged over $`\{U\}`$ configurations, are negligible then. Indeed, to leading order in $`c_k`$ the average probability $`p_{\sigma _k=1}`$ from eq.(7) can be written as

$$p_{\sigma _k=1}c_kW_k(U)\mathrm{\hspace{0.33em}0}\mathrm{𝑖𝑓}c_k0$$
(9)
Expression (9) also suggests that one should try to make $`|W_k|`$ as small as possible, using the freedom one has to shift $`W_k`$ by a constant. This goal should remain compatible however with inequalities (3); otherwise, probability bound violations will appear for $`p_{\sigma _k}`$ and $`P_{acc}`$ in eqs.(7,8).

Actually, the violation of conditions (3) is not completely forbidden. As it was pointed in Ref. , one can address the problem of the lower probability-bound violations by redefining the measure. If the distribution $`P_1[U,\sigma ]`$ in (5) can be negative for some configurations $`\{U,\sigma \}`$, one can effectively simulate with the probability distribution $`P_0[U]\left|P_1[U,\sigma ]\right|`$ instead and include the sign $`sgn(P_1)`$ into the observable expectation value:

$$O=\frac{Osgn(P_1)_{||}}{sgn(P_1)_{||}}$$
(10)
where by $`_{||}`$ we denote the averages with respect to distribution $`P_0[U]\left|P_1[U,\sigma ]\right|`$. Sometimes the admission of very rare sign violations can substantially decrease the probability $`p_{\sigma _k=1}`$. However one should be very careful in using this trick: as the volume of the system increases, one needs an exponential growth of statistics to estimate $`sgn(P_1)_{||}`$ within the same accuracy. In the following we shall always assume fulfillment of the inequalities (3).

After updating the $`U`$ fields one should also update the $`\sigma `$ fields to preserve ergodicity. This requires the calculation of probabilities (7). At this point the reader might say: ”OK, one saves computational effort by not calculating some terms $`W_k`$ in expression (8) while estimating $`P_{acc}`$. Nevertheless, one must calculate these terms when updating the $`\sigma `$ fields! So does one gain anything in the end?” The answer is ’yes’ for the following two reasons:

* The terms $`W_k`$ for which it is reasonable to use the stochastic estimation usually couple many degrees of freedom (this is due to the usual nonlocality of weakly coupled terms, which serve as corrections to more local leading terms in the Hamiltonian). If one uses usual local algorithms (without introducing stochastic $`\sigma `$ variables), one should estimate the term $`W_k`$ each time one updates a degree of freedom which it couples. Contrary to that, if one uses Noisy Monte Carlo, the probabilities (7) should be calculated only once per $`\sigma `$ update.
* The variables $`\sigma _k`$ can be refreshed infrequently, the more so as the associated coupling $`c_k`$ gets smaller. It will be demonstrated in the next Section on a particular example. This slow dynamics of the auxiliary $`\sigma `$ fields does not imply slow dynamics of the physically relevant $`U`$ fields.

Up to now we were quite generic, showing that the Noisy Monte Carlo (NMC) method can be potentially very effective for the variety of theories where the energy (2) is a sum of couplings of decreasing strength. In the next section we illustrate these ideas on a particular example: a 5-loop perturbatively improved SU(2) Yang-Mills model.

## III 5-loop SU(2) gauge theory

We consider a 5-loop SU(2) gauge action in $`4d`$:

$`S`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{5}{}}}c_i{\displaystyle \frac{1}{m_i^2n_i^2}}S_{m_i,n_i}`$ (11)
where the indices $`(m_i,n_i)=(1,1),(2,2),(1,2),(1,3),(3,3)`$ for $`i=1,\mathrm{},5`$ denote the planar, fundamental loops of size $`m\times n`$

$`S_{m_i,n_i}`$ $`=`$ $`{\displaystyle \underset{x,\mu ,\nu }{}}\left(2\mathrm{sgn}(\mathrm{c}_\mathrm{i}){\displaystyle \frac{\mathrm{Tr}}{2}}(\text{ }\text{}+\text{ }\text{})\right)`$ (12)
The Gibbs factor is $`\mathrm{exp}(\frac{\beta }{2}S)`$. Note that in eq.(12) we have arranged the constant term $`2sgn(c_i)`$ to ensure the condition (3) for elementary action terms corresponding to each loop:

$`\{i,\mu ,\nu ,x,U\}:S_{i,\mu ,\nu ,x}{\displaystyle \frac{1}{m_i^2n_i^2}}c_i\left(\mathrm{sgn}(\mathrm{c}_\mathrm{i}){\displaystyle \frac{\mathrm{Tr}}{2}}(\text{ }\text{})\right)0`$ (13)
Using the results of  one can construct a one-parameter set of actions which have no $`𝒪(a^2)`$ and $`𝒪(a^4)`$ corrections

$`c_1`$ $`=`$ $`(1955c_5)/9,c_2=(164c_5)/9`$ (14)
$`c_3`$ $`=`$ $`(64+640c_5)/45,c_4=1/52c_5`$ (15)
Here we take $`c_5=1/20`$ (the same action was used in the context of improved cooling ).

Following the ideas of section II, we estimate the contribution of all loops except the plaquette stochastically. For each loop $`l\{\mu ,\nu ,x\}`$ of sort $`2i5`$ we introduce the auxiliary variable $`\sigma _i(l)=0,1`$; and rewrite the contribution of this loop to Gibbs factor in the form

$$e^{\frac{\beta }{2}S_{i,\mu ,\nu ,x}}=\underset{\sigma _i(l)=0,1}{}[\delta _{\sigma _i(l),0}+\delta _{\sigma _i(l),1}(e^{\frac{\beta }{2}S_{i,\mu ,\nu ,x}}1)]$$
(16)
The resulting distribution of $`\{U,\sigma \}`$ fields is used for the generation of independent $`\{U\}`$ configurations. We shall say that for a given $`\{\sigma \}`$ configuration the loop $`\{l,i\}`$ is “active” if $`\sigma _i(l)=1`$.

Let us describe the updating procedure in the $`\{U,\sigma \}`$ configuration space. Consider first the local updating of the gauge fields $`U`$ when the $`\sigma `$ fields are fixed. The proposal value $`U_{x,\mu }^{new}`$ at a given link $`\{x,\mu \}`$ is generated by heatbath with respect to the measure

$$P_0[U]\mathrm{exp}(\frac{\beta }{2}c_1S_{1,1}[U])$$
(17)
where $`S_{1,1}`$ is the plaquette action (see (12)), and then accepted with probability

$$P_{acc}=\mathrm{min}(1,\frac{\mathrm{P}_1[\mathrm{U}_{\mathrm{x},\mu }^{\mathrm{new}},\sigma ]}{\mathrm{P}_1[\mathrm{U}_{\mathrm{x},\mu }^{\mathrm{old}},\sigma ]})$$
(18)
where

$$P_1[U,\sigma ]=\underset{\begin{array}{c}(l,i)\{x,\mu \}\\ \sigma _i(l)=1\end{array}}{}(e^{\frac{\beta }{2}S_{i,l}[U]}1)$$
(19)
Only active loops which contain the given link $`\{x,\mu \}`$ contribute to the expression in the r.h.s. of (18).

After each $`N_i`$ updates of fields $`U`$ on the entire lattice we update the $`\sigma `$ fields of sort $`i`$. For each loop $`l`$ we assign the values 0,1 to the variable $`\sigma _i(l)`$ with the following probabilities

$$p_{\sigma _i(l)=0}=\mathrm{exp}(\frac{\beta }{2}S_{i,l}[U]);p_{\sigma _i(l)=1}=1\mathrm{exp}(\frac{\beta }{2}S_{i,l}[U])$$
(20)
Due to the absence of interaction between different $`\sigma `$ variables, the probabilities (20) depend only on the gauge configuration, so that $`\sigma `$ variables can be updated independently.

In our simulations we have measured the average values of $`\sigma _i,\mathrm{\hspace{0.33em}2}i5`$ which are listed in Table I. They are quite small, and very close to the perturbative estimate (9). This shows that one can avoid the computation of almost all of the extended ”staples” in the $`U`$ update.

| loop | 1x2 | 1x3 | 2x2 | 3x3 |
| --- | --- | --- | --- | --- |
| $`\sigma `$ | 0.0753 | 0.0199 | 0.0202 | 0.0018 |

Table I: Average value of $`\sigma _i`$ field for each loop of sort $`i`$.

Performing numerical simulations for the 5-loop model (11,15) with auxiliary $`\sigma `$ fields, we were mainly interested in the efficiency of our new NMC algorithm. One can estimate the efficiency of NMC by comparing it with the updating procedures which are commonly used now for the simulation of multiloop actions like (11). In the following we label these usually applied techniques with the collective name ”Usual Monte Carlo” (UMC), to contrast with NMC.

We compare the computer times needed to get the same results with NMC and UMC algorithms as follows. First, we make an analytic estimation of the total computational cost of one update of the $`U`$ fields for both algorithms in units of matrix (link) multiplications. Second, we extract from numerical simulations the integrated autocorrelation times for different observables, in units of $`U`$ update. The computer time needed to estimate any given observable is proportional to the product of the computational cost per update and the autocorrelation time.

For NMC the average computational cost of one update of the $`U`$ fields on the entire lattice is equal to

$$t_U^{NMC}=t_U^{pl}+4V\underset{i=2}{\overset{5}{}}n_{staple}(i)n_{mult}(i)\sigma _i$$
(21)
where $`t_U^{pl}`$ is the cost for generating the proposal configuration with measure (17) (i.e. the update cost for the elementary plaquette action), $`4V`$ is the number of links on the lattice, $`n_{staple}(i)`$ is the number of “staples” which the loops of sort $`i`$ form for each link, $`n_{mult}(i)`$ is the number of matrix multiplications needed to estimate the contribution of one staple of sort $`i`$, and the factor $`\sigma _i`$ accounts for the fact that one needs to calculate the contribution of active loops only. One can easily check that

$$n_{staple}(i)=\frac{3}{2}P_is_i;n_{mult}(i)=P_i$$
(22)
where $`P_i2(m_i+n_i)`$ is the perimeter of loop $`i`$, and $`s_i`$ is a symmetry factor: $`s_i=1`$ for square loops and $`s_i=2`$ for rectangular loops. Then we have

$$t_U^{NMC}=t_U^{pl}+6V\underset{i=2}{\overset{5}{}}P_i^2s_i\sigma _i$$
(23)
On the other hand, the computational cost of one update of the $`\sigma _i`$ fields on the entire lattice is given by

$$t_{\sigma _i}=6Vs_iP_i$$
(24)
Here $`6Vs_i`$ is the number of loops of a given sort on the lattice, and the perimeter $`P_i`$ of the loop appears again as the number of matrix multiplications $`n_{mult}(i)`$ needed to calculate the probabilities (20). Since we update the $`\sigma _i`$ fields only once per each $`N_i`$ updates of the $`U`$ fields, the total computational cost per $`U`$ update for the NMC method is

$$t_{tot}^{NMC}=t_U^{NMC}+\underset{i=2}{\overset{5}{}}\frac{t_{\sigma _i}}{N_i}=t_U^{pl}+6V\underset{i=2}{\overset{5}{}}P_is_i(\frac{1}{N_i}+P_i\sigma _i)$$
(25)
Let us note that the computational cost $`t_{tot}^{UMC}`$ for UMC of one $`U`$ update is approximately equal to the r.h.s. of expression (25) in the limit $`N_i\mathrm{}`$ and $`\sigma _i1`$:

$$t_{tot}^{UMC}=t_U^{pl}+6V\underset{i=2}{\overset{5}{}}P_i^2s_i$$
(26)
Indeed, in the limit when all $`\sigma `$ are set equal to 1 and not updated we recover the usual algorithm (certainly one should correct the expressions (18,19) for $`P_{acc}`$ in this case).

Now we can compare the performance of our NMC algorithm with that of UMC. The naive gain in efficiency from using NMC does not depend on the observable measured, and is equal to the ratio between the computational costs (26) and (25):

$$r_{gain}^{naive}\frac{t_{tot}^{UMC}}{t_{tot}^{NMC}}=\frac{t_U^{pl}+6V_{i=2}^5P_i^2s_i}{t_U^{pl}+6V_{i=2}^5P_is_i(\frac{1}{N_i}+P_i\sigma _i)}$$
(27)
Now one should also take into account the increase of autocorrelation times coming from the introduction of auxiliary variables $`\sigma `$ in the NMC algorithm, so that the real gain is

$$r_{gain}^{real}r_{gain}^{naive}\frac{\tau _{int}^{UMC}}{\tau _{int}^{NMC}}$$
(28)
where $`\tau _{int}^{UMC}`$ and $`\tau _{int}^{NMC}`$ are integrated autocorrelation times for UMC and NMC respectively. Note that $`\tau _{int}^{NMC}`$ is a function of the updating frequencies $`1/N_i`$ of the $`\sigma _i`$ fields. Like $`\tau _{int}`$, the ratio (28) will also depend on the observable.

In Table II we present the autocorrelation times for averaged traces of 6 different loops in units of $`U`$ updates. In the first column we show the results for the Usual Monte Carlo, and in other columns for the NMC algorithm with different frequencies of $`\sigma `$ updates for 1x2, 1x3, 2x2, 3x3 loops. In the last row we present the naive gain (27).

Let us make one useful remark. It is not necessary to keep the same updating frequencies $`1/N_i`$ for all sorts $`i`$ of loops. Actually it is even impractical. The computational cost of $`U`$ update coming from the loop of sort $`i`$ is proportional to the average value of $`\sigma _i`$, which is in turn proportional to the coupling (9). As the coupling decreases, we should expect a reduction of the computational effort for the corresponding terms in the action. That is not the case for the cost of $`\sigma _i`$ update: it does not depend on the coupling and even increases with the nonlocality of the action term (factor $`P_i`$ in expression (24)). In order for the work in the $`\sigma `$ and in the $`U`$ updates coming from loops of sort $`i`$ to remain comparable, one should keep the updating frequencies $`1/N_i`$ proportional to $`\sigma _i`$:

$$\frac{1}{N_i}P_i\sigma _i$$
(29)
Due to the small influence of weakly coupled terms on the dynamics of the system, one can expect only insignificant changes in the autocorrelation behavior as $`N_i`$ increases. These considerations are distinctly demonstrated in Table II, where in two columns we present the results for updating frequencies of $`\sigma `$ fields varying in accordance with (29).

| number of | UMC | 1 | 5 | 5 for 1x2 | 10 | 10 for 1x2 | 20 | 30 | 40 | 50 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`U`$ updates per | no | for | for | 15 for 1x3,2x2 | for | 30 for 1x3,2x2 | for | for | for | for |
| 1 $`\sigma `$ update | $`\sigma `$ | all | all | 105 for 3x3 | all | 210 for 3x3 | all | all | all | all |
| $`\tau _{int}`$ (1x1) | 0.7(1) | 1.9(2) | 2.3(1) | 2.5(2) | 3.1(2) | 3.2(2) | 4.3(4) | 4.5(4) | 3.8(2) | 5.1(4) |
| $`\tau _{int}`$ (1x2) | 0.8(1) | 2.6(3) | 2.8(2) | 3.2(2) | 4.3(4) | 3.9(3) | 5.2(4) | 5.6(5) | 5.7(4) | 7.7(8) |
| $`\tau _{int}`$ (1x3) | 0.8(1) | 2.7(3) | 2.8(2) | 3.2(2) | 4.3(4) | 3.9(3) | 5.1(4) | 5.4(5) | 5.4(4) | 7.4(8) |
| $`\tau _{int}`$ (2x2) | 1.0(1) | 3.4(5) | 3.3(3) | 3.7(3) | 4.7(4) | 4.7(3) | 5.3(4) | 5.8(6) | 5.7(4) | 7.5(8) |
| $`\tau _{int}`$ (2x3) | 1.4(3) | 4.2(7) | 3.8(3) | 4.2(3) | 5.5(5) | 5.7(4) | 5.7(5) | 6.2(6) | 6.3(5) | 7.7(8) |
| $`\tau _{int}`$ (3x3) | 1.8(4) | 5.0(8) | 4.5(5) | 5.0(5) | 5.8(6) | 6.4(5) | 5.9(5) | 6.3(6) | 6.3(5) | 7.5(8) |
| $`r_{gain}^{naive}`$ | 1 | 7.2 | 14.8 | 18.3 | 17.9 | 20.6 | 18.3 | 20.8 | 21.1 | 21.4 |

Table II: Integrated autocorrelation times for average loop traces in units of $`U`$ updates for UMC algorithm (first column) and for NMC algorithm with different frequencies of $`\sigma `$ updates for 1x2, 1x3, 2x2 and 3x3 loops (other columns). The last row presents the naive gain for the NMC algorithm (27).

Table II gives an impressive demonstration of the benefits which come from using the NMC algorithm. The naive gain increases substantially as we decrease the frequencies of $`\sigma `$ updates, while the autocorrelation times grow rather slowly. That is particularly visible for the runs where the updating frequencies for $`\sigma `$ fields are adjusted as per eq.(29). For such runs we can infer that the ’real gain’ $`𝒪(46)`$ in computer time (28) for the observables measured is large enough for a convincing demonstration of the possible advantages coming from using the NMC algorithm.

Let us make a conclusion for this section. We have applied our NMC algorithm for the 5-loop model (11,15). We have shown that with this algorithm a significant gain in efficiency is obtained in comparison with usual updating techniques. Finally we note that the action (11) is a relatively simple one, and one can expect a much greater gain for more complicated highly-improved actions with many nonlocal weakly coupled terms.

## IV Discussion

Let us summarize our algorithm:
$`a)`$ Separate the Hamiltonian (or action) into a dominant term $`c_0W_0`$, to be calculated exactly, and correction terms $`_{k=1}^mc_kW_k`$, to be estimated stochastically.
$`b)`$ Shift the correction terms to guarantee $`c_kW_k0`$.
$`c)`$ Introduce auxiliary local variables $`\sigma _k(l)`$, through identity (6) (here $`l`$ runs through all the elementary ’bonds’ which form $`W_k`$, e.g. loops in gauge theory).
$`d)`$ Update the auxiliary variables $`\sigma _k`$ by heatbath.
$`e)`$ To update the original variables $`U`$, propose a new value $`U^{}`$ sampled from the distribution $`e^{c_0W_0}`$, and accept it with the Metropolis probability $`\mathrm{min}(1,_{\mathrm{k}1;\sigma _\mathrm{k}=1}(\mathrm{e}^{\mathrm{c}_\mathrm{k}\mathrm{W}_\mathrm{k}(\mathrm{U}^{})}1)/(\mathrm{e}^{\mathrm{c}_\mathrm{k}\mathrm{W}_\mathrm{k}(\mathrm{U})}1))`$.

The essential advantage of our algorithm appears in step $`(e)`$: only the terms $`W_k`$ whose associated $`\sigma _k`$ is equal to 1 need to be computed. Since on average, $`\sigma _k`$ goes to zero with $`c_k`$, the computation of almost all correction terms can be avoided.

To avoid simply shifting the cost of the algorithm to step $`(d)`$, we propose to refresh the variables $`\sigma _k`$ infrequently, the more so as the associated coupling $`|c_k|`$ gets smaller. We have pointed out that this introduction of slow dynamics for the $`\sigma _k`$ does not enforce slow dynamics for the system, since $`W_k(l)`$ will fluctuate regardless of the value of $`\sigma _k(l)`$. Our numerical study of Section III confirms this statement.

Let us now speculate on possibilities to use our algorithm to simulate a Hamiltonian with a very large number of terms. A specific example we have in mind is the case of full QCD, where the measure is, for 2 flavors of Wilson quarks:

$$\frac{1}{Z}e^{S_g(U)}det^2(\mathrm{𝟏}\kappa M(U))$$
(30)
where $`S_g`$ is the local gauge action, $`M(U)`$ is a hopping matrix connecting nearest neighbours on a $`4d`$ hypercubic grid, and $`Z`$ normalizes the distribution. The determinant can be turned into $`exp(Tr(Log(\mathrm{𝟏}\kappa M(U))))`$, then the logarithm expanded around 1, giving the loop expansion of the measure above:

$$\frac{1}{Z}e^{S_g(U)2_{l=4}^{\mathrm{}}\frac{\kappa ^l}{l}TrM(U)^l}.$$
(31)
$`TrM(U)^l`$ can be represented as a sum over all closed non-backtracking loops of length $`l`$ on the $`4d`$ hypercubic lattice. The number of types of contributing loops $`n_l`$ is bounded by $`(2d1)^l=7^l`$, because of the branching factor at each hop. Although this upper bound is not saturated, it is clear that the multiplicity of terms of a given length $`l`$ grows exponentially:

$$n_lF_1(l)\alpha ^l;\alpha <7$$
(32)
where $`F_1(l)`$ is a rational function of $`l`$, and $`\alpha ^l`$ is the leading exponential ascend of the number of loops of length $`l`$ in the limit of large $`l`$. For this reason, it seems that sampling numerically the distribution (31) is a disastrous idea: the action contains an infinite number of terms, of exponentially growing multiplicity. Instead, other strategies are being used, based on the transformation of the determinant (30) into a Gaussian integral.

Nevertheless, the coupling $`\frac{\kappa ^l}{l}`$ decreases exponentially as $`l`$ increases. Therefore, the auxiliary variables $`\sigma _l`$ associated in our approach with various loops of length $`l`$ will take value 0 almost always. From (9) one gets

$$\sigma _lc_lF_2(l)k^l\beta ^l$$
(33)
where the exponentially growing factor $`\beta ^l`$ comes from the average trace of Dirac matrices along the loops of length $`l`$, $`F_2(l)`$ is again the rational function.

If one arranges the updating frequencies for each loop $`i`$ as per eq.(29), one can expect that the average computer time needed for estimation of contribution of all loops of length $`l`$ behaves at large $`l`$ as

$$t_ln_ll^2\sigma _lF(l)(\alpha \beta \kappa )^l$$
(34)
(It was pointed above that one should not expect a significant growth of autocorrelation coming from slow dynamics for $`\sigma _l`$.) For $`\kappa <\kappa _{ca}=\frac{1}{\alpha \beta }`$ the computational cost $`t_l`$ decreases exponentially with $`l`$ and the total computational cost of the algorithm $`t=_lt_l`$ converges to some finite value. Our rough estimation from fitting $`n_l`$ and average trace of Dirac matrices in the interval $`4l12`$ gives $`\alpha 5.4;`$ $`\beta 1.4`$, and therefore $`\kappa _{ca}0.13`$.

In the regime $`\kappa <\kappa _{ca}`$ we are in an interesting situation where the influence of very large loops is negligible because their associated coupling in the effective action is extremely small. Therefore, truncating the loop expansion above a certain order will introduce a statistically unobservable bias. Equivalently, one can freeze the associated $`\sigma `$ variables at the value zero, or update them with arbitrarily low frequency. In spite of this extremely (or infinitely) slow dynamic mode of the $`\sigma `$’s, the dynamics of the gauge fields are not affected. Note that the cost of our algorithm grows linearly with the volume $`V`$ of the system. This is better than alternative approaches to the simulation of full QCD: Hybrid Monte Carlo (cost $`V^{5/4}`$) and MultiBoson (cost $`V(LogV)^2`$) . In addition the stepsize, or typical change at each update of a gauge link $`U`$, does not seem restricted a priori for small quark mass, unlike in the two alternative approaches above. Unfortunately the possible high efficiency of our algorithm is counterweighted by its extreme programming complexity.

A less speculative use of our algorithm for full QCD consists of truncating the loop expansion eq.(31) to some order $`l_{max}`$, and to represent the higher orders with the MultiBoson approach . This strategy, called “UV-filtered MultiBoson”, has already been used successfully . However, in Ref. the loop expansion is truncated to its lowest term $`l=4`$, because the exact evaluation of larger loops is too time-consuming. With our stochastic approach, these larger loops can be estimated at low cost. We expect this composite strategy to be particularly efficient.

Acknowledgments: We acknowledge communications with T. DeGrand, M. Ilgenfritz and U. Wenger. T.B. was supported by INTAS under grant 96-0370 and Russian Basic Research Fund under grant 99-01-00190.