text:
# General existence proof for rest frame systems in asymptotically flat space-time

## 1 Introduction

The rest frame equation, or “nice section” equation, for asymptotically flat space-times, is an equation which relates quantities defined at null infinity (scri or $``$). The solutions of this equation are a family of privileged sections of scri, which represent the analog of having a frame of reference.

The physical motivations for the nice section construction have been explained already in the literature. We just mention here that at present the description of astrophysical gravitational radiation in terms of the sources, is carried out through the so called quadrupole radiation formula. The construction of the nice sections provides a necessary tool, that we expect to be useful for the discussion of gravitational radiation and its relation with the notion of multipole moments, including the cases of systems which suffer a non-trivial back reaction due to gravitational radiation, or when the lapse of time considered is relatively long, when compared to the dynamical characteristic time. Those cases can not be dealt with the quadrupole radiation formula, which has been deduced for an instant of time, and for weak fields.

In what follows we will concentrate in the mathematical aspects of the nice section construction, and we will provide a general existence proof of them; which demonstrates that the conjecture that we have enunciated in ref.  is true.

Let us start by recalling the basic structure of scri. Roughly speaking scri is the three dimensional manifold that corresponds to the end points when one takes the limits to infinity of the future directed null geodesic. The topology of a complete scri is $`\times S^2`$, and the natural coordinates are the so called Bondi coordinates; which are labeled by $`(u,\zeta ,\overline{\zeta })`$, where $`(\zeta ,\overline{\zeta })`$ are complex stereographic coordinates of the sphere $`S^2`$, and $`u`$ takes values in $``$. Given an arbitrary real function $`\gamma (\zeta ,\overline{\zeta })`$ of the sphere, we say that $`u=\gamma (\zeta ,\overline{\zeta })`$ defines a section of $``$.

The symmetry group of scri (i.e. the group of coordinate transformations which preserves the asymptotic geometric structure of the space-time) is the Bondi-Metzner-Sach (BMS) group. It is defined by the following relations:

$$\stackrel{~}{u}=K(u\gamma (\zeta ,\overline{\zeta })),\stackrel{~}{\zeta }=\frac{a\zeta +b}{c\overline{\zeta }+d},$$
(1)
with

$$adbc=1,K=\frac{1+\zeta \overline{\zeta }}{|a\zeta +b|^2+|c\zeta +d|^2}$$

where $`a,b,c,d`$ are complex constant, and $`\gamma `$ is a real function on the sphere.

We express the supermomentum at a section $`u=\text{constant}`$ in the following way:

$$P_{lm}=\frac{1}{\sqrt{4\pi }}Y_{lm}(\zeta ,\overline{\zeta })\mathrm{\Psi }(u,\zeta ,\overline{\zeta })𝑑S^2,$$
(2)
where $`dS^2`$ is the surface element of the unit sphere; $`Y_{lm}`$ are the spherical harmonics;

$$\mathrm{\Psi }=\mathrm{\Psi }_2+\sigma \dot{\overline{\sigma }}+ð^2\overline{\sigma };$$
(3)
with $`\mathrm{\Psi }_2`$ and $`\sigma `$ being the leading order asymptotic behavior of the second Weyl tensor component and the Bondi shear respectively; where we are using the GHP notation for the edth operator of the unit sphere, and an over dot means $`/u`$.

Similarly the total Bondi energy-momentum vector at a given section of null infinity can be expressed by

$$P^a=\frac{1}{4\pi }l^a(\zeta ,\overline{\zeta })\mathrm{\Psi }(u,\zeta ,\overline{\zeta })𝑑S^2,$$
(4)
where

$$(l^a)=(1,\frac{\zeta +\overline{\zeta }}{1+\zeta \overline{\zeta }},\frac{\zeta \overline{\zeta }}{i(1+\zeta \overline{\zeta )}},\frac{\zeta \overline{\zeta }1}{1+\zeta \overline{\zeta }}).$$
(5)
In fact the Bondi momentum is a linear combination of the first four component $`(l=0,1)`$ of the supermomentum.

The rest mass $`M`$ of the section is given by

$$M=\sqrt{P^aP_a},$$
(6)
where the indices are moved with the Lorenzian flat metric $`\eta _{ab}`$ at scri.

The BMS group (1) and the supermomentum (2) are the basic ingredients for the definition of the nice section equation.

A rest frame system in special relativity is a coordinate system in which the total four momentum has only temporal component. Given such a system we can perform a translation and we still have a rest frame system.

This idea can be translated to the structure of scri; where the supermomentum takes the role of the momentum in special relativity. With this idea in mind, the nice sections are defined in such a way that the supermomentum has only temporal component. This means that the integrand $`\mathrm{\Psi }`$ at a nice section is a constant, i.e., it does not depend on the angular coordinates $`(\zeta ,\overline{\zeta })`$.

The function $`\mathrm{\Psi }`$ transforms under a BMS transformation in the following way

$$\stackrel{~}{\mathrm{\Psi }}=\frac{1}{K^3}(\mathrm{\Psi }ð^2\overline{ð}^2\gamma ),$$
(7)
where $`\stackrel{~}{\mathrm{\Psi }}=\stackrel{~}{\mathrm{\Psi }}_2+\stackrel{~}{\sigma }\dot{\overline{\stackrel{~}{\sigma }}}+\stackrel{~}{ð}^2\overline{\stackrel{~}{\sigma }}`$.

This means that given an initial arbitrary section, a nice section can be determined by the appropriate BMS transformation such that $`\stackrel{~}{\mathrm{\Psi }}`$ is constant. Therefore we can express the nice section equation by

$$ð^2\overline{ð}^2\gamma =\mathrm{\Psi }(u=\gamma ,\zeta ,\overline{\zeta })+K^3(\gamma ,\zeta ,\overline{\zeta })M(\gamma ),$$
(8)
where the conformal factor $`K`$ can be related to the Bondi momentum by

$$K=\frac{M}{P^al_a},$$
(9)
and $`P^a`$ is evaluated at the section $`u=\gamma `$; which is calculated through the integral (4).

Since the operator $`ð^2\overline{ð}^2`$ occurs very frequently in this work we use the notation $`D=ð^2\overline{ð}^2`$, and we note also that

$$Dð^2\overline{ð}^2=\frac{1}{4}\mathrm{\Delta }^2\frac{1}{2}\mathrm{\Delta }$$
(10)
where $`\mathrm{\Delta }`$ is the Laplacian of the unit sphere.

We say that a supertranslation $`\gamma `$ (or equivalently a section of $``$) satisfies the “nice section” equation if the function $`\stackrel{~}{\mathrm{\Psi }}`$ evaluated in $`\gamma `$ is equal to the total mass $`M`$; this is exactly the content of equation (8).

Note that the nice section equation involves also a Lorenz rotation, as it can be seen by the appearance of the scalar $`K`$ in the equation.

Equation (8) is our fundamental equation, we want to investigate the existence and properties of solutions of this equation.

We say that a real function $`x=x(\zeta ,\overline{\zeta })`$ on the sphere is a translation if $`Dx=0`$. Note that this implies that $`x`$ has an expansion in spherical harmonics with $`l=0,1`$. An arbitrary regular function $`\gamma `$ can be decomposed in

$$\gamma =x+y,$$

where $`y`$ has an expansion with $`l2`$. Given an arbitrary $`x`$, equation (8) is an equation for $`y`$.

The main result of this paper is the following Theorem:

###### Theorem 1.1

If $`\mathrm{\Psi }`$ is a smooth function on scri, the total energy $`P^0`$ is bounded by the constant $`E_0`$, the total mass $`M`$ is bounded from below by $`MM_0`$, $`M_0>0`$; and the gravitational energy density flux $`|\dot{\sigma }|^2\lambda `$; where the constant $`\lambda `$ satisfies

$$\lambda <\frac{\sqrt{27}}{4}(1+(2C_K)^4)^1,$$

and $`C_K`$ is given by

$$C_K=\frac{E_0}{M_0}+\sqrt{\frac{E_0^2}{M_0^2}1};$$

then

(i) For every translation $`x`$ there exists a solution $`y`$ of equation (8), and $`y`$ is a smooth function on the sphere.

(ii) The solutions $`\gamma =x+y(x)`$ are continuous in the 4-parameter translation $`x`$, and if $`x_1`$ and $`x_2`$ are two translations such that the difference, $`x_2x_1`$, corresponds to a future directed time like vector, then

$$\gamma (x_1,\zeta ,\overline{\zeta })<\gamma (x_2,\zeta ,\overline{\zeta }).$$

Since the mass $`M`$ is a decreasing function of $`u`$, the constant $`M_0`$ is the final rest mass as $`u\mathrm{}`$. Note that while $`M_0`$ is Lorentz invariant $`P^0`$ is not; in other words the inequalities depend on the particular Bondi system being used.

It is important to recall that all the hypothesis of the theorem are in terms of physical quantities; namely the total mass $`M`$, the energy $`P^0`$ and the gravitational energy density flux $`|\dot{\sigma }|^2`$. Theorem 1.1 essentially says that the solutions of equation (8) exist and have the expected physical properties, when the gravitational radiation of the space-time is not too high.

As it was pointed out in , there exist another choices for the supermomentum $`\mathrm{\Psi }`$ in the literature. We want to analyze here these other possibilities. One of them is the Geroch supermomentum

$$\mathrm{\Psi }_G=\mathrm{\Psi }_2+\sigma \dot{\overline{\sigma }}+\frac{1}{2}(ð^2\overline{\sigma }\overline{ð}^2\sigma );$$

and another one is the Geroch-Winicour supermomentum111We take this opportunity to correct a mistake in the articles  and , the supermomentum $`\mathrm{\Psi }`$ *is not* the same as the Geroch-Winicour supermomentum, as it is clear from the expression (11)

$$\mathrm{\Psi }_{GW}=\mathrm{\Psi }_2+\sigma \dot{\overline{\sigma }}\overline{ð}^2\sigma .$$
(11)
The function $`\mathrm{\Psi }_G`$ is invariant under supertranslation, then it will not give an equation for $`\gamma `$. Note that $`\mathrm{\Psi }_G`$ vanishes for Minkowski space-time, but equation (8) is not trivial in this space-time; its solutions are the shear-free sections ($`\sigma =0`$). As it is well known, the sections with $`\sigma =0`$ in Minkowski space-time represent a natural way to isolate the Poincare group from the BMS group of scri, and this is precisely the aim of the ‘nice section’ equation in the general case. These considerations rule out $`\mathrm{\Psi }_G`$.

The Geroch-Winicour supermomentum $`\mathrm{\Psi }_{GW}`$ could in principle be used to construct an alternative equation. But the supermomentum $`\mathrm{\Psi }`$ satisfies the following remarkable equation

$$\dot{\mathrm{\Psi }}=|\dot{\sigma }|^2,$$
(12)
which relates in a simple way the time derivative of $`\mathrm{\Psi }`$ with the gravitational radiation flux. This equation allows us to write the hypothesis of the theorems in terms only of the total mass $`M`$, the energy $`P^0`$ and the gravitational energy density flux $`|\dot{\sigma }|^2`$. It might be also possible to prove similar theorems for an equation based on $`\mathrm{\Psi }_{GW}`$, but they would certainly look much more complicated.

In section 2 we prove the existence part (i) of the theorem. In section 3 we prove the part (ii). Finally, for the sake of completeness, we remind in the Appendix the Schauder fixed point theorem, which is the main tool for the existence proof of section 2.

## 2 Existence Proof

The equation (8) is an elliptic, nonlinear, equation on the compact manifold $`S^2`$. The main tool for proving existence of solutions of this type of equations is the Schauder fixed point Theorem 3.2 (see for example page 380 of ref.  for an elementary treatment of a similar equation with the same technique).

Let us define the function $`f`$ as the right-hand side of equation (8)

$$f(x+y,\zeta ,\overline{\zeta })=\mathrm{\Psi }(x+y,\zeta ,\overline{\zeta })+K^3M(x+y,\zeta ,\overline{\zeta }).$$
(13)
The first step is to prove that the function $`f`$ satisfies the appropriate inequality under the hypothesis of Theorem 1.1; this will be proved in the following Lemma.

###### Lemma 2.1

Assume that $`\mathrm{\Psi }`$ is a smooth function on scri, the gravitational radiation flux satisfies $`|\dot{\sigma }|^2\lambda `$ , the total mass $`MM_0>0`$ and the total energy $`P^0E_0`$. Then

$$|f|\lambda |y|+C,$$
(14)
where the constant $`C`$ does not depend on $`y`$.

*Proof:* Take $`u_0`$ an arbitrary, fixed, real number. From equation (12) we have

$$\mathrm{\Psi }(\gamma ,\zeta ,\overline{\zeta })=_{u_0}^\gamma |\dot{\sigma }|^2(u^{},\zeta ,\overline{\zeta })𝑑u^{}+\mathrm{\Psi }(u_0,\zeta ,\overline{\zeta }).$$

We use the hypothesis $`|\dot{\sigma }|^2\lambda `$ to obtain

$$|\mathrm{\Psi }|\lambda |y|+C_1,$$
(15)
where

$$C_1=\underset{S^2}{sup}\{\lambda |xu_0|+|\mathrm{\Psi }(u_0,\zeta ,\overline{\zeta })|\}.$$

Note that the constant $`C_1`$ depends on $`x`$ and $`u_0`$ but not on $`y`$.

Looking at the expression (13) we see that for proving the inequality (14) we need only to find bounds for $`K`$ and $`M`$ in terms of the constant $`M_0`$ and $`E_0`$. The last one follows immediately since $`MP^0`$ and by hypothesis we have $`P^0E_0`$, then

$$ME_0.$$
(16)
Define the velocity vector $`V^a`$ by $`V^a=P^a/M`$. Note that $`V^aV_a=1`$. We use the definition of $`K`$, equation (9), to obtain

$$K=\frac{1}{V^al_a}.$$
(17)
The function $`V^al_a`$ satisfies the following inequality

$$|V^al_a||V^0||V^il_i||V^0|\sqrt{V^iV_i}.$$

We have used that $`l^il_i=1`$, with $`i=1,2,3`$. Then we use that $`(V^0)^2V^iV_i=1`$ to obtain

$$|V^al_a||V^0|\sqrt{(V^0)^21}.$$
(18)
The right-hand side of this inequality is a monotonically decreasing function of $`|V^0|`$. By hypothesis we have the bound

$$V^0=\frac{P^0}{M}\frac{E_0}{M_0}.$$
(19)
We use the bound (19) in equation (18) to obtain

$$|V^al_a|\frac{C_0}{M_0}\sqrt{\frac{C_0^2}{M_0^2}1};$$

and then by equation (17) we get the desired bound for $`K`$:

$$|K|C_K,$$
(20)
where $`C_K`$ is given by

$$C_K=\left(\frac{E_0}{M_0}\sqrt{\frac{E_0^2}{M_0^2}1}\right)^1=\frac{E_0}{M_0}+\sqrt{\frac{E_0^2}{M_0^2}1}.$$

Using the bounds (15), (16) and (20) we obtain the final inequality for the function $`f`$:

$$|f|\lambda |y|+C,$$

where

$$C=C_1+E_0C_K^3.$$

$`\mathrm{}`$

With this Lemma we are in position to prove the following existence Theorem:

###### Theorem 2.1

Assume that $`\mathrm{\Psi }`$ is a smooth function on scri, the total energy $`P^0E_0`$, the total mass $`MM_0>0`$ and the gravitational energy density flux $`|\dot{\sigma }|^2\lambda `$ with

$$\lambda <\sqrt{\frac{27}{4}}.$$

Then for every translation $`x`$ there exists a solution $`y`$ of equation (8), and $`y`$ is a smooth function on the sphere.

*Proof:* Consider the mapping $`A`$ on the space of continuous functions on the sphere $`A:C^0(S^2)C^0(S^2)`$ defined as follows

$$A(y)=D^1f(x+y,\zeta ,\overline{\zeta }),$$

where $`D^1`$ is the inverse of $`D`$. First note that $`A`$ is well defined for every translation $`x`$, since for an arbitrary $`x+y`$ the function $`f`$ is orthogonal to the kernel of $`D`$, and then it is in the domain of $`D^1`$. This is not obvious and it is essential, otherwise equation (8) would be inconsistent. For a proof see .

A solution of equation (8) is a fixed point $`y=A(y)`$ of the mapping $`A`$. To prove that such fixed point exists we will use the Schauder fixed point Theorem 3.2. In the following we will prove that the mapping $`A`$ satisfies the hypothesis of this Theorem.

First we prove that $`A`$ is continuous in the $`C^0`$ norm. By definition

$$DA(y)=f(x+y,\zeta ,\overline{\zeta }),$$

we take the $`L^2`$ norm in both sides of this equation

$$DA(y)_{L^2}=f(x+y,\zeta ,\overline{\zeta })_{L^2}\sqrt{4\pi }f(x+y,\zeta ,\overline{\zeta })_{C^0}.$$

The last inequality follows because $`f`$ is continuous. We use the inequality

$$Dy_{L^2}\sqrt{4\pi }\sqrt{27/4}y_{C^0},$$
(21)
which is nothing but a particular case of the elliptic regularity inequalities for the elliptic operator $`D`$, to obtain

$$A(y)_{C^0}\sqrt{4/27}f(x+y,\zeta ,\overline{\zeta })_{C^0}.$$
(22)
Since $`D^1`$ is a linear operator we obtain also the inequality

$$A(y_1)A(y_2)_{C^0}=D^1(f(y_1)f(y_2))_{C^0}\sqrt{\frac{4}{27}}(f(y_1)f(y_2))_{C^0}.$$
(23)
Since by hypothesis $`f`$ is continuous as a function of $`y`$ then by (23) $`A`$ is also continuous.

Let $`B`$ be the closed ball in $`C^0(S^2)`$ of radius $`C_2`$; $`B`$ is clearly a closed, convex subset of the Banach space $`C^0(S^2)`$. We want to prove that $`A`$ maps $`B`$ into itself when the constant $`C_2`$ is appropriately chosen. Take $`yB`$, using the inequality (22) and Lemma 2.1 we have

$$A(y)_{C^0}\sqrt{\frac{4}{27}}(\lambda |y|+C)\sqrt{\frac{4}{27}}(\lambda C_2+C).$$
(24)
Take $`C_2=\alpha C`$, where $`\alpha `$ satisfies

$$\lambda +\frac{1}{\alpha }\sqrt{\frac{27}{4}};$$

this $`\alpha `$ exists because we have assumed that $`\lambda <\sqrt{\frac{27}{4}}`$. Then from the inequality (24) we obtain

$$A(y)_{C^0}C_2.$$

Therefore we have proved that $`A`$ maps the ball $`B`$ into itself.

It remains to be proved that the set $`A(B)`$ is precompact. The set $`A(B)`$ is a subset of $`H^4(S^2)`$. This is a consequence of the elliptic character of the 4th order operator $`D`$ and the standard elliptic regularity theorems (see for example  )

Since $`H^4(S^2)`$ is compactly imbedded in $`C^0(S^2)`$ (see for example ) it follows that $`A(B)`$ is precompact in $`C^0(S^2)`$.

The hypothesis of Theorem 3.2 are satisfied, and then there exists a solution $`y=y(\zeta ,\overline{\zeta })`$ in $`C^0(S^2)`$.

Since $`y`$ is in $`C^0(S^2)`$ and $`f`$ is smooth, then $`f(x+y,\zeta ,\overline{\zeta })`$ is in $`L^2(S^2)`$. Therefore by the elliptic regularity theorems $`yH^4(S^2)`$; and then, by the Sobolev imbedding Theorem, it is in $`C^2(S^2)`$; but then $`f(x+y,\zeta ,\overline{\zeta })`$ is in $`H^2(S^2)`$. We use elliptic regularity and induction to conclude that $`y`$ is in $`C^{\mathrm{}}(S^2)`$. $`\mathrm{}`$

Note that for the existence proof we only need $`\mathrm{\Psi }`$ to be in $`C^0`$; if we only require this, then by the elliptic regularity theory one would only obtain a $`C^2`$ solution $`y`$.

## 3 Uniqueness and non self-crossing properties of the solutions

Here we want to analyze the dependence of the solution $`y(x,\zeta ,\overline{\zeta })`$ in terms of the translation $`x`$. In the previous section we have shown that given an arbitrary $`x`$ there exists a solution $`y(x,\zeta ,\overline{\zeta })`$, which is smooth in the variables $`(\zeta ,\overline{\zeta })`$. Now we want to prove that under more restrictive conditions this solution is unique. Moreover, we also want to prove that $`y(x,\zeta ,\overline{\zeta })`$ is continuous in $`x`$. This means that the function $`y(x,\zeta ,\overline{\zeta })`$ defines a family of solutions parameterized by $`x`$. We also prove that this family has the following non self-crossing property: if $`x_1`$ and $`x_2`$ are two translations such that the difference $`x_2x_1`$ corresponds to a future directed timelike vector, then

$$\gamma (x_1,\zeta ,\overline{\zeta })<\gamma (x_2,\zeta ,\overline{\zeta }).$$

This implies that the solutions of equation (8) generated by a time translation do not cross among them.

We begin with the following auxiliary Lemma.

###### Lemma 3.1

Assume that the hypothesis of Lemma 2.1 holds. Then

$$f(\gamma _2)f(\gamma _1)_{L^2}\sqrt{4\pi }C_f^{}\gamma _2\gamma _1_{C^0},$$

where

$$C_f^{}=\lambda \left(1+(2C_K)^4\right).$$
(25)
*Proof:* Consider $`f`$ as a mapping $`f:C^0(S^2)L^2(S^2)`$, since $`\mathrm{\Psi }`$ is smooth then $`f`$ is Fréchet differentiable, and by the mean value theorem (see for example  ) we obtain the following bound

$$||f(\gamma _2)f(\gamma _1)||_{L^2}\underset{t(0,1)}{sup}||f^{}(\gamma _1+t(\gamma _2\gamma _1)||_{(C^0,L^2)}||\gamma _2\gamma _1||_{C^0},$$

where $`(C^0,L^2)`$ denotes the operator norm, and $`f^{}`$ is the Fréchet derivative of $`f`$. We have to calculate the bound for $`f^{}`$. The Fréchet derivative $`f^{}`$ is given explicitly by

$$f^{}(\gamma )\delta \gamma =\dot{\mathrm{\Psi }}(\gamma )\delta \gamma +\left(\frac{4P^a}{M}3Kl^a\right)K^3\delta P_a,$$

where

$$\delta P_a=\frac{1}{4\pi }l_a\dot{\mathrm{\Psi }}(\gamma )\delta \gamma 𝑑S.$$

And the operator norm of $`f^{}`$ is

$$f^{}(\gamma )_{(C^0,L^2)}=sup\{f^{}(\gamma )\delta \gamma _{L^2};\delta \gamma _{C^0}=1\}.$$

We use the bounds $`E_0,M_0,\lambda `$ and the bound (20) for $`K`$ to obtain

$$f^{}(\gamma )_{(C^0,L^2)}\sqrt{4\pi }C_f^{},$$

where

$$C_f^{}=\lambda \left(1+(2C_K)^4\right),$$

and we have used the inequality

$$\frac{1}{C_K}KC_K.$$

$`\mathrm{}`$

###### Theorem 3.1

Assume that the hypotheses of Theorem 2.1 hold, and in addition that

$$\lambda <\frac{\sqrt{27}}{4}(1+(2C_K)^4)^1.$$
(26)
Then the solutions $`\gamma =x+y(x)`$ are continuous in the 4-parameter translation $`x`$. Moreover, if $`x_1`$ and $`x_2`$ are two translations such the difference $`x_2x_1`$ corresponds to a future directed time like vector, equivalently $`x_2x_1>0`$; then

$$\gamma (x_1,\zeta ,\overline{\zeta })<\gamma (x_2,\zeta ,\overline{\zeta }),$$

for all $`(\zeta ,\overline{\zeta })`$ in the sphere.

*Proof:* Let $`x_1+y(x_1,\zeta ,\overline{\zeta })`$ and $`x_2+y(x_2,\zeta ,\overline{\zeta })`$ be two solutions of the equation (8); taking the difference we obtain

$$D(y(x_2)y(x_1))=f(x_2+y(x_2),\zeta ,\overline{\zeta })f(x_1+y(x_1),\zeta ,\overline{\zeta }).$$

Let us take the $`L^2`$ norm on both sides of this equation

$$D(y(x_2)y(x_1))_{L^2}=f(x_2+y(x_2),\zeta ,\overline{\zeta })f(x_1+y(x_1),\zeta ,\overline{\zeta })_{L^2},$$

then using Lemma 3.1 and the inequality (21) we obtain

$$\sqrt{27/4}y(x_2)y(x_1)_{C^0}C_f^{}x_2+y(x_2)(x_1+y(x_1))_{C^0}.$$

Since by hypothesis we have inequality (26), it is deduced that $`C_f^{}<\sqrt{27}/4`$ and we obtain

$$y(x_2)y(x_1)_{C^0}\frac{C_f^{}}{\sqrt{27/4}C_f^{}}x_2x_1_{C^0};$$

in this way we have proved that the function $`y(x)`$ is continuous in $`x`$. Note that

$$\frac{C_f^{}}{\sqrt{27/4}C_f^{}}<1,$$

then for $`x_1<x_2`$ we obtain

$$y(x_2)y(x_1)_{C^0}<x_2x_1_{C^0}.$$
(27)
The fact that $`x_2x_1`$ corresponds to a future directed timelike vector, is equivalent to the statement that we can choose a coordinates system (by means of a Lorentz rotation) in which $`x_2x_1`$ is a positive constant. Then $`x_2x_1_{C^0}=x_2x_1`$ in this coordinate system and from the inequality (27) we conclude that

$$\gamma (x_2,\zeta ,\overline{\zeta })>\gamma (x_1,\zeta ,\overline{\zeta }).$$

It is important to remark that after one has proved the last inequality, it is possible to come back to the original Bondi system by the inverse of the original transformation which does not change the inequality. That is, this inequality is also true in the original Bondi system. $`\mathrm{}`$

## Appendix: Schauder Fixed Point Theorem

We use the following version of the Schauder Theorem (see ref.  page 280.)

###### Theorem 3.2

Let $`B`$ be a closed convex set in a Banach space $`V`$ and let $`A`$ be a continuous mapping of $`B`$ into itself such that the image $`A(B)`$ is precompact (i.e. the closure of $`A(B)`$ is compact). Then $`A`$ has a fixed point.

## Acknowledgments

We are deeply indebted to G. Nagy and O. Reula for very illuminating discussion about the Schauder fixed point theorem.

We acknowledge support from Fundación Antorchas, SeCyT-UNC, CONICET and FONCYT BID 802/OC-AR PICT: 00223.


format: text

text:
# Complexity transitions in global algorithms for sparse linear systems over finite fields

## I Introduction

The methods and concepts of statistical physics of disordered systems constitute a very useful tool for the understanding of the onset of computational complexity in randomly generated hard combinatorial problems. Once the optimization problems are translated into zero temperature spin glass problems, one may study the geometrical changes in the space of solutions as symmetry breaking phenomena. In this context one may view the exponential regimes of randomized search algorithms as out-of-equilibrium phases of stochastic processes.

However, combinatorial problems are not always exponentially hard: Problems that can be solved in polynomial time, even in their worst-case realizations compose the so called Polynomial (P) class GAREY . Such problems are often of great practical relevance and are tackled using large scale computations. Examples can be found in all disciplines: In physics, just to make one example, one may study ground states of 2D spin glass like Hamiltonians resorting to a polynomial max-cut algorithm ALAVA . The major application are obviously found in engineering: Examples are design problems (finite elements methods), control theory (convex optimization), coding theory (parity check equations) and cryptography (integer factorization).

Due to the practical relevance of the problems and to the typically large number of variables used for their encoding, that is the size of the problems, it is of basic interest to look at the fine structure of the class P in order to concretely optimize the computational strategies. For instance, in error correcting codes it is crucial to have algorithms that converge in linear time with respect to the number of encoded bits, any power larger than one being considered of no practical interest.

Quite in general, the trade-off between time and memory resources is the guiding criterion which selects the algorithms used in real-world applications. Roughly speaking polynomial algorithms can be divided in different groups depending on the solving strategy they implement. The main groups are local algorithms (e.g. greedy/gradient methods), global algorithms (e.g. Gaussian elimination or Fourier transforms methods), iterative algorithms (e.g. Lanczos method) and parallel algorithms. See Ref. CLR  for a basic introduction to the subject.

In what follows we shall study a prototype problem of the P class, that is the problem of solving large and random sparse systems in some Galois field GF($`q`$). Working in GF($`q`$) is completely equivalent to perform any operation modulo $`q`$.

Firstly, we give a precise analysis of the computational features for non-trivial ensembles of random instances. By a statistical mechanics study, we look into the – symmetry breaking – geometrical structure of the space of solution thereby providing an explanation for the changes in the power law behavior observed in different algorithms. Moreover, we are able to predict and explain in terms of clustering of solutions, the memory catastrophe found in global algorithms such as Gaussian elimination. Such an effect seriously hampers application of this sort of global algorithms in many circumstances, one example being symbolic manipulations. This memory catastrophe induce in turn an even more dramatic increase in CPU time, which make large problems unaffordable above the dynamical threshold $`\gamma _d`$ (see below for its definition).

Secondly, we consider a specific “real-world” application, namely the Integer Factorization problem used in RSA public key cryptography RSA . By a non-trivial mapping of the factoring problem on a sparse linear system modulo 2, endowed with a quite peculiar statistical distribution of matrix elements, we analyze which are the characteristic geometrical properties of solutions that are responsible for the usage of specific algorithms and constitute the possible bottleneck for the near future.

Interestingly enough, the changes in both time or memory requirements during the solution process of sparse systems can be interpreted in physical terms as a dynamical transition at which the phase space of the associated physical systems becomes split into an exponential number of ergodic components. While it is to be expected that local algorithms get stuck by local minima at such phase boundary, it is less obvious to predict which is the counterpart of the dynamical transition in global algorithms, for which polynomial time convergence is guaranteed even for the hardest instances. Indeed the dynamical transition manifests itself as a phase transition in the computational requirements which in turn leads to a slowing down phenomenon that saturates the upper bound for the convergence time. Such a change of scale in memory requirements constitute a serious problem for hardware implementations of large scale simulations.

Hopefully, the paper is written in a style accessible to a cross disciplinary audience.

## II Random Linear systems in GF(2): rigorous results and statistical mechanics analysis

The theory of random equations in finite fields is shared by probability, combinatorics and algebra KOLCHIN .

For the sake of simplicity we limit our statistical mechanics analysis to GF(2) rather than GF($`q`$), the extension to $`q>2`$ being straightforward though technically involved.

As is well known in the context of error correcting codes Sourlas , solving a sparse linear system modulo 2 is equivalent to finding the zero temperature ground states of a class of multiple degree interactions $`p`$-spin models on diluted random graphs.

Let us consider a random linear system in GF(2) in the form $`\widehat{A}\stackrel{}{x}=\stackrel{}{y}\text{mod}\mathrm{\hspace{0.33em}2}`$, where $`\widehat{A}`$ is a 0-1 matrix of dimension $`M\times N`$. For each of its specific choices $`\widehat{A}`$ can be interpreted as the contact matrix of a particular random (hyper-)graph belonging to a specific ensemble. The class of random matrices we shall deal with are defined by the fraction of rows $`a_k`$ with $`k`$ non zero elements. The latter are placed uniformly at random within each row.

We focus on matrices that lead to graphs with an average connectivity value $`k=_ka_kk`$ finite and much less than both $`M`$ and $`N`$. We are interested to the limit of very large matrices, where we can assume $`N,M\mathrm{}`$ with a finite ratio $`\gamma M/N`$.

This is the regime in which a study of the computational cost is important in that it applies directly to large scale computations. In the limit $`N,M\mathrm{}`$ average quantities characterizing the system (e.g. the average fraction of violated equations) are known to be equal to the most probable values (i.e. their probability distribution is strongly peaked self-average ) and therefore single random large systems behave as the average over the ensemble.

We will always assume $`a_1=0`$ at the beginning, since rows with a single one corresponds to trivial equations which can be removed a priori from the set.

The equivalence between linear systems and spin models is quite straightforward. We start from a set of linear equations in GF(2), $`\widehat{A}\stackrel{}{x}=\stackrel{}{y}`$, and we build up a spin Hamiltonian whose ground state energy $`E_{gs}`$ counts the minimal number of unsatisfied equations. In the case where $`E_{gs}=0`$, ground state configurations will correspond to solutions of the original set of linear equations and the zero-temperature entropy will count the number of such solutions.

The construction is done as follows: For every equation, labelled by $`i[1\mathrm{}M]`$, let us define the set of variables $`\stackrel{}{x}`$ entering equation $`i`$ as

$$v(i)\{j[1\mathrm{}N]:A_{ij}=1\}.$$
(1)
With the transformation $`s_j=(1)^{x_j}`$ and $`J_i=(1)^{y_i}`$, we have that every equation can be converted in a term of the Hamiltonian through

$$\underset{j=1}{\overset{N}{}}A_{ij}x_j=y_i\underset{jv(i)}{}x_j=y_i\underset{jv(i)}{}s_j=J_i,$$
(2)
where the multi-spin interaction contain at least 2 spins since we set $`a_1=0`$. Then the Hamiltonian

$$H=\frac{1}{2}\left[M\underset{i=1}{\overset{M}{}}J_i\underset{jv(i)}{}s_j\right],$$
(3)
fits the above requirements and can be used in the analytical treatment.

A better form for the above Hamiltonian can be obtained grouping together $`k`$-spin terms with the same $`k`$, that is

$$H=\frac{1}{2}\left[M\underset{k}{}\underset{i_1<i_2<\mathrm{}<i_k}{}J_{i_1i_2\mathrm{}i_k}s_{i_1}\mathrm{}s_{i_k}\right],$$
(4)
where $`s_i=\pm 1`$ are Ising spins and the couplings $`J_{i_1i_2\mathrm{}i_k}`$ are i.i.d. quenched random variables taking values in $`\{0,\pm 1\}`$. The total number of interactions, that is of terms with $`J0`$, is M, and the energy is zero if and only if all the interactions are satisfied. For each unsatisfied interaction the energy increases by 1.

The fraction of interactions of $`k`$-spin kind is $`a_k`$ and thus the probability of having $`J_{i_1i_2\mathrm{}i_k}0`$ equals $`a_kM/\left(\genfrac{}{}{0pt}{}{N}{k}\right)\gamma a_kk!/N^{k1}`$, while the sign of $`J_{i_1i_2\mathrm{}i_k}`$ depends on the probability distribution of the components of $`\stackrel{}{y}`$,

$$P(J_{i_1i_2\mathrm{}i_k})=\left[1\frac{\gamma a_kk!}{N^{k1}}\right]\delta (J_{i_1i_2\mathrm{}i_k})+\frac{\gamma a_kk!}{N^{k1}}\left[p\delta (J_{i_1i_2\mathrm{}i_k}1)+(1p)\delta (J_{i_1i_2\mathrm{}i_k}+1)\right],$$
(5)
where $`p[0,1]`$ controls the fraction of zeros in $`\stackrel{}{y}`$. As long as the system admits at least one solution, it can always be brought by a gauge transformation in the form with $`p=1\stackrel{}{y}=\stackrel{}{0}`$. This corresponds to have positive or null couplings only, like in a diluted ferromagnetic model.

In order to make a connection between the behavior of solving algorithms and the structure of the matrix $`\widehat{A}`$, we study the geometrical properties of the space of solution, i.e. ground states of (4), as a function of $`\gamma `$ for non-trivial choices of $`\{a_k\}`$. We may have access to the structure of such a space by just performing the $`T=0`$ statistical mechanics analysis of the spin glass model, with control parameter $`\gamma `$.

For $`\gamma `$ large enough, at say $`\gamma _c`$, the system of equations becomes over-determined and some of the equations can no longer be satisfied. This fact is reflected in the ground state energy of the associated spin glass model becoming positive. The interesting aspect of the problem is that, under proper conditions, there appears a clustering phenomenon with macroscopic algorithmic consequences at some intermediate value $`0<\gamma =\gamma _d<\gamma _c`$. We will focus our attention on the latter transition, thus assuming a priori that at least one solution always exist. This allow us to fix $`\stackrel{}{y}\stackrel{}{0}`$ hereafter.

The complete picture of the typical structure of the solution space can be obtained through a replica calculation, following a well tested scheme in diluted systems, as in RIWEZE  and LERIZE . The results of such a calculations have been recently confirmed by a rigorous mathematical derivation MERIZE ; CDMM . We avoid here to repeat standard calculations already presented in RIWEZE ; LERIZE  and extensively reviewed in LEONE , and we directly present the results.

Due to the zero energy condition ($`E_{gs}=0`$ for $`\gamma <\gamma _c`$), the dominance of thermodynamical states is purely to be determined in entropic terms. Defining $`S_0(\gamma )`$ as the logarithm of the number of solutions to $`\widehat{A}\stackrel{}{x}=\stackrel{}{0}`$ divided by $`N`$, we have that

$$S_0(\gamma )=S(m,\gamma )=\mathrm{log}(2)\left[(1m)[1\mathrm{log}(1m)]\gamma \underset{k2}{}a_k(1m^k)\right],$$
(6)
where $`m`$ solves

$$G(m)=1me^{\gamma _{k2}ka_km^{k1}}=0.$$
(7)
When more than one solution to Eq.(7) exist, the one maximazing $`S(m,\gamma )`$ must be chosen.

At fixed $`\{a_k\}`$, one can study the phase diagram as a function of $`\gamma `$. At low enough $`\gamma `$, Eq. (7) has only the trivial solution $`m=0`$ and the system is paramagnetic with entropy $`S(0,\gamma )=\mathrm{log}(2)(1\gamma )`$. Typically a non trivial magnetized solution for the order parameter, $`m^{}>0`$, appears at a value $`\gamma _d`$ such that

$$G(m^{})=0\text{ and }\frac{G(m)}{m}|_{m=m^{}}=0.$$
(8)
This solution becomes entropically favored at a value $`\gamma _c`$ found solving

$$S(0,\gamma _c)=S(m^{},\gamma _c).$$
(9)
The crucial observation is the following. At $`\gamma _d`$, together with the magnetized solution, there appear other spin glass solutions to the saddle-point equation. In particular, it can be shown FLRZ ; MERIZE  that the difference between the paramagnetic and the ferromagnetic entropies,

$$\mathrm{\Sigma }(\gamma )=S(0,\gamma )S(m^{},\gamma ),$$
(10)
gives the configurational entropy of the problem, that is the number of clusters of solutions 111Two solutions belong to the same cluster (resp. to different clusters) if their Hamming distance is $`𝒪(1)`$ \[resp. $`𝒪(N)`$\].. There exist $`\mathrm{exp}[\mathrm{\Sigma }(\gamma )N]`$ well separated clusters \[Hamming distances $`𝒪(N)`$\], each one containing a number $`\mathrm{exp}[S(m^{},\gamma )N]`$ of closed solutions \[Hamming distances $`𝒪(1)`$\].

This clusterizations has two main consequences. Local algorithms for finding solutions running in linear time in $`N`$ stop converging RIWEZE : this is the typical situation for greedy algorithm which get stuck in one of the most numerous local minima at a positive energy.

Global algorithms, which are guaranteed to converge in polynomial time, need to keep track along computation of this complex structure of solutions and a memory linear in $`N`$ turns out to be insufficient, as we will show below.

The simple cases $`\{a_2=1;a_{k2}=0\}`$ and $`\{a_3=1;a_{k3}=0\}`$ are particularly illuminating LERIZE . Self consistency equations for the order parameter and for the ground state entropy can be immediately retrieved from the general formulas (7) and (6).

The first case represents a simple Ising spin model on a random graph. The analysis of thermodynamic phases shows a trivial paramagnetic region at low $`\gamma `$, followed by a second order phase transition at $`\gamma _d=\gamma _c=1/2`$ where both the ground state energy and the magnetization become positive in a continuous way. From the linear system point of view, working with sparse equations with only 2 variables per row is always easy, i.e. algorithms have no slowing down. The solving process progressively fixes variables and smoothly goes on until a complete solution has been found. Indeed, fixing one variable immediately fixes the other one within the equation, and this goes on in a cascade process that prevents the accumulation of too long symbolic memory-taking expressions.

The second case of the 3-spin model was tackled in full detail in RIWEZE ; FLRZ . It shares the general characteristics of all models without, or at most with a very small fraction of, 2-spin terms. In this case a gap between $`\gamma _d`$ and $`\gamma _c`$ opens up and a non vanishing configurational entropy $`\mathrm{\Sigma }`$ appears there. At $`\gamma _c`$ the magnetization jumps to a finite value.

For a general choice of $`\{a_k\}`$, the configurational entropy reads

$$\mathrm{\Sigma }(\gamma )=\mathrm{log}(2)\left[1(1m)[1\mathrm{log}(1m)]+\gamma \underset{k2}{}a_km^k\right],$$
(11)
where $`m`$ is the largest solution to Eq. (7). As discussed in Ref. FLRZ , the values of $`\gamma _d`$ and $`\gamma _c`$ are found as the points where $`\mathrm{\Sigma }(\gamma )`$ first appears with a non zero value and where it reaches zero again.

The algorithmic consequences of having $`\mathrm{\Sigma }(\gamma )>0`$ have been already exposed in Refs. RIWEZE ; FLRZ : For $`\gamma >\gamma _d`$ a glassy state with positive energy arises, which traps any local dynamics, preventing it to converge towards the ground state of zero energy. We conjecture the counterpart on global algorithms, such as Gaussian elimination, to be that the resolution time increases with $`N`$ faster than linear.

In the next section we will check the above conjecture with two different Gaussian elimination algorithms, none of which is able to solve the system in linear time for $`\gamma >\gamma _d`$.

## III Algorithms behaviour

In this section we analyze the performances of a couple of different ‘Gaussian elimination’ algorithms, their difference being in the order equations are solved. We will measure the number of operations and the size of the memory required for the solution of a set of linear equations, that is the complexity for finding all solutions to $`\widehat{A}\stackrel{}{x}=\stackrel{}{y}`$.

We will see that, for a generic ensemble of random problems, any algorithm undergoes an easy/hard transition at a certain $`\gamma `$ value, which can not be pushed beyond the dynamical transition threshold $`\gamma _d`$. In this context we call easy such problems which are solvable with a CPU-time and memory of order $`N`$, and hard those requiring resources scaling with $`N^\alpha `$, where $`\alpha >1`$.

Given a set of $`M`$ linear equations in $`N`$ variables, Gaussian elimination proceeds as follows \[for concreteness we will always work in GF(2)\]: At each step, it takes an equation, e.g. $`x_1+x_2+x_3=y_1`$, solves it with respect to a variable, e.g. $`x_1=x_2+x_3+y_1`$, and then it substitutes variable $`x_1`$ with the expression $`x_2+x_3+y_1`$ in all the equations still unsolved. This procedure gives all the solutions to any set of linear equations in, at most, $`𝒪(N^3)`$ steps and using $`𝒪(N^2)`$ memory. Nevertheless this bounds only holds in the worst case, namely when the matrix $`\widehat{A}`$ is dense. Very often, in actual applications, the matrix is sparse and the algorithm is faster. We define sparse a matrix with $`𝒪(N)`$ ones and dense that with $`𝒪(N^2)`$ ones.

In order to analyze the computational complexity of this problem, and its connections to phase transitions, we focus on a specific ensemble of random problems, generalizations to other ensembles being straightforward. We choose sets of $`M=\gamma N`$ linear equations, each one containing exactly $`k=3`$ of the $`N`$ variables, taking values in GF(2). Thus the connectivity of a variable, defined as the number of equations this variable enters in, takes values from a Poissonian distribution of mean $`3\gamma `$.

For very large $`N`$, that is in the thermodynamical limit, we are interested in how the complexity changes with $`\gamma `$. Moreover, for a fixed $`\gamma `$ such that the problem is hard, we would like to know when (in terms of the running-time $`t`$) and why the algorithm becomes slower and slower.

The running-time $`t`$ is measured as the number of equations already solved, normalized by $`N`$, and thus takes values in $`[0,\gamma ]`$. $`\widehat{A}_t`$ is the matrix representing the set of equations after $`tN`$ steps, and it has the form shown in Fig. 1. See Fig. 2 for the actual shape of $`\widehat{A}_t`$ in a specific case with 1024 equations in 1024 variables. For ease of simplicity, we have reordered the variables and the equations of the system, such that, at the $`i`$-th step, we solve the $`i`$-th equation with respect to $`x_i`$. With this choice the left part of the matrix $`\widehat{A}_t`$ has ones on the diagonal and zeros below. The right part can be naturally divided in an upper part $`U`$ and a lower one $`L`$. The density of ones in the $`L`$ matrix — let us call it $`\rho (t;\gamma )`$ — is uniform and depends on the initial $`\gamma `$, the time $`t`$ and the algorithm used for solving the linear system. The density of ones in the $`U`$ part is not uniform and varies from row to row, as shown in Fig. 1 with gray tones. For continuity reasons the density at the $`m`$-th row of $`U`$ is exactly $`\rho (m/N;\gamma )`$. Then $`U`$ is sparse or dense depending on whether $`L`$ is. Defining $`k(t;\gamma )=\rho (t;\gamma )N(1t)`$ the average number of ones per row in $`L`$, we have that a sparse (resp. dense) matrix corresponds to having a finite $`k`$ (resp. $`\rho `$).

At each time step, the number of operations required are directly related to the density of the matrix $`\widehat{A}_t`$ and thus to that of $`L`$. More specifically, solving with respect to the variable in the upper left corner of $`L`$, the number of operations is proportional to the number of ones in the first row of $`L`$, i.e. $`k(t;\gamma )`$, times the number of rows of $`L`$ having a one in the first column, i.e. $`\rho (t;\gamma )N(\gamma t)`$, and thus equals

$$k(t;\gamma )\rho (t;\gamma )N(\gamma t)=k^2\frac{\gamma t}{1t}=N^2\rho ^2(\gamma t)(1t).$$
(12)
Then, if the matrix $`L`$ is sparse a finite number of operations per step is enough, while $`𝒪(N^2)`$ operations are required when $`L`$ is dense. Integrating over time $`t[0,\gamma ]`$, we have that the total complexity is given by

$$N_0^\gamma \frac{\gamma t}{1t}k^2(t;\gamma )𝑑t=N^3_0^\gamma (\gamma t)(1t)\rho ^2(t;\gamma )𝑑t.$$
(13)
Since the function $`\rho (t;\gamma )`$ is continuous in $`t`$, we conclude that

$`\begin{array}{c}\rho (t;\gamma )1/N\\ k(t;\gamma )\text{ finite}\end{array}\text{, for all }t[0,\gamma ]\}`$ $`\rho _{max}(\gamma )=0`$ $`\{\begin{array}{c}\text{CPU time}N\\ \text{Memory}N\end{array}`$ (18)
$`\begin{array}{c}\rho (t;\gamma )\text{ finite}\\ k(t;\gamma )N\end{array}\text{, for some }t[0,\gamma ]\}`$ $`\rho _{max}(\gamma )>0`$ $`\{\begin{array}{c}\text{CPU time}N^3\\ \text{Memory}N^2\end{array}`$ (23)
where

$$\rho _{max}(\gamma )=\underset{N\mathrm{}}{lim}\underset{t[0,\gamma ]}{\mathrm{max}}\rho (t;\gamma )$$
(24)
is the order parameter signaling the onset of the hard regime.

Having found the relation between the density of ones in $`L`$ and the computational complexity we are interested in, we can now run the algorithms and measure the density $`\rho (t;\gamma )`$. The easy/hard transition should manifest itself with $`\rho _{max}(\gamma )`$ becoming different from zero.

### III.1 Simplest Gaussian elimination

Let us start with the simplest algorithm, which solves the equations in the same (random) order they appear in the set and with respect to a randomly chosen variable. In this very simple case, one can easily show that the complexity for solving a set of linear equations with initial parameter $`\gamma =\gamma _0`$ is exactly the same as for solving a larger system with $`\gamma >\gamma _0`$ up to time $`t=\gamma _0`$. For this reason, in this case the function $`\rho (t;\gamma )`$ does not depend on $`\gamma `$ and can be calculated once for all the relevant $`\gamma `$ values.

Moreover, it is known RIWEZE  that this algorithm, in the limit of very large $`N`$, keeps the matrix sparse for all $`\gamma <2/3`$.

In Fig. 3 we show the function $`\rho (t)`$ for many large $`N`$ values. The dotted-dashed line is a guide to the eyes and it should not be too much different from the thermodynamical limit: It goes through the two points ($`\gamma =2/3`$ and $`\gamma =0.918`$) where $`\rho (t)`$ must vanish and coincide with numerical data in the region, where data for different sizes seem to be quite close to the asymptotic shape.

In the thermodynamical limit, the algorithm keeps the matrix sparse for times $`t2/3`$ and so it undergoes an easy/hard transition at $`\gamma =2/3`$: As long as $`\gamma 2/3`$, $`\rho _{max}(\gamma )=0`$, while $`\rho _{max}(\gamma )>0`$ for $`\gamma >2/3`$. As we will see below the location of the transition depends on the algorithm used and, in this case, does not correspond to any underlying thermodynamical transition.

We note en passant that the $`\gamma `$ value where the $`L`$ matrix becomes sparse again seems to correspond to the critical point $`\gamma _c=0.918`$ RIWEZE ; MERIZE ; CDMM  (marked with a vertical line in Fig. 3). An explanation to this observation will be given in a forthcoming publication. It implies that the value of the critical point $`\gamma _c`$, which is relevant e.g. in the XORSAT model XOR-SAT  in theoretical computer science, could be obtained also by solving differential equations for $`\rho (t)`$.

### III.2 Smart Gaussian elimination

Now we turn to a more clever Gaussian elimination algorithm, which works as follows: At each time step, it chooses the variable $`x`$ having the smallest connectivity in $`L`$, i.e. that corresponding to the less dense column of $`L`$, and solves with respect to $`x`$ any of the equations where $`x`$ enters in.

Clearly, in this case, the dynamics and thus the density of ones in $`L`$ depend on the initial $`\gamma `$ value: A smaller $`\gamma `$ implies that for a longer time we can choose variables of connectivity 1, which do not increase the average number of ones per row in $`L`$. It can be rigorously shown MERIZE ; CDMM  that this procedure keeps the density of the $`L`$ matrix constant, $`\rho (t;\gamma )=\rho (0;\gamma )`$, for times smaller than $`t^{}=\gamma (1m^3)`$, where $`m`$ is the largest solution to $`1m=\mathrm{exp}(3\gamma m^2)`$. The last equation is exactly Eq. (7) with $`\{a_3=1,a_{k3}=0\}`$.

Running the algorithm for different $`\gamma `$ values we obtain the densities reported in the main panel of Fig. 4. For $`\gamma <\gamma _d=0.818`$ the density remains $`𝒪(1/N)`$ all along the run, while for $`\gamma >\gamma _d`$ there is a time when the density becomes finite and the problem hard to handle.

In order to better show what happens around $`t^{}`$, we have plotted in the inset of Fig. 4 the mean number of ones per row, $`k(t)`$. It is clear that for $`\gamma <\gamma _d`$ this number remains constant, since one can solve the system choosing only variables of connectivity 1, not altering the $`L`$ matrix. On the contrary, for $`\gamma \gamma _d`$ there is a time $`t^{}(\gamma )`$ when variables of connectivity 1 terminate, and the algorithm has to start making substitutions in $`L`$, thus increasing the density of ones.

Then $`\gamma _d`$ marks the onset of computational hardness, both in memory and CPU time. One may object that also this value for the easy/hard transition may depend on the particular algorithm. Note, however, that a completely different linear algorithm described in Ref. RIWEZE  (which firstly works with high-connectivity variables) seems to work up to $`\gamma _d`$. Moreover, as seen in the previous section, we have analytically found that at $`\gamma _d`$ a transition takes place, which drastically changes the structure of the solutions space, and so we argue that any algorithm running in linear time can converge only up to $`\gamma _d`$. Indeed is shown in MERIZE  that solutions spontaneously form clusters for $`\gamma >\gamma _d`$ and this particular structure requires a larger memory to be stored.

## IV The RSA cryptosystem and factorization

In this section we shall validate the above scenario on a concrete application, namely integer factorization problems arising in the RSA cryptosystem. Such problems allow for a non-trivial mapping onto huge linear systems in GF(2) with a rather peculiar structure of the underlying contact matrix. In order to be as self-contained as possible, we firstly give a short review of the problem and the methodology (a detailed description of the RSA cryptosystem can be found in RSA ).

The only known method for breaking RSA implies factorization of the private key, which consists in a natural number which is the product of two big prime numbers, $`n=pq`$, with $`p`$ and $`q`$ approximately of the same size $`\sqrt{n}`$. Keys currently used in applications are numbers $`n`$ ranging from 1024 bits (309 decimal digits) to 2048 bits (617 digits) length.

The first attempt at a massive parallel factorization was the RSA129 (129 digits, 428 bits) challenge, solved in 1994 with the *quadratic sieve* (QS) algorithm. More recently, in August, 1999 the RSA155 challenged was solved using the *general number field sieve* (*GNFS*) algorithm. This has forced to abandon the 512-bit (155 digits) length for sensitive information security.

There are now several sub-exponential algorithms for solving the factorization problem, the faster of which is *GNFS. QS* and *GNFS* share the same structure, consisting of two phases: a first one in which a big (the size depending mostly on the size of $`n`$) linear system in GF(2) is produced, and a second one in which this system is solved.

Although the first phase is definitely more costly, the solving phase (which affect this paper) takes a respectable part of the total time and memory requirement. Especially as numbers get bigger this becomes a limitation, because the fastest solving methods used employ a sole workstation, with the consequent memory restriction. Moreover, in recent factorizations a new filtering phase has been placed between the previous two, in which pieces of the system (specifically columns of the $`\{0,1\}`$matrix) get discarded in order to simplify the solving phase, effectively transferring part of the total time from the second phase to the first one.

In this section we will study statistical characteristics of linear systems produced by the QS algorithm. Next subsection is dedicated to a schematic description of QS.

### IV.1 The QS algorithm

For a nice description of the QS algorithm see CRYPTO . Synthetically, QS works at follows. It builds a list of integer numbers $`\{y_i\}_{iI}`$ such that:

* $`y_ix_i^2(\text{mod}n)`$ for some $`x_i`$ and $`y_ix_i`$;
* $`y_i`$ is completely factorizable in a given (relatively small) subset of $`B`$ primes called the factor-base.

This is called the *sieving* phase. The algorithm then searches a subset $`JI`$ of elements of the list such that $`_{iJ}y_i=z^2`$ is a square (*solving* phase). Once found, $`z^2x^2(\text{mod}n)`$ (here $`x=_{iJ}x_i`$) and this implies that $`n`$ divides $`(x+z)(xz)`$ and then $`gcd(xz,n)`$ will likely (further trials will increase the probability) be a non-trivial factor of $`n`$.

#### IV.1.1 Sieving

In order to find element pairs $`x_i`$,$`y_i`$ such that $`y_ix_i^2(\text{mod}n)`$ we can use the polynomial $`y=f(x)=x^2n`$ and evaluate it at different values of $`x`$, keeping only values of $`y`$ which completely factorize between the first $`B`$ primes (the factor-base). The sieving will allow us to do this efficiently.

The idea is that, given $`p`$, it is easy to find which are the values of $`f(x)`$ which are divisible by $`p`$, because $`p`$ divides $`f(x)`$ if and only if $`f(x)=x^2n0(\text{mod}p)`$ and this is a quadratic equation in GF(p), having at most 2 solutions. These solutions are nothing but the square roots of $`n`$ modulo $`p`$ (if they exist).

This has a first consequence, i.e. that a prime $`p`$ will not divide $`f(x)`$ if $`n`$ is not a square $`\text{mod}p`$ independently of the value of $`x`$. So if we can detect these primes, we can eliminate them directly from our set of primes. Detecting them is very easy: Using Fermat’s little theorem, we know that

$$n^{p1}1(\text{mod}p),$$
(25)
assuming that $`p`$ do not divide $`n`$ (which is trivially a reasonable assumption, anyway, because we are searching a divisor of $`n`$). If $`p`$ is an odd prime, i.e. not 2 (all $`n`$ are a squares $`(\text{mod}\mathrm{\hspace{0.17em}2})`$), then calling $`m=n^{\frac{p1}{2}}`$ we have that $`m^21(\text{mod}p)`$, so $`m\pm 1(\text{mod}p)`$. This $`m`$ will prove to be handy.

If $`ns^2(\text{mod}p)`$ then $`ms^{p1}1(\text{mod}p)`$. Conversely, if $`m1`$ then $`n`$ is a square modulo $`p`$ (not proven here). The number $`m`$ is called the Lagrange symbol and can be computed efficiently in one of the firsts stages of the algorithm. Useful primes (those with $`m=1`$) are roughly a random half of all the first $`B`$ primes. So now we will keep only this half and redefine “the first $`B`$ primes” as “the first $`B`$ primes with $`m=1`$”.

Computing the square root modulo $`p`$ is a bit more difficult than knowing that it exists, but can also be done efficiently. For instance, the easiest case is when $`\frac{p+1}{4}`$ is an integer, then $`\left(n^{\frac{p+1}{4}}\right)^2=n^{\frac{p+1}{2}}mn(\text{mod}p)`$. As $`m=1`$ (or else there is no solution) then $`\pm n^{\frac{p+1}{4}}(\text{mod}p)`$ are the required square roots.

Once we have computed the two solutions $`f(x_p^{1,2})0(\text{mod}p)`$, then adding $`p,2p,3p,\mathrm{}`$ to them we will obtain *all* $`x`$ such that $`f(x)`$ is divisible by $`p`$.

The sieve idea is to initialize an array with values of $`f(x)`$ for consecutive $`x[[\sqrt{n}],[\sqrt{n}]+M]`$ indexed by $`x`$, and then for each $`p`$ in our factor base to divide the corresponding arithmetic progression of $`\{f(x_p^{1,2}+kp),k=1,\mathrm{}\}`$ by $`p`$. At the end those values which are completely factored between the primes in the factor base will become $`1`$ (Well, not exactly. Some of them can have multiple times the same prime factor. But we can set up a threshold instead of $`1`$ below which we consider the number completely factored. We can recheck afterwards). We take those values and put their factorization in an array

$$\begin{array}{ccc}& \begin{array}{ccc}f(x_1)& \mathrm{}& f(x_m)\end{array}& \\ \begin{array}{c}p_1\\ \mathrm{}\\ p_B\end{array}& \left[\begin{array}{ccc}\alpha _1^{(1)}& \mathrm{}& \alpha _1^{(m)}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \alpha _B^{(1)}& \mathrm{}& \alpha _B^{(m)}\end{array}\right]& (\text{mod}\mathrm{\hspace{0.17em}2})\end{array}$$

#### IV.1.2 Solving

The *solving* phase is conceptually simple: A solution of the homogeneous linear system $`\widehat{A}v=0`$ is a $`\{0,1\}`$vector $`v`$ which represent correctly the subset $`J`$, in the sense that $`v_i=1`$if and only if $`iJ`$.

### IV.2 The matrix ensemble

#### IV.2.1 Correlations

We have implemented the simplest QS described in CRYPTO  in order to analyze the output matrix ensemble. We attempted to look for correlations in the presence/absence of different primes in the set of divisors of the variables $`y_i`$. Specifically we checked that there is virtually no correlation between rows of the matrix: We have taken one such output matrix (resulting from the factorization of a product of two $`20`$ digits primes) and computed the covariance between the corresponding spin variables $`s_1,s_2`$ of two rows $`r_1,r_2`$, the averages being taken along different columns,

$$s_1s_2s_1s_2.$$

Once repeated for all $`r_1<r_2`$, we found that all pairs have correlations in the interval $`0\pm 0.06`$, a proportion of $`0.9999`$ pairs having correlations in $`0\pm 0.02`$.

#### IV.2.2 Dependence on “factorization hardness”

We then examined dependence of the resulting distributions of ones per row on the “factorization hardness” of the number $`n`$. Typically (depending on the algorithm) the complexity of factorization depends on the size of the smallest prime divisor of $`n`$ 222This is why in RSA we choose $`n=pq`$ with $`p,q\sqrt{n}`$.: For instance, *trial division* ends in exactly this amount of steps. It was conjectured that this would be reflected in the structure of the output matrix.

We have constructed 25 numbers $`n`$ with different factor sizes (from now on, factor type 10+10+10 will mean a 30 digit number constructed as a product of tree 10-digit primes) organized as follows:

* 5 of type 20+20
* 5 of type 10+30
* 5 of type 13+13+13
* 5 of type 10+10+10+10
* 5 of type 5+5+5+5+5+5+5+5

All 25 numbers differed between them in less than a 0.01%. We then made QS compute the factorization matrices, with a factor base of size 1500. This value for the size of the factor base has been chosen experimentally in order to minimize the sieving phase duration. The resulting matrices were of size $`1500\times 1510`$ and were then post-processed in order to remove rows and columns with a single 1. The final size is thus reduced of about 200 columns and rows. The resulting distributions of ones per row are plotted in Fig. 5, showing very little variations. They can be very well described by a unique distribution, which is substantially a power law with some little deviations in the range of type-2 and type-3 rows. The best fit in the region $`X>3`$ gives an exponent $`2.2`$.

Our conclusion is that statistical properties of the resulting matrix do not depend on the factorization hardness. The bottleneck for factorizing a large hard number is mainly determined by the time required by QS to generate the matrix, which indeed strongly depends on the size of the smallest factor. In the rest of the paper we will analyze the solving phase, assuming the factorization matrix to have uncorrelated rows and the number of ones per row to be a random variable extracted from distribution in Fig. 5. These two assumptions have been experimentally verified.

### IV.3 Linear solving methods

Plain standard Gaussian elimination execution time is cubic in the size of the matrix (our matrices are almost square). Fortunately, we can pack $`32`$ matrix entries in a single $`4`$-byte word, and then the sum operation is implemented as the low-cost bit-wise logical XOR operation, saving a factor 32 in time.

As also the matrix is very sparse, instead of keeping in memory all of it, we can memorize only the position of 1’s. This forbid us to use the factor-32 trick, but allows us to do the first steps very quickly. At some time in the Gaussian elimination process (typically more than half of the process), the remaining (non eliminated) part will be very dense, and then it will be convenient to switch to the standard method above. This is what was done in the solving phase of RSA129.

Another option is to use in one of the stages an iterative algorithms, like the discrete Lanczos. The Lanczos method has the advantage of having a stable $`𝒪(N^2)`$ total time for a sparse matrix, but finds only one solution (or a prefixed quantity in the block-Lanczos variant) instead of all of them. For factorization this is not a problem, because we need only a few solutions to have a reasonable chance. This is the method that was used in the solving phase of RSA155.

### IV.4 Power law distributed $`\{a_k\}`$: Phase diagram and comparison with real application data

The previous analysis leads to the construction of matrices whose density of non zero entries follows quite well a power law distribution with light deviations due to rows with a small number of ones and a cutoff, $`k_{max}`$, of some hundreds. Then we use the following distribution in the analytical treatment:

$`a_2`$ $`=`$ $`a,`$ (26)
$`a_k`$ $`=`$ $`ϵk^s\text{for }3kk_{max},`$ (27)
where $`ϵ`$ is a normalizing factor equal to $`(1a)/_{k=3}^{k_{max}}k^s`$. The factorized integers considered in the previous section lead to an exponent $`s2.2`$ and to a non zero support up to $`k_{max}200`$. The choice of keeping $`a_2`$, and only $`a_2`$, as an independent parameter is dictated by the very difference in the physical behavior of 2-spin terms and $`k`$-spin terms with $`k>2`$.

The study of the phase diagram in the control parameter $`\gamma `$ for choices of $`a`$, $`s`$ and $`k_{max}`$ retrieved from real data reveals a non trivial behavior.

In Fig. 6 we show the phase diagram for $`s=2.2`$ and $`k_{max}=200`$. Only part of the entire phase diagram ($`a[0,1]`$, $`\gamma [0,1]`$) is shown for clarity. The lines further go on smoothly outside the drawn portion.

If $`a`$ is high enough, we are in the rightmost region I of the phase diagram, where algorithms smoothly find solutions to the system and do not undergo any critical slowing down. Indeed, crossing the bold iperbole $`\gamma =1/(2a)`$ given by the condition $`G(m)/m|_{m=0}=0`$, the system undergoes a continuous transition in the order parameter $`m`$, representing the fraction of variables taking the same value in all the solutions. The problem of finding solutions is always easy, as for the case $`\{a_2=1;a_{k2}=0\}`$ explained before.

Decreasing $`a`$ we meet a first intermediate region II, where the birth of a meta-stable non-trivial saddle-point solution at $`\gamma =\gamma _d(a)`$ is given by the solution of Eq. (8). However, algorithms should not be much affected by this meta-stable state, because the system starts magnetizing continuously before, crossing the bold line. Increasing $`\gamma `$ up to the critical value $`\gamma _c(a)`$ one meets a first order transition, where the magnetization, that was already non-zero, undergoes a further jump.

The second central region III shows an inversion between $`\gamma _d(a)`$ and the bold line $`1/(2a)`$. These two intermediate regions have not been exhaustively studied yet, because real data all fall in the leftmost one. The shape of the central part of the phase diagram is very sensitive to the choice of the control parameters $`s`$ and $`k_{max}`$, as shown in Fig. 7.

The $`\gamma _c(a)`$ curve in the second and third regions is found solving

$$S(m^{},\gamma _c)=S(m_{},\gamma _c),$$
(28)
where $`m_{}`$ is the smallest positive solution to $`G(m)=0`$, which corresponds to the magnetization of the ferromagnetic state arisen from the second order transition (bold line). The points of crossing showing the onset of different regions, from right to left, are found respectively as: $`\frac{G(m)}{m}=0\&S(m^{},\gamma )=S(m_{},\gamma )`$, $`\frac{G(m)}{m}=0\&\gamma =\frac{1}{2a}`$ and $`S(m^{},\gamma )=S(0,\gamma )\&\gamma =\frac{1}{2a}`$.

The leftmost part IV shows the typical behavior described in RIWEZE . Increasing $`\gamma `$ the system never reaches the continuous transition on the bold line, but it undergoes a first dynamical transition at $`\gamma _d(a)`$ and second thermodynamical one at $`\gamma _c(a)`$, found via Eq. (28) with $`m_{}=0`$ since we are still below the second order transition line. Configurational entropy is non-zero between $`\gamma _d(a)`$ and $`\gamma _c(a)`$, and solving algorithms are affected by it. There are typically other spinodal lines in the phase diagram, but they always correspond to sub-optimal solutions, and were, therefore, not shown in the picture.

In real data the fraction of 2-variables equations is typically of the order of 0.2 and $`\gamma 1`$. So we always work deep into phase IV where, during the solving procedure, the system undergoes a first dynamical transition, that corresponds to a slowing down of the solving algorithms, before finding solutions.

## V Conclusions

We have analyzed the behaviour of different type of polynomial algorithms in the solutions of large-scale linear systems over finite fields. The connection between memory requirements and clustering phase transitions as been made clear on both artificially generated problem as well as on a “real-world” applications. While the role of the dynamical glass transition in local search algorithm was already well known (trapping in local minima), we have provided a clear example of the role of such type of glass transition in global dynamical processes which are guaranteed to converge to the global optimum in some polynomial time. The memory catastrophe found is such cases constitutes a concrete limitation for the performance of single-machine programs.


text:
# Phase transition and critical behaviour of the d=3 Gross-Neveu model

## Abstract

A second order phase transition for the three dimensional Gross-Neveu model is established for one fermion species $`N=1`$. This transition breaks a parity-like discrete symmetry. It constitutes its peculiar universality class with critical exponent $`\nu =0.63`$ and scalar and fermionic anomalous dimension $`\eta _\sigma =0.31`$ and $`\eta _\psi =0.11`$, respectively. We also compute critical exponents for other $`N`$. Our results are based on exact renormalization group equations.

An understanding of systems with many fermionic degrees of freedom is one of the big challenges in statistical physics. Due to the anticommuting nature of the variables numerical simulations are not straightforward—analytical methods are crucially needed. One typically has to solve a functional integral for a $`d`$ dimensional system with Grassmann variables. Approximate solutions for “test models” would be of great value. The Gross-Neveu (GN) model  is one of the simplest fermionic models. In three dimensions a discrete symmetry forbids a mass term unless it is spontaneously broken. In the symmetric phase the GN model is therefore a realization of a statistical system of gapless fermions. For a large number $`N`$ of fermion species it is known  that a second order phase transition separates the symmetric phase from an ordered phase where the symmetry is spontaneously broken and the fermions become massive. Using methods based on an exact renormalization group equation  a second order transition for $`N2`$ was confirmed. We know, however, of no previous work which clarifies the existence and nature of the phase transition in the simplest model with only one fermion species. The model with one fermion species is inaccessible to lattice simulations due to the fermion doubling problem and the $`1/N`$ expansion is not expected to give reasonable results for $`N=1`$. The case $`N=1`$ is also of special interest since an order parameter $`\overline{\psi }_j\psi ^j0`$ leads to a ground state which does not admit any discrete symmetry involving the reflection of all coordinates, in contrast to the models with $`N2`$.

In this letter we improve the exact renormalization group approach and establish a second order phase transition for $`N=1`$. We also compute the critical exponents. This is important beyond a possible relevance for real physical systems: the GN model constitutes a peculiar universality class due to the presence of massless fermions at the critical point. Just as the $`O(N)`$-Heisenberg models for bosons, the GN model could in the future become a benchmark for our understanding of critical systems in presence of fermions.

The GN model describes $`N`$ fermionic fields with local quartic interaction. Here $`\psi _j,j=1\mathrm{}N`$, are irreducible representations of the group $`O(d)`$ including parity reflections, i. e. $`2^{d/2}`$ component Dirac spinors for $`d`$ even and $`2^{(d1)/2}`$ for $`d`$ odd. The classical Euclidean action

$$S=\mathrm{d}^dx\left\{\overline{\psi }_j(x)i/\psi ^j(x)+\frac{\overline{G}}{2}\left(\overline{\psi }_j(x)\psi ^j(x)\right)^2\right\}.$$
(1)
is symmetric under a coordinate reflection $`\psi (x)\psi (x),\overline{\psi }(x)\overline{\psi }(x)`$. (We note that $`\psi `$ and $`\overline{\psi }`$ are independent variables in an Eucidean formulation.) A nonvanishing expectation value of $`\overline{\psi }_j\psi ^j`$ spontaneously breaks this symmetry. If the spinors $`\psi _j`$ contain more than one irreducible representation of the rotation (or Lorentz) group $`SO(d)`$ we can find alternative definitions of the coordinate reflection where $`\overline{\psi }_j\psi ^j`$ remains invariant. In particular, this is realized for Dirac spinors in even dimensions where a mass term couples two irreducible representations $`\psi _L`$ and $`\psi _R`$. One may then define a standard parity transformation under which $`\overline{\psi }\psi `$ is invariant. In this case, however, $`\overline{\psi }_R\psi _L0`$ breaks a discrete chiral symmetry $`\psi _L\psi _L,\overline{\psi }_L\overline{\psi }_L`$. We will concentrate here on one two-component spinor in three dimensions. For $`N=1`$ the above “parity transformation” is the only possible choice of coordinate reflections and $`\overline{\psi }_j\psi ^j0`$ spontaneously breaks this symmetry.

Depending on the value of the coupling $`\overline{G}`$ the model will be in the symmetric phase or exhibit spontaneous symmetry breaking (SSB) with a nontrivial expectation value of the order parameter $`G\overline{\psi }_j\psi ^j`$. We will describe a space dependent fermion bilinear $`i\overline{G}\overline{\psi }_j(x)\psi ^j(x)`$ by a real scalar field $`\sigma (x)`$ such that SSB is indicated by nonzero homogeneous $`\sigma 0`$. By a partial bosonization we express the GN model (1) as an equivalent Yukawa model with

$$S_\sigma =\mathrm{d}^dx\left\{\overline{\psi }_ji/\psi ^j+i\sigma \overline{\psi }_j\psi ^j+\frac{1}{2\overline{G}}\sigma ^2\right\}.$$
(2)
The equivalence of the partition function can be seen by performing the Gaussian $`\sigma `$-integration ($`\overline{\eta }\psi `$ means $`\mathrm{d}^dx\overline{\eta }_j(x)\psi ^j(x)`$ etc.)

$`Z[\overline{\eta },\eta ]`$ $`=`$ $`{\displaystyle 𝒟\sigma 𝒟\psi 𝒟\overline{\psi }\mathrm{exp}(S_\sigma [\sigma ,\psi ,\overline{\psi }]+\overline{\eta }\psi +\eta \overline{\psi })}`$ (3)
$`=`$ $`{\displaystyle 𝒟\psi 𝒟\overline{\psi }\mathrm{exp}\left(\overline{\psi }i/\psi +\overline{\eta }\psi +\eta \overline{\psi }\frac{\overline{G}}{2}(\overline{\psi }\psi )^2\right)}`$ (5)
$`\times {\displaystyle }𝒟\sigma \mathrm{exp}({\displaystyle \frac{1}{2\overline{G}}}(\sigma +i\overline{G}\overline{\psi }\psi )^2),`$
where the last factor yields an irrelevant constant.

A powerful tool for non-perturbative examinations are exact renormalization group equations for the effective action $`\mathrm{\Gamma }`$, i. e. the generating functional of the 1PI Green’s functions . Starting with the classical action $`S`$ at the UV cutoff $`\mathrm{\Lambda }`$, we obtain a type of coarse-grained free energy $`\mathrm{\Gamma }_k`$ by integrating out the quantum fluctuations with momenta larger than a given scale $`k`$. Eventually, we reach the macroscopic thermodynamic potential  $`\mathrm{\Gamma }`$ at $`k0`$. The IR cutoff $`k`$ is implemented in a smooth way by introducing a mass-like term $`\mathrm{\Delta }S_k`$ into the classical action which gives extra masses to modes with momenta smaller than $`k`$. In the limit $`k0`$ the IR cutoff is absent, and fluctuations at all scales have been taken into account. At $`k\mathrm{\Lambda }`$ all fluctuations are suppressed and $`\mathrm{\Gamma }_\mathrm{\Lambda }`$ approaches $`S`$.

This procedure should be explained more precisely. For notation purposes we combine the fermionic and (real) bosonic fields to a column vector $`\chi =(\sigma ,\psi ,\overline{\psi }^T)`$, and the row $`𝒥`$ contains all the external sources $`𝒥=(J,\overline{\eta },\eta ^T)`$. We start with the generating functional for the connected correlation functions in the presence of the IR cutoff:

$$W_k[𝒥]=\mathrm{ln}𝒟\chi \mathrm{exp}\left(S_\sigma [\chi ]\mathrm{\Delta }S_k[\chi ]+𝒥\chi \right)$$
(6)
For vanishing $`\mathrm{\Delta }S_k`$ (at $`k0`$) this matches exactly with the free energy of the Gross-Neveu model (1). The effective action is then defined via a modified Legendre transformation by

$`\mathrm{\Gamma }_k[\mathrm{\Phi }]:=W_k[𝒥[\mathrm{\Phi }]]+𝒥[\mathrm{\Phi }]\mathrm{\Phi }\mathrm{\Delta }S_k[\mathrm{\Phi }],`$ (7)
which depends on the expectation values of the fields $`\mathrm{\Phi }=\chi `$.

The infrared cutoff takes the form $`\mathrm{\Delta }S_k[\mathrm{\Phi }]=\frac{1}{2}\mathrm{\Phi }^TR_k\mathrm{\Phi }`$, and $`R_k`$ is a matrix:

$$R_k(p,q)=\left(\begin{array}{ccc}R_{kB}(q)& 0& 0\\ 0& 0& R_{kF}^T(q)\\ 0& R_{kF}(q)& 0\end{array}\right)(2\pi )^d\delta ^d(pq)$$

with $`R_{kB}(q)=Z_{\sigma ,k}q^2r_{kB}(q)`$ and $`R_{kF}(q)=Z_{\psi ,k}q/r_{kF}(q)`$. ($`Z_{\sigma ,k}`$ and $`Z_{\psi ,k}`$ are wavefunction renormalizations.) With these definitions, an exact renormalization group equation for the scale dependence of $`\mathrm{\Gamma }_k`$ can be found :

$$_t\mathrm{\Gamma }_k=\frac{1}{2}\mathrm{STr}\left\{(\mathrm{\Gamma }_k^{(2)}+R_k)^1_tR_k\right\}$$
(8)
where $`t=\mathrm{ln}(k/\mathrm{\Lambda })`$. The super-trace runs over momenta and all internal indices and provides appropriate minus signs for the fermionic sector. The heart of the flow equation is the fluctuation matrix

$$\left(\mathrm{\Gamma }_k^{(2)}(p,q)\right)_{ab}:=\frac{\stackrel{\stackrel{}{}}{\delta }}{\delta \mathrm{\Phi }_a^T(p)}\mathrm{\Gamma }_k\frac{\stackrel{\stackrel{}{}}{\delta }}{\delta \mathrm{\Phi }_b(q)}.$$
(9)
Together with $`R_k`$ it represents the exact inverse propagator at the scale $`k`$.

Equation (8) is an exact but complicated functional differential equation. There is no way around some approximation by truncating the most general form of $`\mathrm{\Gamma }_k`$. We work here in the lowest order of a systematic derivative expansion where $`\mathrm{\Gamma }_k`$ contains a scalar potential, kinetic terms and a Yukawa coupling. In momentum space it is given by ($`dq=\mathrm{d}^dq/(2\pi )^d`$):

$`\mathrm{\Gamma }_k[\sigma ,\psi ,\overline{\psi }]={\displaystyle }\mathrm{d}^dxU_k(\sigma )+{\displaystyle }\mathrm{d}q\{{\displaystyle \frac{Z_{\sigma ,k}}{2}}\sigma (q)q^2\sigma (q)`$ (10)
$`Z_{\psi ,k}\overline{\psi }_j(q)q/\psi ^j(q)+{\displaystyle }\mathrm{d}pi\overline{h}_k\overline{\psi }_j(p)\sigma (pq)\psi ^j(q)\}`$ (11)
The connection between $`\overline{G}`$ and $`\overline{h}`$ becomes clear if we rescale $`\sigma `$ in (2) to $`\overline{h}\sigma `$ and set $`\overline{G}=\overline{h}^2/m_\sigma ^2`$, $`m_\sigma `$ denoting the boson mass.

Using a truncation of the effective action causes the limit $`k0`$ to depend on the precise form of the cutoff functions. In order to take control of this we have used two different choices, an exponential  and a linear  cutoff:

$$yr_B^{exp}(y)=\frac{y}{\mathrm{exp}(y)1},yr_B^{lin}(y)=(1y)\mathrm{\Theta }(1y),$$

where $`y=q^2/k^2`$, and $`r_F(y)`$ is chosen in both cases such that $`y(1+r_B)=y(1+r_F)^2`$. We introduce renormalized, dimensionsless quantities $`h_k^2=Z_\sigma ^1Z_\psi ^2k^{d4}\overline{h}_k^2`$, $`\stackrel{~}{\rho }=\frac{1}{2}Z_\sigma k^{2d}\sigma ^2`$, $`u_k=U_kk^d`$, and we use $`u_k^{}=\frac{u_k}{\stackrel{~}{\rho }}`$ etc.

We obtain a set of evolution equations for the effective parameters of the theory by inserting (11) in the exact flow equation (8). The evolution equation of the effective scalar potential $`u_k`$ is found by evaluating $`\mathrm{\Gamma }_k^{(2)}`$ for a constant scalar background field:

$`_tu_k`$ $`=`$ $`du_k+(d2+\eta _\sigma )\stackrel{~}{\rho }u_k^{}`$ (12)
$`+`$ $`2v_d\left\{l_0^d(u_k^{}+2\stackrel{~}{\rho }u_k^{\prime \prime };\eta _\sigma )N^{}l_0^{(F)d}(2h_k^2\stackrel{~}{\rho };\eta _\psi )\right\},`$ (13)
where we use $`N^{}=d_\gamma N`$, with $`d_\gamma `$ the dimension of the $`\gamma `$ matrices, and $`v_d^1=2^{d+1}\pi ^{d/2}\mathrm{\Gamma }(d/2)`$. The definition of the threshhold functions $`l_n^d`$ and $`l_n^{(F)d}`$ as well as $`l_{n_1,n_2}^{(FB)d}`$, $`m_{n_1,n_2}^d`$, $`m_{2/4}^{(F)d}`$, and $`m_{n_1,n_2}^{(FB)d}`$ used below (see (24)–(27)) is given in the appendix of . These functions contain all the momentum integrations. For the linear cutoff the integrations can be performed analytically—an enormous advantage for the subsequent numerics. We obtain:

$`l_n^d(\omega ;\eta _\sigma )={\displaystyle \frac{2(\delta _{n,0}+n)}{d}}\left(1{\displaystyle \frac{\eta _\sigma }{d+2}}\right){\displaystyle \frac{1}{(1+\omega )^{n+1}}}`$ (14)
$`l_n^{(F)d}(\omega ;\eta _\psi )={\displaystyle \frac{2(\delta _{n,0}+n)}{d}}\left(1{\displaystyle \frac{\eta _\psi }{d+1}}\right){\displaystyle \frac{1}{(1+\omega )^{n+1}}}`$ (15)
$`l_{n_1,n_2}^{(FB)d}(\omega _1,\omega _2;\eta _\psi ,\eta _\sigma )={\displaystyle \frac{2}{d}}{\displaystyle \frac{1}{(1+\omega _1)^{n_1}(1+\omega _2)^{n_2}}}\times `$ (16)
$`\times \left\{{\displaystyle \frac{n_1}{1+\omega _1}}\left(1{\displaystyle \frac{\eta _\psi }{d+1}}\right)+{\displaystyle \frac{n_2}{1+\omega _2}}\left(1{\displaystyle \frac{\eta _\sigma }{d+2}}\right)\right\}`$ (17)
$`m_{n_1,n_2}^d(\omega _1,\omega _2;\eta _\sigma )={\displaystyle \frac{1}{(1+\omega )^{n_1}(1+\omega )^{n_2}}}`$ (18)
$`m_2^{(F)d}(\omega ;\eta _\psi )={\displaystyle \frac{1}{(1+\omega )^4}},m_4^{(F)d}(\omega ;\eta _\psi )={\displaystyle \frac{1}{(1+\omega )^4}}+`$ (19)
$`+{\displaystyle \frac{1\eta _\psi }{d2}}{\displaystyle \frac{1}{(1+\omega )^3}}\left({\displaystyle \frac{1\eta _\psi }{2d4}}+{\displaystyle \frac{1}{4}}\right){\displaystyle \frac{1}{(1+\omega )^2}}`$ (20)
$`m_{n_1,n_2}^{(FB)d}(\omega _1,\omega _2;\eta _\psi ,\eta _\sigma )=\left(1{\displaystyle \frac{\eta _\sigma }{d+1}}\right){\displaystyle \frac{1}{(1+\omega _1)^{n_1}(1+\omega _2)^{n_2}}}`$ (21)
The anomalous dimensions $`\eta _\sigma `$ and $`\eta _\psi `$ are defined as

$$\eta _\sigma (k)=_t\mathrm{ln}Z_{\sigma ,k},\eta _\psi (k)=_t\mathrm{ln}Z_{\psi ,k}.$$
(22)
Taking second derivatives of (8) with respect to the fields by means of (9), we obtain evolution equations for the exact inverse boson and fermion propagators. Expanding the momentum dependence at $`q^2=0`$ for $`\sigma =\sigma _{0k}`$ at the potential minimum yields equations for $`\eta _\sigma `$ and $`\eta _\psi `$:

$`\eta _\sigma (k)`$ $`=`$ $`8{\displaystyle \frac{v_d}{d}}\{\kappa _k(3u_k^{\prime \prime }+2\kappa _ku_k^{\prime \prime \prime })^2m_{4,0}^d(u_k^{}+2\kappa _ku_k^{\prime \prime },0;\eta _\sigma )`$ (23)
$`+N^{}`$ $`h_k^2`$ $`[m_4^{(F)d}(2h_k^2\kappa _k;\eta _\psi )2h_k^2\kappa _km_2^{(F)d}(2h_k^2\kappa _k;\eta _\psi )]\}`$ (24)
$`\eta _\psi (k)`$ $`=`$ $`8{\displaystyle \frac{v_d}{d}}h_k^2m_{1,2}^{(FB)d}(2h_k^2\kappa _k,u_k^{}+2\kappa _ku_k^{\prime \prime };\eta _\psi ,\eta _\sigma )`$ (25)
Here $`\kappa _k`$ denotes the position of the minimum of $`u_k(\stackrel{~}{\rho })`$ and all derivatives of the potential are evaluated at $`\stackrel{~}{\rho }=\kappa _k`$. Dividing the evolution equation of the fermion propagator by $`\sigma `$ and evaluating it at zero momentum and $`\sigma =\sigma _{0k}`$ yields the evolution equation of $`h_k^2`$:

$`_th_k^2`$ $`=`$ $`(2\eta _\psi +\eta _\sigma +d4)h_k^2`$ (26)
$`+`$ $`8v_dh_k^4l_{1,1}^{(FB)d}(2h_k^2\kappa _k,u_k^{}+2\kappa _ku_k^{\prime \prime };\eta _\psi ,\eta _\sigma )`$ (27)
We study the flow of the effective potential in two different ways that correspond to two different truncations of $`\mathrm{\Gamma }_k`$. First, we expand the evolution equation of the potential (13) in a Taylor series up to third order in $`\stackrel{~}{\rho }`$ around its minimum and evolve the resulting coupled ordinary differential equations . In a second, more involved approach, the full equation for $`u_k^{}`$ is discretized on a grid, yielding coupled ordinary differential equations (one for each grid point) which are solved simultaneously. We use the equation for $`u_k^{}`$ instead of the one for $`u_k`$ for numerical reasons.

The only critical parameter in the theory is the four-fermion coupling $`\overline{G}_{\mathrm{\Lambda }cr}`$. This parameter has to be tuned in order to be near the phase transition. The flow equations (13)–(27) are evolved from an UV cutoff scale $`\mathrm{\Lambda }`$ to $`k0`$. We concentrate on the three-dimensional case with one fermion species ($`d=3`$ and $`N=1`$). The initial values of the parameters are chosen such that $`\mathrm{\Gamma }_\mathrm{\Lambda }=S_\sigma `$: $`Z_{\sigma \mathrm{\Lambda }}=10^{10}0`$, $`Z_{\psi \mathrm{\Lambda }}=1`$, $`\overline{h}_\mathrm{\Lambda }^2=c\mathrm{\Lambda }`$ and $`u_\mathrm{\Lambda }^{}(\stackrel{~}{\rho })=e:=(Z_{\sigma \mathrm{\Lambda }}\overline{G}_\mathrm{\Lambda }\mathrm{\Lambda })^1`$ for all $`\stackrel{~}{\rho }`$. The GN model corresponds to $`c=1`$ whereas for $`c1`$ we investigate Yukawa type theories in the same universality class.

Very close to the phase transition we find scaling solutions for all evolving parameters corresponding to vanishing beta functions. This behaviour indicates a second order phase transition. The symmetric regime is characterized by $`u_k^{}(0)>0`$ as well as $`\kappa _k=0`$, while after spontaneous symmetry breaking $`u_k^{}(0)`$ becomes negative and $`\kappa _k>0`$. Starting in the symmetric regime the system evolves into the SSB regime and reaches the scaling solution. Further down the flow, it evolves back into the symmetric regime for $`\overline{G}_\mathrm{\Lambda }<\overline{G}_{\mathrm{\Lambda }_{cr}}`$ or it remains in the SSB regime for $`\overline{G}_\mathrm{\Lambda }>\overline{G}_{\mathrm{\Lambda }_{cr}}`$. This behaviour contrasts the one for $`N2`$ where the scaling solution is located in the symmetric regime.

Figure 1 shows the critical flow of $`\kappa _k`$ and $`u_k^{}(0)`$ for two different values of $`c`$, namely $`c=1`$ (full and dashed line, resp.) and $`c=10^{10}`$ (thin dashed lines). It nicely demonstrates that only the beginning of the flow is affected by the choice of the non-critical parameters. All universal quantities like critical exponents and mass gaps are independent of $`c`$. The polynomial expansion of the potential works well only for small enough $`c`$. The value of the critical coupling $`\overline{G}_{\mathrm{\Lambda }cr}`$ depends, however, on the choice of $`c`$. In the limit $`c1`$ it converges towards a constant value, which is reported in table I.

We calculate the critical exponents characterizing a second order phase transition. In order to test the reliability of the numerical algorithms, we determine the exponents in the symmetric as well as in the ordered phase. The scale dependent renormalized boson mass $`m_{\sigma R}^2(k)`$ is defined as the curvature of the potential at its minimum. In the SSB regime it is given by $`m_{\sigma R}^2(k)=2k^2\kappa _ku_k^{\prime \prime }(\kappa _k)`$, in the symmetric regime by $`m_{\sigma R}^2(k)=k^2u_k^{}(0)`$. In the ordered phase the running of $`m_{\sigma R}(k)`$ essentially stops once $`k`$ becomes much smaller than $`m_{\sigma R}(k)`$. The situation is slightly more difficult in the symmetric phase since the fermions are massless. Their fluctuations induce a scale dependence of $`Z_{\sigma ,k}`$ even for very small momenta $`k`$: $`\eta _\sigma 1`$ for $`k0`$ in contrast to $`\eta _\sigma 0`$ for $`k0`$ in the ordered phase. To get rid of this problem we define the renormalized boson mass at some fixed small ratio $`r_c=k/\overline{m}_{\sigma R}`$  as

$$\overline{m}_{\sigma R}^2=k_c^2\left(u_{k_c}^{}(0)u_{k_c}^{cr}(0)\right),k_c=r_c\overline{m}_{\sigma R}.$$
(28)
where $`u_k^{cr}(0)`$ denotes $`u_k^{}(0)`$ on the critical trajectory. In the numerical calculations we have used a ratio $`r_c=0.01`$, but our results do not depend on $`r_c`$ for $`r_c0.1`$. Thus we define the critical exponent $`\nu `$

$`\nu ={\displaystyle \frac{1}{2}}\underset{\delta e0}{lim}{\displaystyle \frac{\mathrm{ln}\overline{m}_{\sigma R}^2(\delta e)}{\mathrm{ln}\delta e}},`$ (29)
where $`\delta e=|ee_{cr}|`$. The exponent $`\gamma `$ is defined as usual, since the unrenormalized boson mass $`m_\sigma ^2=Z_\sigma m_{\sigma R}^2`$ is not affected by the fluctuations of the massless fermions:

$`\gamma =\underset{\delta e0}{lim}{\displaystyle \frac{\mathrm{ln}m_\sigma ^2(\delta e)}{\mathrm{ln}\delta e}}`$ (30)
In the ordered phase both exponents are defined as usual, using $`m_{\sigma R}^2`$ instead of $`\overline{m}_{\sigma R}^2`$ for the definition of $`\nu `$. The critical exponent $`\beta `$ is defined as

$$\beta =\frac{1}{2}\underset{\delta e0}{lim}\frac{\mathrm{ln}\sigma _0^2}{\mathrm{ln}\delta e},$$
(31)
with $`\sigma _0=lim_{k0}\sigma _{0k}`$. The exponents $`\eta _\sigma `$ and $`\eta _\psi `$ for the critical correlation functions are computed by taking the values of the scale dependent anomalous dimensions $`\eta _\sigma (k)`$ and $`\eta _\psi (k)`$ at the scaling solution (sec. 4.2 of ). Table I lists our results for the critical exponents. We find a good match of the values in the two different phases. Besides, we have checked the index relations $`\gamma =\nu (2\eta _\sigma )`$ and $`\beta =\frac{1}{2}\nu (d2+\eta _\sigma )`$ which are well fulfilled. The dependence of the exponents on the cutoff functions $`r_B`$ and $`r_F`$ as well as on the truncation can also be seen from the table. The latter one seems to be weaker when using the linear cutoff. The error on the exponent could be larger than the difference between the results for different truncations of the potential and different cutoffs. This issue could be investigated by extending the truncation (11), e. g. by including a quartic fermion interaction $`(\overline{\psi }\psi )^2`$ along the lines discussed in .

Since the fluctuations generate bosonic and fermionic masses, one might be interested in the resulting gaps. They are proportional to the order parameter $`\rho _0=\frac{1}{2}Z_\sigma \sigma _0^2`$:

$$m_{\sigma R}=\mathrm{\Delta }_\sigma \rho _0,m_{\psi R}=\mathrm{\Delta }_\psi \rho _0$$
(32)
The gaps $`\mathrm{\Delta }_\sigma `$ and $`\mathrm{\Delta }_\psi `$ are shown in table I.

| Truncation | full eq. for $`u_k^{}`$ | | Taylor expansion | |
| --- | --- | --- | --- | --- |
| Cutoff | lin | exp | lin | exp |
| $`\nu _{symm}`$ | 0.621 | 0.640 | 0.623 | 0.633 |
| $`\gamma _{symm}`$ | 1.051 | 1.077 | 1.053 | 1.062 |
| $`\nu (2\eta _\sigma )`$ | 1.051 | 1.076 | 1.054 | 1.064 |
| $`\nu _{ssb}`$ | 0.620 | 0.637 | 0.622 | 0.632 |
| $`\gamma _{ssb}`$ | 1.050 | 1.071 | 1.053 | 1.062 |
| $`\beta _{ssb}`$ | 0.406 | 0.420 | 0.407 | 0.417 |
| $`\frac{\nu }{2}(1+\eta _\sigma )`$ | 0.405 | 0.420 | 0.407 | 0.417 |
| $`\eta _\sigma `$ | 0.308 | 0.319 | 0.308 | 0.319 |
| $`\eta _\psi `$ | 0.112 | 0.114 | 0.112 | 0.113 |
| $`\mathrm{\Delta }_\sigma `$ | 16.0 | 17.6 | 16.8 | 18.1 |
| $`\mathrm{\Delta }_\psi `$ | 14.2 | 14.9 | 14.4 | 15.2 |
| $`\overline{G}_{\mathrm{\Lambda }cr}\mathrm{\Lambda }`$ | 43.13 | 26.68 |  |  |

TABLE I. Critical exponents and mass gaps for $`N=1`$ in three dimensions

We have also calculated the critical exponents in the three-dimensional GN model for $`N>1`$. They are compared with results from other methods in table II. For $`N=2,4`$ and $`12`$ critical exponents have been calculated in the $`1/N`$-expansion to $`O(1/N^2)`$  and anomalous dimensions to $`O(1/N^3)`$ using conformal techniques . Monte-Carlo methods have been used to calculate critical exponents for $`N=4`$ which are compared with results from the $`ϵ=4d`$ expansion to $`O(ϵ^2)`$. Our values for the Taylor expansion with exponential cutoff (b) agree well with  where the precise numerical implementation was different.

The largest discrepancy between different approaches concerns the values of the anomalous dimensions. This is similar as for the bosonic $`O(N)`$ models and is expected to improve substantially in the next order in the derivative expansion . The overall picture is, however, quite satisfactory and lends support to the validity of our method based on the exact renormalization group equation. We believe that our finding of a second order phase transition for $`N=1`$ is quite robust.

| N | 2 | 4 | 12 |
| --- | --- | --- | --- |
|  | 0.927<sup>a</sup> | 1.018<sup>a</sup> | 1.018<sup>a</sup> |
| $`\nu `$ | 0.962<sup>b</sup> | 1.016<sup>b</sup> | 1.011<sup>b</sup> |
|  | 0.738<sup>c</sup> | 0.903<sup>c</sup> | 1.007<sup>c</sup> |
|  |  | 0.948<sup>d</sup> |  |
|  |  | 1.00(4)<sup>e</sup> |  |
|  | 0.525<sup>a</sup> | 0.756<sup>a</sup> | 0.927<sup>a</sup> |
| $`\eta _\sigma `$ | 0.554<sup>b</sup> | 0.786<sup>b</sup> | 0.935<sup>b</sup> |
|  | 0.635<sup>c</sup> | 0.776<sup>c</sup> | 0.914<sup>c</sup> |
|  |  | 0.763<sup>d</sup> | 0.913<sup>f</sup> |
|  |  | 0.754(8)<sup>e</sup> |  |
|  | 0.071<sup>a</sup> | 0.032<sup>a</sup> | 0.0087<sup>a</sup> |
| $`\eta _\psi `$ | 0.067<sup>b</sup> | 0.028<sup>b</sup> | 0.0057<sup>b</sup> |
|  | 0.105<sup>f</sup> | 0.044<sup>f</sup> | 0.0124<sup>f</sup> |
| $`\mathrm{\Delta }_\sigma `$ | 15.1<sup>a</sup> | 12.0<sup>a</sup> | 5.5<sup>a</sup> |
|  | 17.3<sup>b</sup> | 14.0<sup>b</sup> | 6.4<sup>b</sup> |
| $`\mathrm{\Delta }_\psi `$ | 10.8<sup>a</sup> | 7.5<sup>a</sup> | 3.3<sup>a</sup> |
|  | 11.4<sup>b</sup> | 7.9<sup>b</sup> | 3.4<sup>b</sup> |

TABLE II. Critical exponents for various $`N`$ in three dimensions
(a) linear cutoff, full equation for $`u_k^{}`$
(b) exponential cutoff, Taylor expansion of $`u_k^{}`$
(c) $`1/N`$-expansion (Padé-Borel resummed)
(d) $`ϵ`$-expansion (Padé-Borel resummed)
(e) Monte-Carlo simulations
(f) conformal techniques