text:
# Super Yang-Mills on the lattice with domain wall fermions

## 1 Introduction

It is believed that super-symmetric (SUSY) field theories may play an important role in describing the physics beyond the Standard Model. Non-perturbative studies of these theories are of great interest. First-principles numerical simulations may be able to provide additional information and confirmation of existing analytical calculations. Typically first principles numerical simulations of field theories are done within the framework of the lattice regulator. A host of results have been produced in this way for many field theories and most notably QCD. Several SUSY theories can also be formulated on the lattice and be studied numerically. To be more specific consider the problems of putting a SUSY theory on the lattice (see for example ):

1) Since space-time is discrete only a subgroup of the Poincaré group survives and as a result SUSY is broken. This problem is not severe and is of the same nature as in QCD. The symmetry breaking operators that are allowed by the remaining symmetries are irrelevant. One can calculate at several lattice spacings $`a`$ and then take the $`a0`$ limit. No fine tuning is needed.

2) If the SUSY theory under consideration involves scalar fields one can have scalar mass terms that break SUSY since typically they are not forbidden by some symmetry. Since these operators are relevant fine tuning will be needed in order to cancel their contributions. The four-dimensional $`𝒩=1`$ Super Yang-Mills (SYM) theory does not involve scalars and therefore it does not have this problem.

3) A naive regularization of fermions results in the well known doubling problem . For each fermion species in the four-dimensional continuum 16 are generated on the lattice with total chirality of zero. This results in the wrong number of degrees of freedom and therefore breaks SUSY. However, this problem may be possible to treat as in QCD. This is the case for $`𝒩=1`$ SYM.

One possible way to remove the unwanted fermion degrees of freedom is to add an irrelevant operator (Wilson term ) that gives them heavy masses of the size of the cutoff. This term unavoidably breaks the chiral symmetry  and as a result a gluino mass term is no longer forbidden. Since such a term is relevant, fine tuning of the bare fermion mass is necessary as the continuum limit is approached in order to cancel its contribution. Although fine tuning is not a welcomed property this method makes it possible to recover the continuum target theory.

Therefore, it is possible to simulate numerically the $`𝒩=1`$ SYM theory using existing lattice “technology” since all three difficulties can be circumvented. This observation was made some time ago . In particular, it was argued that, using a standard lattice gauge theory action with a pure gauge Wilson plaquette term and Wilson fermions in the adjoint representation, numerical simulations could be done. Pioneering work using these methods has already produced very interesting numerical results . Also, for proposed lattice tests of SYM see . For a supersymmetric formulation on the lattice using Kogut–Susskind  fermions see .

There are two unwelcomed difficulties in using Wilson fermions. The first has already been mentioned and it is the need for fine tuning. The second is of a technical nature. It turns out that the Pfaffian resulting from the fermionic integration is not positive definite  at finite lattice spacing. However, it does become positive definite as the continuum limit is approached and therefore as a “cure” only the absolute value of the Pfaffian is used . However, this introduces non-analyticities that may make the approach to the continuum limit difficult.

Both of these difficulties can be brought under control by using an alternative fermion lattice regulator, domain wall fermions (DWF). The use of DWF in supersymmetric theories has been explored in the very nice work of . The methods in this paper are along the lines of these references. Domain wall fermions were introduced in , were further developed in  and in . They provide a new way for treating the unwelcomed chiral symmetry breaking that is introduced when the fermion doubler species are removed. Here a variant of this approach will be used . For reviews on the subject please see  and references therein. DWF have already been used for numerical simulations of the two flavor dynamical Schwinger model , dynamical QCD , quenched QCD , as well as for simulations of 4-Fermi models . The use of DWF in supersymmetric theories has also been explored in a different fashion in . Furthermore, the use of overlap  type fermions has been explored in , and the use of other related types of fermions has been explored in .

In the lattice DWF formulation of a vector-like theory the fermionic fields are defined on a five dimensional space-time lattice using a local action. The fifth direction can be thought of as an extra space-time dimension or as a new internal flavor space. The gauge fields are introduced in the standard way in the four dimensional space-time and are coupled to the extra fermion degree of freedom in a diagonal fashion. The key ingredient is that the boundary conditions of the Dirac operator along the fifth direction are taken to be free. As a result, although all fermions are heavy, two chiral, exponentially bound surface states appear on the boundaries (domain walls) with the plus chirality localized on one wall and the minus chirality on the other. The two chiralities mix only by an amount that is exponentially small in $`L_s`$, where $`L_s`$ is the number of lattice sites along the fifth direction, and form a Dirac spinor that propagates in the four-dimensional space-time with an exponentially small mass. Therefore, the amount of chiral symmetry breaking that is artificially induced by the regulator can be controlled by the new parameter $`L_s`$. In the $`L_s\mathrm{}`$ limit the chiral symmetry is exact, even at finite lattice spacing, so there is no need for fine-tuning.

For the first time the approach to the chiral limit has been separated from the approach to the continuum limit. Furthermore, the computing requirement is linear in $`L_s`$. This is to be contrasted with traditional lattice fermion regulators where the chiral limit is approached only as the continuum limit is taken, a process that is achieved at a large computing cost. Specifically, because of algorithmic reasons, the computing cost to reduce the lattice spacing by a factor of two grows by a factor of $`2^{810}`$ in four dimensions. Therefore, the unique properties of DWF provide a way to bring under control the systematic chiral symmetry breaking effects using today’s supercomputers.

The purpose of this paper is two-fold. First, the techniques for performing a numerical simulation of the full $`𝒩=1`$ SU(2) SYM theory using DWF are collected and it is demonstrated that they work as expected by performing numerical simulations of the full theory. Second, the gluino condensate is measured. It is expected that a non-zero gluino condensate must form . However, there are also arguments that the theory has a phase where a gluino condensate does not form . In the numerical simulations performed here it is found that a non-zero gluino condensate is sustained in the limit of zero gluino mass. This result is at a finite lattice spacing and therefore SUSY is still broken albeit by irrelevant operators.

It must be emphasized that in this work due to limited computer resources no attempt has been made to extrapolate to the continuum limit. It is possible that in this limit the gluino condensate may vanish. Future work using larger computer resources could calculate the gluino condensate at several lattice spacings and extract the continuum value. But even then, it will never be possible to numerically prove that the finite lattice spacing theory is not separated from the continuum theory by a phase transition. This problem is not particular to the case at hand but is in the nature of numerical investigations. They can provide strong evidence but not unquestionable proof. A well known case with similar problems relates to the question of confinement and chiral symmetry breaking in QCD.

This paper is organized as follows. In section 2, the DWF lattice formulation of $`𝒩=1`$ SU(2) SYM is presented. In section 3, analytical considerations relating to the gluino mass, the Ward identities and the effects of topology in the patterns of chiral symmetry breaking are given. The numerical methods used in the simulations are discussed in section 4. The numerical results are presented in section 5 and the paper is concluded in section 6.

## 2 Lattice formulation

In this section, the $`𝒩=1`$, SU(2) SYM lattice action and operators are presented. The approach is similar in spirit as to the case of Wilson fermions . The DWF formulation for this theory is identical to  and . It is presented below for the convenience of the reader and in order to establish notation.

The $`𝒩=1`$, SU(2) SYM theory is an SU(2) gauge theory with Majorana fermions in the adjoint representation. As such, the fermionic path integral results in the analytic square root of the corresponding Dirac determinant. This then is the Pfaffian of an antisymmetric matrix that has the same determinant as the Dirac operator. On the lattice, the Dirac operator can be defined using Wilson’s approach as in  or the DWF approach as in  and .

The partition function is:

$$Z=[dU][d\mathrm{\Psi }][d\mathrm{\Phi }]e^S$$
(1)
$`U_\mu (x)`$, $`\mu =1,2,3,4`$ is the four-dimensional gauge field in the fundamental representation, $`\mathrm{\Psi }(x,s)`$ is a (real) five-dimensional Majorana field in the adjoint representation and $`\mathrm{\Phi }(x,s)`$ is a (real) five-dimensional bosonic Pauli Villars (PV) type field with the same indices as the Majorana field. $`x`$ is the coordinate in the four-dimensional space-time box with extent $`L`$ along each of the four directions. The boundary conditions along these directions are taken to be periodic for all fields. The coordinate of the fifth direction is $`s=0,1,\mathrm{},L_s1`$, where $`L_s`$ is the size of that direction and is taken to be an even number. The action $`S`$ is given by:

$$S=S_G(U)+S_F(\mathrm{\Psi },U)+S_{PV}(\mathrm{\Phi },U).$$
(2)
$`S_G(U)`$ is the pure gauge part and is defined using the standard single plaquette action of Wilson:

$$S_G=\beta \underset{p}{}(1\frac{1}{2}\mathrm{ReTr}[U_p])$$
(3)
where $`\beta =4/g^2`$ and $`g`$ is the gauge coupling.

The fermion part $`S_F(\mathrm{\Psi },U)`$ is given by:

$$S_F=\underset{x,x^{},s,s^{}}{}\overline{\mathrm{\Psi }}(x,s)D_F(x,s;x^{},s^{})\mathrm{\Psi }(x^{},s^{})$$
(4)
where $`D_F`$ is the DWF Dirac operator in the form of . Specifically it is:

$$D_F(x,s;x^{},s^{})=\delta (ss^{})D/(x,x^{})+D/^{}(s,s^{})\delta (xx^{})$$
(5)
$`D/(x,x^{})`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu =1}{\overset{4}{}}}\left[(1+\gamma _\mu )V_\mu (x)\delta (x+\widehat{\mu }x^{})+(1\gamma _\mu )V_\mu ^{}(x^{})\delta (x^{}+\widehat{\mu }x)\right]`$ (6)
$`+`$ $`(m_04)\delta (xx^{})`$
$$D/^{}(s,s^{})=\{\begin{array}{cc}P_R\delta (1s^{})m_fP_L\delta (L_s1s^{})\delta (0s^{})\hfill & s=0\hfill \\ P_R\delta (s+1s^{})+P_L\delta (s1s^{})\delta (ss^{})\hfill & 0<s<L_s1\hfill \\ m_fP_R\delta (0s^{})+P_L\delta (L_s2s^{})\delta (L_s1s^{})\hfill & s=L_s1\hfill \end{array}$$
(7)
$$P_{R,L}=\frac{1\pm \gamma _5}{2}$$
(8)
where $`V`$ is the gauge field in the adjoint representation. It is related to the field in the fundamental representation by (see for example ):

$$[V_\mu (x)]_{a,b}=2\mathrm{T}\mathrm{r}[U_\mu ^{}(x)T^aU_\mu (x)T^b]$$
(9)
and

$$V_\mu (x)=V_\mu ^{}(x)=[V_\mu ^1(x)]^T$$
(10)
where $`T^a=\frac{1}{2}\sigma ^a`$ with $`\sigma ^a`$ the Pauli matrices. In the above equations $`m_0`$ is a five-dimensional mass representing the “height” of the domain wall and it controls the number of light flavors in the theory. In order to get one light species in the free theory one must set $`0<m_0<2`$ . The parameter $`m_f`$ explicitly mixes the two chiralities and as a result it controls the bare fermion mass of the four-dimensional effective theory. The dependence of the bare fermion mass on $`m_0`$ and $`L_s`$ is discussed in section 3.1.

The fermion field $`\overline{\mathrm{\Psi }}`$ is not independent but is related to $`\mathrm{\Psi }`$ by the equivalent of the Majorana condition for this 5-dimensional theory :

$$\overline{\mathrm{\Psi }}=\mathrm{\Psi }^TCR_5$$
(11)
where $`R_5`$ is a reflection operator along the fifth direction and $`C`$ the charge conjugation operator in Eucledean space which can be set to:

$$C=\gamma _0\gamma _2.$$
(12)
Therefore, the fermion action can also be written as:

$$S_F=\underset{x,x^{},s,s^{}}{}\mathrm{\Psi }^T(x,s)M_F(x,s;x^{},s^{})\mathrm{\Psi }(x^{},s^{})$$
(13)
where

$$M_F(x,s;x^{},s^{})=CR_5D_F(x,s;x^{},s^{})$$
(14)
is an antisymmetric matrix as can be easily checked . As a result the fermionic integral gives the anticipated Pfaffian:

$$[d\mathrm{\Psi }]e^{S_F}=\mathrm{Pf}(M_F).$$
(15)
Because $`det(CR_5)=1`$ one also has that $`det(M_F)=det(D_F)`$ and therefore:

$$\mathrm{Pf}(M_F)=\sqrt{det(D_F)}.$$
(16)
The Pauli-Villars action $`S_{PV}`$ is designed to cancel the contribution of the heavy fermions . Viewing the extra dimension as an internal flavor space  one can see that there are $`L_s1`$ heavy fermions with masses near the cutoff and one light fermion. The PV subtraction subtracts the $`L_s`$ heavy particles. As was pointed in  this amounts to a “double” regularization of the light degree of freedom, first by the lattice and then by the PV field. The form of the PV subtraction used here is as in  and is given by:

$$S_{PV}=\underset{x,x^{},s,s^{}}{}\mathrm{\Phi }^T(x,s)M_F[m_f=1](x,s;x^{},s^{})\mathrm{\Phi }(x^{},s^{}).$$
(17)
The integral over the PV fields results in:

$$[d\mathrm{\Phi }]e^{S_{PV}}=\frac{1}{\mathrm{Pf}(M_F[m_f=1])}.$$
(18)
Green functions in this work are measured using four-dimensional fermion fields constructed from five-dimensional fermion fields using the projection prescription :

$`\chi (x)`$ $`=`$ $`P_R\mathrm{\Psi }(x,0)+P_L\mathrm{\Psi }(x,L_s1)`$
$`\overline{\chi }(x)`$ $`=`$ $`\overline{\mathrm{\Psi }}(x,L_s1)P_R+\overline{\mathrm{\Psi }}(x,0)P_L.`$ (19)
In the $`L_s\mathrm{}`$ limit of the theory these operators directly correspond to insertions in the overlap of appropriate creation and annihilation operators .

Using eq. 11 and 19 the Majorana condition on the four-dimensional fermion field is:

$$\overline{\chi }=\chi ^TC.$$
(20)
Because this is the correct condition for a four-dimensional field one can see that the definition in eq. 11 not only produces an antisymmetric fermion matrix $`M_F`$ but is also consistent with the projection prescription in eq. 19 as expected.

## 3 Analytical Considerations

In this section some analytical considerations are presented. In the $`𝒩=1`$ SYM theory, a gluino mass term is the only relevant operator that can break supersymmetry and is also the only relevant operator that can break (at the classical level) the $`U(1)_A`$ symmetry. Therefore, the two symmetries are intimately related to the mechanisms that can introduce a bare gluino mass term. These mechanisms depend on the “extra” regulator parameters $`m_0`$ and $`L_s`$. This is discussed below. Next the fate of the $`U(1)_A`$ chiral symmetry and the effects of topology are presented. The chiral and supersymmetric Ward identities are derived in the last subsection.

### 3.1 The “extra” DWF parameters

DWF introduce two extra parameters, the size of the fifth direction $`L_s`$ and the domain wall height or five-dimensional mass $`m_0`$. These two parameters together with the explicit mass $`m_f`$ control the bare fermion mass $`m_{\mathrm{eff}}`$. In the free theory one finds :

$$m_{\mathrm{eff}}=m_0(2m_0)\left[m_f+(1m_0)^{L_s}\right],0<m_0<2.$$
(21)
In the interacting theory one would expect that $`m_0`$ as well as its range of values will be renormalized. From the above equation one can see that for the free theory the value of $`m_0=1`$ is optimal in the sense that finite $`L_s`$ effects do not contribute to $`m_{\mathrm{eff}}`$. In the interacting theory one would expect that there is no such “optimal value” since, in a heuristic sense, $`m_0`$ will fluctuate. For a more detailed analysis please see . Then one would like $`L_s`$ to be large enough so that the second term in eq. 21 will be small allowing for simulations at reasonably small masses and/or for dependable extrapolations to the $`m_f0`$, $`L_s\mathrm{}`$ limit.

The effects of finite $`L_s`$ on the chiral symmetry can be best understood in the overlap formalism . In that formalism a transfer matrix $`T`$ along the extra direction is constructed. Because the gauge fields are not changing along that direction the product of transfer matrices simply results in $`T^{L_s}`$. For $`L_s=\mathrm{}`$ this is a projection operator that projects the reference vacuum state to a ground state. The fermion determinant is then the overlap of the reference vacuum state with that ground state. In  it was shown that, as a lattice gauge field configuration changes, from say the zero topological sector to sector one, an eigenvalue (or a degenerate set of eigenvalues) of the corresponding Hamiltonian $`H`$ changes sign. As a result, the filling level of the ground state becomes different from that of the reference vacuum state. Then the overlap is zero indicating the presence of an exact zero mode. This remarkable property is maintained to a good degree even at finite $`L_s`$ as was found in . Unfortunately, this property is also the reason for most of the difficulties with DWF. As the eigenvalue of the Hamiltonian $`H`$ changes sign it crosses zero. In such a configuration the transfer matrix has an eigenvalue equal to one and therefore even at $`L_s=\mathrm{}`$ there is no decay along the extra direction, the two chiralities do not decouple, and chiral symmetry can not be restored. Fortunately, configurations for which $`H`$ has an exact zero eigenvalue (for a given $`m_0`$) are of measure zero  and therefore are of no consequence. However, configurations in their neighborhood are not of measure zero and such configurations will exhibit very slow decay rates. Therefore, in order to restore chiral symmetry, very large values of $`L_s`$ may be needed. Since one would expect that the neighborhoods of such configurations are suppressed closer to the continuum limit this problem should become less severe as that limit is taken. This has been observed in numerical simulations of the Schwinger model , of full QCD , and of quenched QCD .

In the region where it makes sense to parameterize these effects by a residual mass in an effective action it has been found that:

$$m_{\mathrm{eff}}=c_0m_f+m_{\mathrm{res}},m_{\mathrm{res}}=c_1\mathrm{exp}(c_2L_s)$$
(22)
where for dynamical QCD at the currently accessible lattice spacings the decay is found to be $`c_20.02`$ . For quenched QCD the situation is better because current computing resources can simulate lattices with smaller lattice spacing. There, a value of $`c_20.1`$ is found . Also in these studies the value of $`c_2`$ was a weakly changing function of $`m_0`$ indicating that for practical purposes there is no optimal value of $`m_0`$.

In the case of the $`𝒩=1`$ SYM SU(2) theory the Hamiltonian corresponding to the five dimensional transfer matrix has eigenvalues that are doubly degenerate because the fermion fields are in the adjoint representation . Therefore when there is a “topology” change two eigenvalues will have to cross through zero (as compared to one for fundamental fermions). This may make this theory harder to study than QCD in the sense that larger $`L_s`$ values may be required. On the other hand, since no massless Goldstone particles are expected, the sensitivity of the spectrum on $`L_s`$ may be considerably milder. In any case, in this paper the only fermionic observable that will be discussed is the gluino condensate. This quantity is known to approach its $`L_s\mathrm{}`$ limit with faster decay rates than the ones in $`m_{\mathrm{res}}`$ (for a discussion and results for full QCD see ; there the decay rate for the chiral condensate was about five times faster than that for $`m_{\mathrm{res}}`$).

As was discussed above, the range of $`m_0`$ is renormalized by the interactions. It has been found that as the lattice spacing increases and one moves away from the continuum limit this range shrinks in size and for currently accessible spacings in QCD that range is about $`[1.4,2.0]`$. As one moves even farther away from the continuum limit this range can shrink to zero and then it will not be possible to have light DWF modes . However, it must be emphasized that for as long as the range of allowed values of $`m_0`$ is not of zero size the overlap formalism, although it does not specify how it is approached, guarantees the existence of the $`L_s\mathrm{}`$ limit. In this work, $`m_0=1.9`$ and, as it will be shown in section 5, the behavior of the gluino condensate vs. $`L_s`$ is consistent with an exponential ansatz.

### 3.2 Chiral symmetry and topology

Fermions in the adjoint representation of the SU(N) gauge group have a Dirac operator with index:

$$2N\nu $$
(23)
where $`\nu `$ is the winding of the background field configuration. Classical instantons have integer winding and they cause condensation of operators with $`2N`$ Majorana fermions. This results in the breaking of the $`U(1)_A`$ chiral symmetry down to the $`Z_{2N}`$ symmetry by the corresponding anomaly. The remaining $`Z_{2N}`$ symmetry may break spontaneously down to $`Z_2`$ . Mechanisms for this further breaking have been explored for example in  where instantons and fractionally charged objects such as torons  or caloron monopole constituents  were investigated as the source of this symmetry breaking.

Since in a toroidal geometry fractional winding numbers are possible , the partition function of the full theory can be expressed as:

$$Z(\theta )=\underset{\nu }{}e^{i\nu \theta }Z_\nu ,\nu =0,\pm 1/N,\pm 2/N,\mathrm{}$$
(24)
where $`\theta `$ is the vacuum angle and $`Z_\nu `$ is the partition function on the sector with winding $`\nu `$. For the theory with a soft breaking by a mass $`m_f`$ the interplay of the volume and mass in the formation of the gluino condensate has been analyzed in . The reader is referred to that reference for a very nice presentation on the subject. Assuming a mass gap is present in the theory the authors of  show that non-zero contributions to the gluino condensate $`\overline{\chi }\chi `$ come almost exclusively from the $`\nu =1/N`$ sector if $`m_f\times V\times \overline{\chi }\chi 1`$. On the other hand, if $`m_f\times V\times \overline{\chi }\chi 1`$ all sectors contribute to a non-zero condensate.

The above considerations result in an unusual picture. If the infinite volume limit is taken (followed by the massless limit) it is possible that a gluino condensate will form due to spontaneous breaking of the discrete symmetry $`Z_{2N}`$ down to $`Z_2`$. On the other hand, at a finite volume and zero mass a gluino condensate can form due to the presence of fractional winding configurations. Since the volume is finite, this can not be the result of spontaneous symmetry breaking. Instead, it is similar to symmetry breaking due to topological effects as, for example, in one flavor QCD. As pointed above the size of $`m_f\times V\times \overline{\chi }\chi `$ controls which “scenario” takes place.

On the lattice there is no clear definition of topology. The path integral over the SU(N) group space generates configurations of all possible windings. In order for the lattice theory to be able to reproduce phenomena that relate to topology it is essential that the lattice Dirac operator obeys the index theorem in a statistical sense. This is highly non–trivial since it is obviously related to the doubling problem. Traditional fermions (Wilson or staggered) do not exhibit exact zero modes at finite lattice spacing. On the other hand, as mentioned in section 3.1, DWF at $`L_s=\mathrm{}`$ have exact zero modes and at finite $`L_s`$ have robust zero modes to a good approximation . An approximate form of the index theorem has been found to be obeyed for fundamental fermions in the overlap formulation in quenched SU(2)  and in quenched SU(3) .

The index of adjoint fermions in the overlap formulation in quenched SU(2) has been studied in . In that work it was pointed out that the overlap Dirac operator for adjoint fermions in the SU(2) gauge group is necessarily even–valued. Then the question posed by the authors of  was whether or not all even values are realized or only values that are multiples of four are present. The latter case corresponds to configurations with instantons. The former case corresponds to fractional winding numbers. Configurations with fractional winding were found and their presence persisted as the lattice spacing was decreased

In this paper DWF are used at finite $`L_s`$ and therefore some of the clarity present in the $`L_s=\mathrm{}`$ case will be lost. However, the full theory (including the fermion determinant) is studied here. Furthermore, it is interesting to see if at a small volume and zero mass the gluino condensate still forms and if it does to what extent its value is due to zero mode effects. The numerical results are presented in section 5.

### 3.3 Ward identities

As discussed in the introduction and in section 3, the DWF formulation of the $`𝒩=1`$ SU(2) SYM theory at the $`L_s\mathrm{}`$ limit is expected to preserve the $`U(1)_A`$ chiral symmetry (at the classical level) and break supersymmetry only by irrelevant operators. Since the DWF formulation contains many more fields than the continuum theory, one may naturally wonder what are the SUSY transformations in terms of these fields. In particular, while the continuum theory has a single Majorana fermion the DWF lattice theory contains $`L_s`$ Majorana fermions and $`L_s`$ corresponding PV fields. Since all these fields, except for one Majorana fermion, have masses near the cutoff, one can expect that the SUSY transformations should only transform the gauge field and the light Majorana fermion represented by the boundary field $`\chi `$ of eq. 19. Similarly, the chiral symmetry transformations should only involve the field $`\chi `$. However, one should expect that this choice of SUSY and chiral transformations is not unique. For example, see  for a different choice of QCD chiral transformations that involve all fermion fields in one half of the fifth direction transforming vectorially and all fermions in the other half also transforming vectorially but with opposite charge. That choice could also be appropriate here for the chiral transformations, but it may make the SUSY ones more complicated.

As a first step in deriving the Ward identities, the fermionic part of the action in eq. 4 is rewritten in terms of the boundary field $`\chi `$:

$$S_F=S_{F_0}+S_{F_\chi }$$
(25)
where $`S_{F_0}`$ does not depend on the field $`\chi `$ and

$`S_{F_\chi }={\displaystyle \underset{x,x^{}}{}}`$ $`[\overline{\chi }(x)D/_N(x,x^{})\chi (x^{})m_f\overline{\chi }(x)\delta (x,x^{})\chi (x^{})`$ (26)
$`\overline{\chi }(x)B(x,x^{})\varphi (x^{})\overline{\varphi }(x)B(x,x^{})\chi (x^{})]`$
where

$`\varphi (x)`$ $`=`$ $`P_R\mathrm{\Psi }(x,L_s1)+P_L\mathrm{\Psi }(x,0)`$
$`\overline{\varphi }(x)`$ $`=`$ $`\overline{\mathrm{\Psi }}(x,0)P_R+\overline{\mathrm{\Psi }}(x,L_s1)P_L`$ (27)
are the “wrong” projected fields in the sense that they are defined on the opposite wall from where the corresponding light mode is localized. If indeed there is localization one would expect that in the $`L_s\mathrm{}`$ limit these fields will have no overlap with the light mode. The operator $`D/_N`$ is the naive part of the four-dimensional Wilson operator in eq. 6 and $`B`$ is the symmetry breaking part ($`B`$ is the equivalent of $`B`$ in ):

$$D/_N(x,x^{})=\frac{1}{2}\underset{\mu =1}{\overset{4}{}}\gamma _\mu \left[V_\mu (x)\delta (x+\widehat{\mu }x^{})V_\mu ^{}(x^{})\delta (x^{}+\widehat{\mu }x)\right]$$
(28)
$$B(x,x^{})=(5m_0)\delta (xx)\frac{1}{2}\underset{\mu =1}{\overset{4}{}}[V_\mu (x)\delta (x+\widehat{\mu }x^{})+V_\mu ^{}(x^{})\delta (x^{}+\widehat{\mu }x)].$$
(29)
These operators have the following properties:

$`\{D/_N,\gamma _5\}`$ $`=`$ $`0,D/_N^{}=D/_N`$
$`[B,\gamma _5]`$ $`=`$ $`0,B^{}=B^T=B`$ (30)
First the Ward identity corresponding to the $`U(1)_A`$ symmetry is derived. The symmetry transformations are:

$`\delta _A\chi (x)`$ $`=`$ $`i\alpha (x)\gamma _5\chi (x)`$
$`\delta _A\overline{\chi }(x)`$ $`=`$ $`i\alpha (x)\overline{\chi }(x)\gamma _5`$ (31)
where $`\alpha (x)`$ is an infinitesimal real number and $`\delta _A`$ symbolizes the change under the chiral transformation. Then the Ward identity is:

$$\mathrm{\Delta }_\mu J_\mu (x)𝒪(y)=2m_fJ_5(x)𝒪(y)+2J_B(x)𝒪(y)+i\delta _A𝒪(y)$$
(32)
where the backward difference is defined as $`\mathrm{\Delta }_\mu f(x)f(x)f(x\mu )`$. The currents are:

$$J_\mu (x)=\overline{\chi }(x)\gamma _5\gamma _\mu V_\mu (x)\chi (x+\mu )$$
(33)
$$J_5(x)=\overline{\chi }(x)\gamma _5\chi (x)$$
(34)
$$J_B(x)=\underset{y}{}\overline{\chi }(x)\gamma _5B(x,y)\varphi (y).$$
(35)
If in the above Ward identity $`𝒪(y)=J_5(y)`$ one gets

$$\mathrm{\Delta }_\mu J_\mu (x)J_5(y)=2m_fJ_5(x)J_5(y)+2J_B(x)J_5(y)2\overline{\chi }(y)\chi (y)$$
(36)
In this identity the term with $`J_B`$ will be responsible for producing the ABJ anomaly in the $`L_s\mathrm{}`$ limit. On the other hand, if $`L_s`$ is kept finite this term is similar to the one for Wilson fermions which, besides producing the ABJ anomaly, also produces a mass redefinition. For an analysis of QCD with DWF at finite $`L_s`$ see .

As mentioned earlier these chiral transformations are different than the ones in . If the transformations relevant for a non-singlet current in QCD were done on the fields $`\overline{\chi }`$, $`\chi `$, one obtains a Ward identity exactly as in  but with the currents $`A_\mu ^a(x)`$ and $`J_{5q}^a(x)`$ replaced with:

$$A_\mu ^a(x)=\frac{1}{2}\left[\overline{\chi }(x)\gamma _5\gamma _\mu \lambda ^aU_\mu (x)\chi (x+\mu )+\overline{\chi }(x+\mu )\gamma _5\gamma _\mu \lambda ^aU_\mu ^{}(x)\chi (x)\right]$$
(37)
$$J_{5q}^a(x)=\frac{1}{2}\underset{y}{}\left[\overline{\chi }(x)\gamma _5\lambda ^aB(x,y)\varphi (y)+\overline{\varphi }(y)\gamma _5\lambda ^aB(y,x)\chi (x)\right].$$
(38)
The derivation of the SUSY Ward identity is similar to the one for Wilson fermions. One can use the existing calculations for Wilson fermions  to elucidate the differences between the two formalisms. Here the derivation in  will be followed. The symmetry transformations are as in  and commute with parity.

The change of the pure gauge action with respect to the transformation of the gauge field is of course the same. In terms of the symmetry breaking part of the Ward identity it contributes a term denoted below by $`X_2(x)+X_3(x)`$ where $`X_2`$, $`X_3`$ are as in . This term breaks SUSY because of the explicit breaking of the Lorentz symmetry. Using improved pure gauge lattice actions can alleviate the effects of this breaking. Such an improvement is not considered here.

The change of the fermion and Pauli-Villars parts of the DWF action with respect to the transformation of the gauge fields produces terms for all $`L_s`$ slices. In particular the variation of the fermion matrix $`D_F`$ of eq. 5 with respect to the gauge field is:

$$\delta _UD_F(x,x^{};s,s)(m_f)=\delta (ss^{})\delta _UD/(x,x^{})$$
(39)
One sees that $`\delta D_F`$ is independent of $`m_f`$ and is diagonal in the fifth direction. Furthermore $`\delta D/(x,x^{})`$ is the same as the variation of the Wilson operator. Therefore, this variation contributes to the symmetry breaking part of the Ward identity the terms:

$$X_4^F(x)=\underset{s}{}X_4^F(x,s)$$
(40)
and

$$X_4^{PV}(x)=\underset{s}{}X_4^{PV}(x,s)$$
(41)
where $`X_4^F(x,s)`$ is as $`X_4`$ in  except that the four-dimensional Wilson fermion fields that have their spin indices contracted are replaced by $`\overline{\mathrm{\Psi }}(x,s)`$, $`\mathrm{\Psi }(x,s)`$ while the other Wilson fermion field is replaced by $`\chi (x)`$. Similarly $`X_4^{PV}(x,s)`$ is as $`X_4`$ in  except that the four-dimensional Wilson fermion fields that have their spin indices contracted are replaced by the Pauli-Villars fields $`\mathrm{\Phi }^T(x,s)CR_5`$, $`\mathrm{\Phi }(x,s)`$, the other Wilson fermion field is replaced by $`\chi (x)`$ and the sign of the second term in $`X_4`$ is minus instead of plus due to the commutativity of the Pauli-Villars fields.

The change of the action with respect to the fermion field transformations can be partially deduced from the corresponding Wilson fermion calculation. Since this transformation only involves the action $`S_{F_\chi }`$ in eq. 26, one can observe that the first two terms of that action are identical with the action of naive fermions (Wilson fermions with $`r=0`$). These will contribute identical terms as the $`r=0`$ part of the Wilson action. They contribute to the divergence of the SUSY current and to the mass term of the Ward identity given below. Finally, the transformation of the last term of the action $`S_{F_\chi }`$ in eq. 26 is easy to calculate and is denoted by $`X_1(x)`$:

$$X_1(x)=P_{\rho \sigma }^a(x)\sigma _{\rho \sigma }B^{ab}(x,x^{})\varphi ^b(x^{}).$$
(42)
This term is closely related to $`X_1`$ of .

The Ward identity is:

$$\mathrm{\Delta }_\mu S_\mu (x)𝒪(y)=m_fD_s(x)𝒪(y)+X_S(x)𝒪(y)\delta _S𝒪(y)$$
(43)
where the supersymmetric current $`S_\mu `$ and the quantity $`D_S`$ are as in . The symmetry breaking term $`X_S(x)`$ is also similar to the one in :

$$X_S(x)=X_1(x)+X_2(x)+X_3(x)+X_4^F(x)X_4^{PV}(x)$$
(44)
As mentioned above the symmetry breaking term $`X_2(x)+X_3(x)`$ is due to the breaking of Lorentz symmetry by the lattice. The $`X_4^F(x)`$ and $`X_4^{PV}(x)`$ terms break the symmetry as in Wilson fermions. These terms do not cancel each other exactly <sup>1</sup><sup>1</sup>1We thank Y. Shamir for pointing this out to us. However, one would expect large cancellations of heavy modes. The terms in $`X_4^F(x,s)`$ that are proportional to the Wilson parameter involve fields that couple to the light modes by an amount that is exponentially small in $`L_s`$. One would expect these terms to be nearly canceled by the corresponding Pauli-Villars terms resulting in exponentially small contributions. The remaining terms that involve fields away from the relevant domain walls should also yield similar cancellations. As a result the only terms that should make significant contributions should be the ones that involve fields of the “correct” chirality near the domain walls. These few terms would couple to the light modes and be further regularized by the corresponding Pauli-Villars terms. Clearly this analysis of cancellations is heuristic. A detailed calculation using for example perturbation theory or transfer matrix methods would be interesting but it is beyond the scope of this paper.

Finally, the symmetry breaking term $`X_1(x)`$ involves the field $`\varphi (x)`$ that is expected to have no overlap with the light mode in the $`L_s\mathrm{}`$ limit. If $`L_s`$ is finite then DWF are similar to Wilson fermions and an analysis as in  should indicate that this term is responsible for the same mass redefinition as the one in the chiral Ward identity.

## 4 Numerical methods

As can be seen from section 2 the $`𝒩=1`$ SU(2) SYM theory can be simulated as a theory with 0.5 flavors of Dirac fermions in the adjoint representation. An efficient and popular algorithm that can be used to simulate any number of flavors is the hybrid molecular dynamics R (HMDR) algorithm of . Because of the Grassmann nature of fermions these algorithms need to invert the matrix $`D_F`$ of eq. 5. That matrix is not Hermitian. This is a problem since some of the more efficient inversion algorithms require the matrix to be Hermitian. However, because:

$$\gamma _5RD_FR\gamma _5=D_F^{}$$
(45)
one has that:

$$det[D_F]^2=det[D_FD_F^{}].$$
(46)
Then one can invert the Hermitian matrix $`D_FD_F^{}`$ and then use the HMDR algorithm to take the appropriate power so that the desired number of flavors is simulated. This method is adopted here and the 0.25 power is taken in order to go from a theory with two Dirac fermions to a theory with one Majorana fermion. In other words, the fermion determinant that is used in the simulation is:

$$(det[D_FD_F^{}])^{0.25}=|det[D_F]|^{0.5}=det[D_F]^{0.5}$$
(47)
where in the last equality use was made of the fact that for non-negative $`m_f`$ $`det[D_F]`$ is also non-negative . This approach was also taken in  for Wilson fermions. For an approach that uses Wilson fermions and the multibosonic algorithm  instead of the HMDR algorithm see . However, as mentioned earlier in the case of Wilson fermions the last equality in eq. 47 is not true for all gauge field configurations.

The HMDR algorithm uses molecular dynamics methods in order to produce the correct statistical ensembles. Because the molecular dynamics step size $`\delta \tau `$ is finite discretization errors are introduced. There are two ways one can deal with this problem. One is to simulate at various values of $`\delta \tau `$ and then extrapolate to $`\delta \tau =0`$. Another method is to use $`\delta \tau `$ small enough so that the errors are negligible when compared with the statistical errors.

In order to ensure this, one can simulate the two Dirac flavor theory at the same parameters and same $`\delta \tau `$. For the two flavor theory, one has a local action and therefore, at the end of the evolution, one can employ a Metropolis accept-reject step. Then the finite $`\delta \tau `$ errors are “converted” to a non-ideal acceptance rate and in effect they are reflected in the final statistical errors. This is the exact hybrid Monte Carlo $`\mathrm{\Phi }`$ (HMC$`\mathrm{\Phi }`$) algorithm of . Therefore the acceptance rate is an indication of the size of the finite $`\delta \tau `$ errors in the HMD integration. By simulating the two Dirac flavor theory with (HMC$`\mathrm{\Phi }`$) one can set $`\delta \tau `$ so that the acceptance rate is high, say $`90\%`$. Since the coefficient of the finite $`\delta \tau `$ errors is proportional to the number of flavors one would expect that for 0.25 flavors the errors would be small and at the few percent level.

The only fermion observable measured in this work is the gluino condensate. By inserting appropriate source terms as in  the gluino condensate was measured as the trace of $`D_F^1`$ with spin and fifth-direction indices restricted as dictated by eq. 19 The trace was calculated using a standard stochastic method. All inversions in this work were done using the conjugate gradient (CG) algorithm. An even-odd preconditioned form of the matrix $`D_F^{}D_F`$ was inverted. For more details on the numerical algorithms and methods employed to DWF simulations see .

## 5 Numerical results

### 5.1 Simulation parameters

In all simulations the domain wall height was chosen to be $`m_0=1.9`$. As mentioned in the previous section, the finite $`\delta \tau `$ errors were kept to the few percent level by using a small $`\delta \tau `$. For all simulations the step size was set to $`\delta \tau =0.01`$ and the trajectory length to $`\tau =0.5`$. In order to confirm that this choice introduces finite step size errors that are small compared to the statistical errors an HMC$`\mathrm{\Phi }`$ simulation for two Dirac flavors was run for $`L_s=12`$ and $`m_f=0.04`$. It produced an acceptance rate of $`90\%`$ suggesting that the finite $`\delta \tau `$ errors of the $`0.5`$ flavor theory are small. Furthermore, an HMDR simulation was also run for two Dirac flavors using the exact same parameters. The value of the gluino condensate obtained from these two simulations was the same within statistical errors.

The CG stopping condition for all simulations was set to $`10^6`$ for the evolution and to $`10^8`$ for the calculation of $`\overline{\chi }\chi `$. The number of CG iterations varied between $`100`$ for $`m_f=0.08`$, $`L_s=12`$ and $`250`$ for $`m_f=0.0`$, $`L_s=24`$.

The $`8^4`$ volume simulations were done with $`\beta =2.3`$. The value of $`\beta `$ was chosen so that one is not close to the point where the box size becomes too small and a thermal transition takes place, but also not too deep in the strong coupling regime where the finite $`L_s`$ effects become severe. The transition point of the $`N_t=8`$ quenched theory is at $`\beta =2.5115(40)`$ . In figure 1 the magnitude of the fundamental Wilson line $`|W|`$ measured in quenched simulations in an $`8^4`$ volume is plotted vs. $`\beta `$. In the quenched theory this is an order parameter. As can be seen from that figure, a rapid crossover takes place around $`\beta =2.5`$. In the same figure the value of $`|W|`$ from a simulation of the dynamical theory at $`\beta =2.3`$ is also shown (cross). The quenched and dynamical values are very similar indicating that at $`\beta =2.3`$ the dynamical theory is in a phase that “confines” fundamental sources. Therefore, the box size is large enough to avoid finite temperature effects that would of course spoil SUSY. Using the quenched theory as a guide the $`4^4`$ simulations were done at $`\beta =2.1`$ since the quenched transition at $`N_t=4`$ is known to take place at $`\beta =2.2986(6)`$ . At $`\beta =2.1`$ the lattice spacing is larger than at $`\beta =2.3`$. However, the lattice sizes are small and do not allow a reliable measurement of the lattice spacing. According to , $`\beta =2.12.3`$ is in the beginning of the weak coupling regime. Then if one uses the weak coupling form in  one finds that the lattice spacing at $`\beta =2.1`$ is about a factor of two larger than the one at $`\beta =2.3`$.

In order to estimate the necessary number of thermalization sweeps two simulations were run on an $`8^4`$ lattice at $`\beta =2.3`$, $`L_s=12`$ and $`m_f=0.04`$. The first simulation used an initial configuration with all gauge links set to the identity (ordered) and the other used an initially random configuration (dis-ordered). The evolutions in “computer time” are shown in figure 2. As can be seen, the two ensembles converged after about $`100`$ sweeps. This number of thermalization sweeps was then used in all other simulations which were started from an ordered initial configuration. The number of measurements after thermalization for all simulations is about $`200`$ with measurements done in every trajectory. The gluino condensate was measured with a single “hit” stochastic estimator.

### 5.2 The gluino condensate at the chiral limit

In order to be able to extrapolate to the chiral limit, corresponding to $`L_s\mathrm{}`$ and $`m_f=0`$, the mass $`m_f`$ and the size of the fifth direction $`L_s`$ was varied. The results of all simulations are given in tables 1 and 2 in the appendix. Three different methods were used to analyze the data and calculate the gluino condensate in the chiral limit.

I. For fixed $`L_s`$, the data for $`m_f=0.08,0.06,0.04,0.02`$ were fit to a function:

$$b_0+b_1m_f.$$
(48)
This functional form is valid provided $`m_f`$ is small enough. Otherwise, higher order terms must also be included. The data and fits are shown in figure 3 and the results of the fits are given in table 3. Then the extrapolated values $`b_0`$ were fit to a form:

$$c_0+c_1\mathrm{exp}(c_2L_s).$$
(49)
This functional form is approximate but it is expected to be valid close enough to the continuum and has been found to be consistent in simulations of the Schwinger model  and of QCD even at relatively large lattice spacings (see for example ). The data and fit is shown in figure 4 and the results of the fit are given in table 4.

II. For fixed $`m_f`$ the data for $`L_s=12,16,20,24`$ were fit to the form of eq. 49. The data and fits are shown in figure 5 and the results of the fit are given in table 4. Then the extrapolated values $`c_0`$ were fit to the form of eq. 48. The data and fit is shown in figure 6 and the results of the fit are in table 3.

III. Additional simulations were done for $`m_f=0`$ and $`L_s=12,16,20,24`$. The data were fit to the form of eq. 49. The data and fits are shown in figure 7 and the results of the fit are in table 4.

The $`m_f0`$ and $`L_s\mathrm{}`$ extrapolated values of the gluino condensate for each one of the above three methods are summarized in table 5. As can be seen, all values are consistent within the statistical errors. This suggests that the systematic errors inherent to the limited statistics and to fits onto functions that represent the data only for a limited range are small. Furthermore, it suggests that the fitting functions used are consistent (please see subsection 5.4 for more discussion on the validity of these fitting functions).

### 5.3 The telltale signals of topology in numerical simulations

In order to investigate the issues discussed in section 3.2 the gluino condensate was also calculated in a smaller $`4^4`$ lattice volume at $`\beta =2.1`$. It was measured only for $`m_f=0`$ and method III above was used to extrapolate to the $`L_s\mathrm{}`$ limit. The data and fit are shown in figure 8 and the fit results are given in table 4. The $`8^4`$ data from figure 7 are presented again in this figure to aid comparison. The value has decreased indicating that scaling is violated. However, without more simulations at other lattice spacings and volumes one can not conclude much from this result. The $`\beta =2.1`$ coupling is in the strong coupling region and furthermore the $`4^4`$ lattice volume is rather small.

However, it is interesting to notice that the parameter $`V\times \overline{\chi }\chi _{L_s\mathrm{}}8.4`$ (a factor of 12 coming from the normalization of $`\overline{\chi }\chi `$ has been included). Since $`m_f=0`$ the effective mass $`m_{\mathrm{eff}}`$ gets its value from finite $`L_s`$ effects. As $`L_s`$ is increased $`m_{\mathrm{eff}}`$ becomes small. From analysis of $`m_{\mathrm{eff}}`$ in strong coupling QCD  one would roughly guess that $`m_{\mathrm{eff}}<0.1`$. Then $`[m_{\mathrm{eff}}\times V\times \overline{\chi }\chi _{L_s\mathrm{}}]<1`$. In that case, the analysis of  can be followed and one would expect the value of the condensate in the $`4^4`$ lattice to be mostly supported by configurations with total winding of $`\pm 1/2`$. Indeed, this can be seen from figure 9. In that figure the evolutions in “computer time” are shown. The “spikes” in the evolution are apparent and they become more pronounced and less frequent as $`L_s`$ is increased (and in effect $`m_{\mathrm{eff}}`$ is decreased). This is exactly how the effect of zero modes for small $`[m_{\mathrm{eff}}\times V\times \overline{\chi }\chi _{L_s\mathrm{}}]`$ would present itself in a numerical simulation of the dynamical theory. As the fermion mass is made smaller, $`\overline{\chi }\chi `$ is expected to receive most of its value from sectors with winding $`\pm 1/2`$. However, in these sectors the fermion determinant is very small because of the zero mode. Since the probability for the algorithm to generate a gauge field configuration is proportional to the fermion determinant one would expect that these sectors will be visited less and less frequently as the effective mass is decreased. When these sectors are visited the value of $`\overline{\chi }\chi `$ will be very large (spikes) in order to compensate for the infrequent sampling. In this way the presence of the zero mode in the observable “balances” the presence of the zero mode in the determinant. As the mass is made smaller one would have to increase the size of the statistical sample in order to include enough of these increasingly “rare” but very large fluctuations. For similar results in the Schwinger model and QCD see .

A histogram of the values of $`\overline{\chi }\chi `$ is presented in figure 10 (solid line). For small $`L_s`$ the effective mass is larger and $`\overline{\chi }\chi `$ is distributed with a symmetric looking distribution around the mean value. However, for $`L_s=40`$ the effective mass is smaller and the distribution has a more pronounced “tail” towards larger values. In order to investigate this further numerical simulations at exactly the same parameters, but without the fermion determinant (quenched theory) were done. The histograms from these simulations are shown in the same figure for comparison (dotted lines). One can observe that the absence of the fermion determinant had the effect of shifting the distributions to higher values. This is expected since configurations with small eigenvalues that produce larger values of $`\overline{\chi }\chi `$ are not suppressed anymore and are produced more frequently. Also, one can observe that the number of configurations with $`\overline{\chi }\chi `$ larger than $`0.007`$ that appeared as spikes in figure 9 have now increased in number. These observations lend support to the presence of small near-zero eigenvalues. Furthermore, configurations with fractional topological charge have already been found in quenched SU(2) simulations at similar couplings . It would be very interesting to calculate the index for the configurations of figure 9 using the methods of  and see to what extent there is a correlation between fractional topological charge and the observed spikes. This correlation should be exact at $`L_s\mathrm{}`$ but it will be obscured at finite $`L_s`$ by the presence of non-zero $`m_{\mathrm{eff}}`$. This investigation is beyond the scope of this work.

Furthermore, it should also be noted that on the $`8^4`$ lattice there are no visible spikes up to $`L_s=24`$. This can be seen in figure 11. Presumably this is because the product $`[m_{\mathrm{eff}}\times V\times \overline{\chi }\chi _{L_s\mathrm{}}]`$ is probably much larger than in the $`4^4`$ lattice. Again, this statement is not exact since the value of $`m_{\mathrm{eff}}`$ was not measured.

These results are consistent with the discussion in section 3.2. However, since even with $`m_f=0`$, an $`L_s`$ extrapolation is essentially an extrapolation from non zero masses these results are not necessarily the results of a simulation at exactly $`L_s=\mathrm{}`$. It is still possible that if such a simulation were done one could have found that the gluino condensate is zero. This could happen since in a finite volume and zero mass the effects of spontaneous symmetry breaking are absent and the zero mode effects alluded to above may not be sufficient to sustain a non zero vacuum expectation value. This type of simulation is possible and can be done using the overlap formalism  or exact Neuberger fermions . However, if one is to maintain exact chiral symmetry these methods will demand large computing resources.

### 5.4 The fine print

Perhaps the largest uncertainty in the analysis presented in the previous subsections has to do with the assumption of exponential decay as in eq. 49. For small enough lattice spacings and large enough $`L_s`$ this behavior is expected to be true. All data presented in this work were well represented by this ansatz. However, as with any numerical investigation, one can never completely disprove all other possibilities. While such an exercise over all possible functions would clearly be fruitless there are some alternative forms that may be reasonable to consider since they are based on analytical considerations.

Far from the continuum limit, the approach to the chiral limit may become power law  or even completely disappear . In order to explore the possibility of power law behavior the $`m_f=0`$ data for the $`8^4`$ and $`4^4`$ volumes were fit to the form:

$$d_0+d_1L_s^{d_2}.$$
(50)
The results of the fit are given in table 6 (the fits are not presented in any of the figures). As can be seen from that table the $`\chi ^2/\mathrm{dof}`$ of these fits is significantly larger than the one of the corresponding exponential fits to the same data.

Another possibility is decay to zero with two different exponential decay rates. Such a behavior was found to be consistent with investigations of the two flavor Schwinger model  for a quantity that is expected to vanish in the chiral limit. There it was argued that the fast decay rare is due to fluctuations within a given topological sector while the slow decay rate is due to the presence of topology changing configurations. Therefore, for $`m_f=0`$ one could try to fit the largest three $`L_s`$ points to the form:

$$e_0\mathrm{exp}(e_1L_s).$$
(51)
The results of the fit for the $`8^4`$, $`\beta =2.3`$, $`m_f=0`$, and $`L_s=16`$, $`20`$, $`24`$ points as well as for the $`4^4`$, $`\beta =2.1`$, $`m_f=0`$, and $`L_s=24`$, $`32`$, $`40`$ points are shown in table 7. The fit to the $`8^4`$ data is acceptable. However, the fit to the $`4^4`$ data has a rather large $`\chi ^2/\mathrm{dof}`$. Because this fit is for larger $`L_s`$ than the $`8^4`$ fit one would expect that if there were a second exponential decaying to zero its effect would be more pronounced in the $`4^4`$ fit. Therefore, the large $`\chi ^2/\mathrm{dof}`$ of the $`4^4`$ fit suggests that the presence of a second exponential decaying to zero is not likely. This could be made more precise if simulations with larger $`L_s`$ values for the $`8^4`$ and $`4^4`$ lattices were done. However, such simulations are beyond the computing resources of this project. Also, the analysis in  suggests functional forms with more parameters. It would be interesting to fit to these forms but that would require more data points and better statistics both of which are also beyond the computing resources of this project.

Finally, the SUSY breaking by the irrelevant terms may have non-negligible effects at the lattice spacings studied here. Although it was found that the chiral condensate is non–zero at the chiral limit in two lattice spacings, this is not enough to estimate its value in the continuum limit.

## 6 Conclusions

The formulation of $`𝒩=1`$, SU(2) Super Yang-Mills theory on the lattice with domain wall fermions (DWF) has several advantages over more traditional lattice fermion regulators. Even at non-zero lattice spacing the chiral limit can be taken by letting $`L_s\mathrm{}`$, where $`L_s`$ is the number of sites along the fifth auxiliary direction. Since in that limit there is no gluino mass term, supersymmetry is broken only by irrelevant operators and there is no need for fine tuning. Also, in that limit the theory has exact zero modes on non-trivial topological backgrounds.

However, even at finite $`L_s`$, where numerical simulations are done, these properties are maintained to a good degree allowing extrapolations to the $`L_s\mathrm{}`$ limit. Furthermore, the Pfaffian resulting from the integration of Majorana fermions is positive definite at finite $`L_s`$, non-zero lattice spacing and for any background gauge field configuration. As a result, one can unambiguously interpret it as a probability measure to be used by the numerical simulation for importance sampling. This property also allows the use of standard numerical algorithms where any number of flavors $`N_f`$ can be simulated. By contrast, Wilson fermions have this positivity property only at the continuum limit.

In this work, the full $`𝒩=1`$, SU(2) Super Yang-Mills theory was numerically simulated on the lattice using DWF. The gluino condensate $`\overline{\chi }\chi `$ was measured. These simulations did not present any unexpected technical difficulties.

A finite value of $`L_s`$ breaks chiral symmetry and induces a small gluino mass. In addition, an explicit gluino mass $`m_f`$ was used to provide extra control. Several $`m_f`$ and $`L_s`$ values were used (all corresponding to positive gluino mass) and the value of $`\overline{\chi }\chi `$ was extrapolated to the chiral limit using three different methods. All methods gave consistent results indicating small systematic effects and suggesting that the functions used for the fits are consistent. These simulations were done on a lattice with $`8^4`$ lattice sites.

Additional simulations on a lattice with $`4^4`$ lattice sites but approximately double the lattice spacing were done. Again, extrapolations to the chiral limit gave a non-zero $`\overline{\chi }\chi `$. In this lattice $`\left[\mathrm{mass}\times \mathrm{volume}\times \overline{\chi }\chi \right]<1`$. Then analytical considerations suggest that the value of $`\overline{\chi }\chi `$ must come mostly from topological sectors with fractional topological charge of $`\pm 1/2`$. Indeed, as the mass was made smaller unusually large values (spikes) were observed in the statistical sample of $`\overline{\chi }\chi `$ indicating the singular contribution of these sectors.

The spectrum of the theory is of great interest but it was not possible to measure on the small lattices considered here. Also, the gluino condensate was measured only on two different lattice spacings and therefore it was not possible to extrapolate to the continuum limit where comparisons with analytical results would be possible. Future work could explore these very interesting topics.

## Acknowledgments

This research was supported in part by NSF under grant # NSF-PHY96-05199 (J. Kogut and P. Vranas). All numerical simulations were done on an 18 GFlops part of the QCDSP supercomputer at Columbia University and on the 6 GFlops QCDSP supercomputer at Ohio State University. We would like to thank N. Christ and R. Mawhinney for providing us with the Columbia University resource and G. Kilcup for providing us with the Ohio State University resource. Also, we would like to thank Y. Shamir and T. Schaefer for useful comments. P. Vranas would like to thank D.B. Kaplan, A. Kovner, A. Nyffeler, and E. Weinberg for enlightening discussions.

## 7 Appendix

In this appendix the various tables are presented.


format: text

text:
# Four neutrino oscillation analysis of atmospheric neutrino data and application to long baseline experiments

## 1 Four neutrino analysis of atmospheric neutrinos

Four neutrino mixing schemes have caught much interest, since they are the simplest scenario which accounts for the solar and atmospheric neutrino problems and the LSND data in the framework of neutrino oscillations. To reconcile the data at 90%CL of LSND, Bugey and CDHSW one has to have two degenerate massive states ($`m_1^2m_2^2m_3^2m_4^2`$). For simplicity I assume $`\mathrm{\Delta }m_{21}^2=\mathrm{\Delta }m_{}^2`$, $`\mathrm{\Delta }m_{43}^2=\mathrm{\Delta }m_{\text{atm}}^2`$ and I adopt the notation in  for the $`4\times 4`$ MNS matrix:

$`\left(\begin{array}{c}\nu _e\\ \nu _\mu \\ \nu _\tau \\ \nu _s\end{array}\right)=\left(\begin{array}{cccc}U_{e1}& U_{e2}& U_{e3}& U_{e4}\\ U_{\mu 1}& U_{\mu 2}& U_{\mu 3}& U_{\mu 4}\\ U_{\tau 1}& U_{\tau 2}& U_{\tau 3}& U_{\tau 4}\\ U_{s1}& U_{s2}& U_{s3}& U_{s4}\end{array}\right)\left(\begin{array}{c}\nu _1\\ \nu _2\\ \nu _3\\ \nu _4\end{array}\right).`$ (13)
For the range of the $`\mathrm{\Delta }m^2`$ suggested by the LSND data at 90%CL, which is given by 0.2 eV$`{}_{}{}^{2}<\mathrm{\Delta }m_{\text{LSND}}^2<`$ 2 eV2 when combined with the data of Bugey E776 and KARMEN2, the constraint by the Bugey data gives $`|U_{e3}|^2+|U_{e4}|^21`$. Also $`|U_{s3}|^2+|U_{s4}|^2`$ has to be very small if one demands that the number $`N_\nu `$ of effective neutrinos in Big Bang Nucleosynthesis (BBN) be less than four. In this case the solar neutrino deficit is explained by $`\nu _e\nu _s`$ oscillations with the Small Mixing Angle (SMA) MSW solution and the atmospheric neutrino anomaly is accounted for by $`\nu _\mu \nu _\tau `$. However, some people have given conservative estimate for $`N_\nu `$ and if their estimate is correct then $`|U_{s3}|^2+|U_{s4}|^21`$ may no longer hold. Recently the Superkamiokande group has reported their result on the solar neutrino experiment which indicates that the SMA and the Vacuum Oscillation solutions are disfavored. Meanwhile Giunti, Gonzalez-Garcia and Peña-Garay have analyzed the solar neutrino data in the four neutrino scheme without BBN constraints. Their updated results show that the SMA solution exists for $`0c_s<0.8`$ ($`c_s|U_{s1}|^2+|U_{s2}|^2`$), while the Large Mixing Angle (LMA) and LOW solutions survive only for $`0c_s<0.4`$ and $`0c_s<0.2`$, respectively. In this talk I will present some of the updated results of my work on the four neutrino oscillation analysis of the Superkamiokande atmospheric neutrino data and I will briefly give some implications to long baseline experiments, assuming that the solar neutrino deficit is solved by the LMA solution. For details of the analysis and the references see . In the analysis of atmospheric neutrinos, the effect of $`\mathrm{\Delta }m_{}^2`$ is negligible, so I assume $`\mathrm{\Delta }m_{21}^2=0`$ and $`U_{e3}=U_{e4}=0`$ for simplicity. To avoid contradiction with the CDHSW data, I will take $`\mathrm{\Delta }m_{32}^2`$=0.3eV2 as a reference value. The best fit to the atmospheric neutrino data for 1144 days is obtained for $`\mathrm{\Delta }m_{43}^2=2.0\times 10^3\text{eV}^2`$, $`(\theta _{24},\theta _{34},\theta _{23})=(45^{},30^{},20^{})`$, $`\delta _1=45^{}`$ and the allowed region at 90%CL is obtained (The allowed region for 1144 day data is almost the same as for 990 day data). As a sample let me show the result for the case of $`\theta _{24}=45^{},\delta _1=90^{}`$. The shadowed area in Fig.1 (a) is the 90%CL allowed region projected on the $`(\theta _{34},\theta _{23})`$ plane for various values of $`\mathrm{\Delta }m_{43}^2`$. If I demand that the solar neutrino deficit be solved by the LMA solution, then $`c_s<0.4`$ which is depicted in Fig.1 (b). Combining the results on the atmospheric neutrinos and the solar neutrinos, the allowed region becomes the shadowed area in Fig.1 (c). It turns out that if I require $`c_s<0.4`$ then the solution prefers relatively large $`\theta _{23}`$ for any value of $`\delta _1`$.

## 2 Application to long baseline experiments

First implication of the present scheme is the deficit which is expected to be seen at the K2K experiments. Ignoring the matter effect, I have for the neutrino energy $`E`$ 1GeV and the path length $`L`$=250km

$`P(\nu _\mu \nu _\mu )=12|U_{\mu 2}|^2(1|U_{\mu 2}|^2)`$
$`4|U_{\mu 3}|^2|U_{\mu 4}|^2\mathrm{sin}^2\left({\displaystyle \frac{\mathrm{\Delta }m_{\text{atm}}^2L}{4E}}\right),`$ (14)
where $`|U_{\mu 2}|=|c_{24}s_{23}|`$ is not necessarily small. The behavior of the ratio $`r`$#(CC + NC events with oscillations)/ #(CC + NC events without oscillations) as a function of $`\mathrm{\Delta }m_{\text{atm}}^2`$ is given in Fig.2. Because of the contribution of the first term on the RHS of (14) the ratio is lower for the present four neutrino scheme than for the standard two flavor case with the same $`\mathrm{\Delta }m_{\text{atm}}^2`$.

Second implication is the possible measurement of CP violation at the JHF experiment which is expected to start from 2006. Since there is no strong constraint on the mixing angles $`\theta _{24}`$, $`\theta _{34}`$, $`\theta _{23}`$, CP violation in the channel $`\nu _\mu \nu _s`$ could be large. If we measure the neutral current $`\pi ^0`$ production for $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ beams and compare them, then the absolute value of the ratio

$`R{\displaystyle \frac{\frac{N_\nu (\delta _1)}{N_{\overline{\nu }}(\delta _1)}|_{\text{dat}}\frac{N_\nu (\delta _1=0)}{N_{\overline{\nu }}(\delta _1=0)}|_{\text{MC}}}{\frac{N_\nu (\delta _1)}{N_{\overline{\nu }}(\delta _1)}|_{\text{dat}}+\frac{N_\nu (\delta _1=0)}{N_{\overline{\nu }}(\delta _1=0)}|_{\text{MC}}}}`$
could be significantly larger than the statistical fluctuation ($``$ 0.05 for $`10^{21}`$ POT) (The yields is given in Table 1 for an optimistic set of the oscillation parameters ($`\delta _1=90^{}`$, $`\theta _{24}=40^{}`$, $`\theta _{23}=30^{}`$, $`\theta _{34}=60^{}`$, $`\mathrm{\Delta }m_{\text{atm}}^2=1.6\times 10^3`$eV2). The ratio $`R`$ is depicted in Fig.3 where the same set of the oscillation parameters as in Table 1 and the Wide Band Beam is assumed in calculations).

## 3 Conclusions

I have presented my result on the four neutrino oscillations of the Superkamiokande atmospheric neutrino data without assuming the BBN constraints. By combining the analysis on the solar neutrino data by Giunti et al. and assuming the LMA solar solution, I found that there is relatively large contribution of $`\mathrm{\Delta }m_{\text{LSND}}^2`$ in the atmospheric neutrino oscillations. It suggests that the number of events at K2K is less than the one for the standard two flavor scenario for the same $`\mathrm{\Delta }m_{\text{atm}}^2`$ and that CP violation could be measured at JHF after running for several years with $`10^{21}`$POT/yr.

## Acknowledgments

I would like to thank T. Nakaya for pointing out the NC $`\pi ^0`$ channel to me to measure CP violation. This research was supported in part by a Grant-in-Aid for Scientific Research of the Ministry of Education, Science and Culture, #12047222, #10640280.


text:
# Self-similar spherical collapse with non-radial motions

## 1 introduction

In the hierarchical scenario for structure formation, the initial density field is Gaussian with a fluctuation amplitude that decreases with scale (Peebles 1980). Non-linear gravitational evolution then causes matter to aggregate into bound virialized objects (halos) which are believed to harbor galaxies, and galaxy groups and clusters. The mass distribution in halos can be inferred from a variety of observations. These include observations of rotation curves of spiral galaxies (e.g., Persic et. al. 1996), motions of satellite galaxies (Zaritsky & White 1994), lensing distortions of background galaxy images by the potential wells of massive halos (e.g., Bartelmann & Schneider 1999), and velocity dispersion and X-ray maps of galaxy clusters (e.g., Bahcall & Fan 1998, Sarazin 1986). Nevertheless, no single observational method probes the mass profile over the entire extent of the halo. Therefore a conclusive analysis of the observations must rely on an assumed model for the halo mass profile. This is important for example in estimating the total masses of galaxy clusters. Further, halo profiles in the inner regions depend on the type of dark matter (e.g., Moore 1994). A comparison between observed and model mass profiles in the context of a given cosmological model can provide important information on the nature of dark matter.

We lack a complete theory of nonlinear gravitating systems. So a detailed study of nonlinear collapse into bound objects must rely on N-body simulations. The simulations seem to be converging on the shape of the mass profile in the halo inner regions (Navarro, Frenk & White 1997, Moore et. al. 1998, Jing 1999, Klypin et. al. 2000). As of yet no satisfactory analytic explanation of the simulations results has been suggested (but see Weinberg 2000, Syer & White 1998). In a generic collapse the infalling matter is in irregular clumps undergoing dynamical friction and tidal stripping as they sink towards the center (e.g., Avila-Rees et. al. 1998, Nusser & Sheth 1999). Moreover, a collapsed object today might have gone through a merger with another object of comparable mass, a process that might have an effect on its current mass profile. Because of these factors, the focus of analytic studies has been the collapse of perturbations with special configurations. An important step forward in these studies came with the realization that the potential in the inner collapse regions varies very little over the dynamical time scale (Gunn 1977), admitting an adiabatic invariant for the motion of a particle. Using adiabatic invariance, Fillmore & Goldreich (1984, hereafter FG84, see also Bertschinger 1985) have shown that in spherical symmetry an initial density perturbation $`\mathrm{\Delta }M/MM^ϵ`$ develops an asymptotic mass profile $`Mr`$ near the origin for $`0<ϵ2/3`$ and $`Mr^{3/(1+3ϵ)}`$ for $`ϵ>2/3`$. FG84 restricted their analysis to particles collapsing on purely radial orbits. Spherical collapse with non-radial particle motions has also been considered in the context of adiabatic invariance. Gurevich & Zybin (1988a,b) developed a formalism based on the theory of adiabatic capture (e.g., Lifshitz & Pitaevskii 1981), to estimate the mass profile in a spherical collapse including non-radial motions. Although their formalism can be useful to study the collapse of various initial conditions, Gurevich & Zybin applied it only to the collapse near the center of a smoothed initial density maxima. Ryden & Gunn (1987) included non-radial motions in their application of adiabatic invariance in a numerical scheme to study the evolution of a spherical peak in a Gaussian density field (Bardeen et. al. 1986, Hoffman & Shaham 1985).

White & Zartizky (1991) described the $`Mr`$ profile for $`ϵ2/3`$ in the purely radial spherical collapse in terms of the crowding of orbits as particles pass through the center. They conjectured that if particles had non-radial orbits then the crowding effect is avoided and the scaling $`Mr^{3/(1+3ϵ)}`$ $`ϵ>2/3`$ would also be valid for $`0<ϵ2/3`$. Here we examine this conjecture in detail by generalizing the analysis of FG84. We consider two schemes for assigning angular momenta to the collapsing particles.

The outline of this paper is as follows. In section 2 we write the equations of motion and the relevant initial conditions and present two schemes for assigning angular momenta to particles. In section 3 we discuss self-similarity and adiabatic invariance in the context of non-radial motions. In section 4 we derive the asymptotic profiles for various cases. We conclude with a summary and discussion in section 5.

## 2 The equations

We consider the evolution of an isolated spherical positive density perturbation in a flat universe made of collision-less matter only (density parameter $`\mathrm{\Omega }=1`$). We describe the perturbation by a large number of equal mass particles having angular momenta in random directions such that the mean angular momentum at any point in space is zero. Spherical symmetry is then preserved and the angular momentum of each particle is conserved as the system evolves. It is often convenient here to think of a particle as a a spherical shell obeying the equations of motion of a point particle moving in the gravity field of the perturbation with conserved angular momentum. We write now the equations of motion governing the evolution of a shell with a given angular momentum per unit mass, $`L`$. Denote by $`r`$ be the distance from the symmetry center, and by $`M(<r,t)`$ the mass contained inside $`r`$ at time $`t`$. The trajectory $`r(t)`$ of the shell as a function of time is then given by the solution to

$$\frac{\mathrm{d}^2r}{\mathrm{d}t^2}=\frac{GM(<r,t)}{r^2}+\frac{L^2}{r^3},$$
(1)
where the mass $`M(<r,t)`$ is determined from the distribution of shells at time $`t`$. Deferring the discussion of the way shells (particles) are assigned angular momenta, these equations determine the evolution of the system given a set of initial conditions specified at an early time $`t_i0`$. We assume that the initial radial velocity of a shell at $`r_i`$ is equal to the Hubble expansion velocity at $`t_i`$, that is,

$$\frac{\mathrm{d}r_i}{\mathrm{d}t}=\frac{2}{3t_i}r_i.$$
(2)
We express the initial mass distribution in terms of the relative mass excess $`\delta _i\mathrm{\Delta }M/M`$ interior to an initial radius $`r_i=[3M/4\pi \rho _b(t_i)]^{1/3}`$ where $`\rho _b(t_i)=(6\pi Gt)^2`$ is the mean background density $`t_i`$. We take the scale free form

$$\delta _i=\frac{\mathrm{\Delta }M}{M}=\left(\frac{M}{M_0}\right)^ϵ,$$
(3)
where $`ϵ>0`$. A perturbation with $`ϵ>1`$ can be realized by placing a point mass at the center of a void with local density contrast $`M^ϵ`$ (cf. Chuzhoy & Nusser 2000 ). In cosmology the initial perturbations are small so we restrict the analysis to shell far enough from the center so that $`\delta _i1`$.

The general solution to the equations can be obtained by numerical integration where the mass distribution is continuously updated from the new positions of the shells at each time step. However the equations allow an analytic solution for the motion of a shell that has not crossed any other so that $`M[<r(t)]=const`$. We will see later that for some choices of the angular momenta, this solution is applicable in the outer regions of the collapse. The solution can be written in the parametric form (e.g., Landau & Lifshitz 1976),

$$r=r_0(1e\mathrm{cos}\eta );t=t_0(\eta e\mathrm{sin}\eta ),$$
(4)
where

$$e^2=1+\frac{2EL^2}{G^2M^2},r_0=\frac{L^2}{GM}\frac{1}{1e^2},t_0=\sqrt{\frac{r_0^3}{GM}}..$$
(5)
According to this solution, the shell expands until it reaches a maximum expansion radius (turnaround radius) $`r_{}=r_0(1+e)`$, at time $`t_{}=\pi t_0`$. The solution ceases to be valid shortly after maximum expansion when the shell crosses other shells returning from the inner regions after passing through their maximum expansion sometime ago. The expressions for $`r_{}`$ and $`t_{}`$ for purely radial motion can be recovered by taking the limit $`L0`$ in (5). In this limit we find (cf. Peebles 1980, FG84)

$`r_{}`$ $`=`$ $`r_0(1+e)={\displaystyle \frac{L^2}{GM(1e)}}`$ (6)
$``$ $`{\displaystyle \frac{GM}{E}}=r_i\delta _i^1=\left[{\displaystyle \frac{3M_0}{4\pi \rho _b(t_i)}}\right]^{1/3}\left({\displaystyle \frac{M}{M_0}}\right)^{1/3+ϵ}`$ (7)
$`t_{}`$ $`=`$ $`{\displaystyle \frac{\pi }{2\sqrt{2}}}{\displaystyle \frac{GM}{|E|^{3/2}}}={\displaystyle \frac{3\pi }{4}}\left[{\displaystyle \frac{1}{6\pi G\rho _b(t_i)}}\right]^{1/2}\delta _i^{3/2}`$ (8)
$`=`$ $`{\displaystyle \frac{3\pi }{4}}t_i\delta ^{3/2}={\displaystyle \frac{3\pi }{4}}t_i\left({\displaystyle \frac{M}{M_0}}\right)^{3ϵ/2}.`$ (9)
The solution (4) also describes the motion of a particle with non-vanishing angular momentum in an attractive force field $`1/r^2`$. We remark here that by taking the angular momentum to zero we do not recover the motion of a particle with zero angular momentum. In the limit $`L0`$ the particle has zero radial velocity near the center oscillates between $`r=0`$ and a maximal radius, $`r_{}`$. A particle with $`L=0`$ on the other hand has an infinite radial velocity at $`r=0`$ and oscillates between $`r_{}`$ and $`r_{}`$.

### 2.1 Two schemes for assigning angular momentum

We will consider two schemes for assigning angular momenta to particles. In the first scheme (scheme A) angular momenta are assigned at the initial time in such a way that no additional scale is introduced in the collapse. The angular momentum of each particles is conserved in the subsequent evolution of the perturbation. In the second scheme (scheme B) a particle acquires its angular momentum when it is at maximum expansion, prior to that it is assumed to have a purely radial motion. The angular momentum of a particle in scheme B is $`\sqrt{GM(<r_{})r_{}}`$ per unit mass (White & Zartizky 1991), so no additional physical scale is introduced.

In scheme A the energy per unit mass of a shell can be written in terms of the initial radius $`r_i`$ and mass excess $`\delta _iM_i`$ inside the shell at the initial time, $`t_i`$,

$$E=\frac{L^2}{2r_i^2}G\frac{\delta _iM_i}{r_i},$$
(10)
so to make all energy components scale similarly we choose

$$L^2=2\alpha G\delta _iM_ir_i=2\alpha G\left[\frac{3}{4\pi \rho _b(t_i)}\right]^{1/3}M_0^{4/3}\delta _i^{4/3ϵ},$$
(11)
where $`1>\alpha >0`$ insuring that the energy $`E=(\alpha 1)G\delta M_i/r_i`$ is negative and the shell is bound. This choice for $`L`$ does not introduce any scale in the initial conditions.

The eccentricity, $`e`$, corresponding to the motion before shell crossing is, according to (5), given by

$$e^2=1+4\alpha (\alpha 1)\delta _i^2.$$
(12)
Since in cosmological perturbations the initial density contrast is small, the last relation implies that the eccentricity is very close to unity. Therefore, according to (9), the turnaround time of the shell is $`t_{}=t_i\delta ^{3/2}r_i^{3ϵ}`$ and the turnaround radius is $`r_{}=r_i\delta _i^1r_i^{1+3ϵ}`$. This scaling means that inner shells collapse earlier than outer shells and the solution (4) is always valid outside the the radius of the shell at maximum expansion at the current time, $`t`$ (hereafter, the current turnaround radius).

In scheme B, a shell moves with zero angular momentum until it reaches its turnaround radius, where it is assigned an angular momentum per units mass according to

$$L^2=^2GM_{}r_{}=^2G\left[\frac{3}{4\pi \rho _b(t_i)}\right]^{1/3}M_0^{4/3}\delta ^{4/3+ϵ},$$
(13)
where $``$ is a constant value for all shells. To insure that a shell remains bound and subsequently collapses towards the center we demand $`<1`$. Note that the angular momentum assigned to a particle at the initial time in scheme A is $`\delta _i^2M_{}r_{}`$.

Shortly after reaching its turnaround radius, a shell starts its oscillations in the collective gravitational potential well generated by shells that have passed their turnaround radii at earlier times. Angular momentum introduces an effective repellent force, $`L^3/r^3`$, preventing shells from reaching the center. We will see that in scheme A angular momentum affects the density profile only inside a radius equal to the initial radius of the shell at the current turnaround radius. For cosmological initial conditions this radius tends to zero since $`t_i0`$. So angular momentum does not play a role in fixing the density profile of the evolved perturbation. In this case the results of FG84 for collapse of particles on purely radial orbits remain valid. In scheme B, angular momentum can prevent particles from penetrating a significant fraction, depending on $``$, of the current turnaround radius.

## 3  Self-similarity and adiabatic invariance

The initial conditions are scale free, the two schemes for assigning angular momenta do not introduce any additional scale in the problem. So the only characteristic scale in the evolution of the perturbation is the scale of non-linearity which at any time $`t`$ can be taken as the current turnaround radius, $`r_t(t)`$. The mass scale corresponding to the current turnaround radius is $`M_t=M(<r_t)`$. So the mass distribution $`M(r,t)`$ must satisfy the self-similarity condition

$$M(r,t)=M_t(r/r_t),$$
(14)
where $``$ is a function of $`r/r_t`$ only. Substituting $`t_{}=t`$ in (9) we find

$`r_t`$ $`=`$ $`\left[{\displaystyle \frac{3M_0}{4\pi \rho _b(t_i)}}\right]^{1/3}\left({\displaystyle \frac{t}{t_i}}\right)^{2/3+2/(9ϵ)}`$ (15)
$`M_t`$ $`=`$ $`M_0\left({\displaystyle \frac{t}{t_i}}\right)^{2/(3ϵ)}.`$ (16)
Assume that in some region the mass distribution can be approximated as

$$M=k(t)r^\gamma ,k(t)=k_0t^s,$$
(17)
where $`k_0,s`$, and $`\gamma `$ are constants111For collapse described by a finite number of particles, the form (17) is meant to describe the mass averaged over a few particle oscillations within $`r`$.. We will determine $`\gamma `$ in the next section by generalizing the analysis of FG84 for collapse with no angular momentum (see also Zaroubi & Hoffman 1992). However, the following useful constraint

$$s+\frac{2}{3ϵ}\frac{2}{3}\gamma \left(1+\frac{2}{3ϵ}\right)=0,$$
(18)
can readily be found by substituting the asymptotic form (17) for the mass in the self-similarity condition (14) and using (16).

Throughout the rest of the paper we will assume that the potential well near the center varies adiabatically (Gunn 1977, FG84). This means that a shell near the center makes many oscillations before the potential changes significantly (Gunn 1977, FG84). In a slowly varying potential the action variables associated with the motion of a particle is invariant. The angular action variable is the angular momentum which, thanks to spherical symmetry, is conserved independent of whether or not the potential changes adiabatically. The radial action variable associated with the motion of a particle in the potential well of (17) is

$`J_r`$ $`=`$ $`2\pi {\displaystyle _{r_b}^{r_a}}dr\left({\displaystyle \frac{\mathrm{d}r}{\mathrm{d}t}}\right)`$ (19)
$`=`$ $`{\displaystyle dr\left[2\left(E\frac{\kappa (t)}{\gamma 1}r^{\gamma 1}\right)\frac{L^2}{r^2}\right]^{1/2}},`$ (20)
where $`r_a`$ and $`r_b`$ are, respectively, the apocenter and pericenter of the particle’s orbit. We prove now that the invariance of $`J_r`$ implies that the ratio $`\xi r_b/r_a`$ is constant with time. This important property will simplify the determination of $`\gamma `$ in the next section. We prove it in the following way. The energy $`E`$ and $`\kappa `$ can be expressed in terms of $`r_a`$ and $`r_b`$ as

$`{\displaystyle \frac{2\kappa }{\gamma 1}}`$ $`=`$ $`L^2{\displaystyle \frac{r_a^2r_b^2}{r_a^{\gamma 1}r_b^{\gamma 1}}},`$ (21)
$`2E`$ $`=`$ $`L^2{\displaystyle \frac{r_a^2r_b^2}{r_a^{\gamma 1}r_b^{\gamma 1}}}r_a^{\gamma 1}+{\displaystyle \frac{L^2}{r_a^2}}.`$ (22)
Therefore

$$J_r=2\pi L_\xi ^1du\left[\left(\xi ^21\right)\left(\frac{1u^{\gamma 1}}{1\xi ^{\gamma 1}}\right)+\left(1u^2\right)\right]^{1/2}$$
(23)
The angular momentum $`L`$ is conserved, so the invariance of $`J_r`$ means that $`\xi `$ is constant. Since the only special time in the particle history is its turnaround time $`t_{}`$, the invariance of $`\xi `$ implies that

$$\frac{r_{a,b}}{r_{}}=\left(\frac{t}{t_{}}\right)^q,$$
(24)
near the center. The index $`q`$ can be expressed in terms of $`\gamma `$ and $`s`$ by using (24 ) in (21). This yields

$$s=q(\gamma +1).$$
(25)
## 4 The asymptotic behavior

Assuming the asymptotic form (17) we now turn to estimating the exponents $`\gamma `$, $`s`$, and $`q`$. Let $`P(r,r_i,t)=_r^{r_a}𝑑r(\mathrm{d}r/\mathrm{d}t)^1/_{r_b}^{r_a}𝑑r(\mathrm{d}r/\mathrm{d}t)^1`$ the fraction of time a particle with initial radius $`r_i`$ spends inside the radius $`r`$ at the present time. The function $`P`$ is well defined if the mass within $`r`$ changes very little during one oscillation period of a particle inside $`r`$. The mass $`M(<r,t)`$ can be written in terms of $`P`$ as

$$M(<r,t)=_0^{M_t}\mathrm{d}M_iP(r,r_i,t),$$
(26)
where $`M_i=4\pi \rho _b(t_i)r_i^3/3`$. By substituting (14) and (17) in this last relation we obtain

$$\left(\frac{r}{r_t}\right)^\gamma =_0^{M_t}\frac{\mathrm{d}M_i}{M_t}P(r,r_i).$$
(27)
Following FG84 we write $`P`$ in terms of $`ur/r_a`$. By definition, $`P(u)=0`$ for $`u<\xi `$ and $`P(u)=1`$ for $`u>1`$. For $`\xi u1`$, we write $`P(u)=I(u)/I(1)`$, where

$$I(u)=_\xi ^udv\left[\left(\xi ^21\right)\left(\frac{1v^{\gamma 1}}{1\xi ^{\gamma 1}}\right)+\left(1v^2\right)\right]^{1/2}.$$
(28)
Using the relations (9) and (24) we express the mass, $`M_i`$, in terms of $`u`$. So

$$\left(\frac{r}{r_t}\right)^{\gamma p}=\frac{1}{p}_{r/r_t}^{\mathrm{}}duu^{(1+p)}P(u,\xi ),$$
(29)
where

$$p=\frac{6}{2+3(23q)ϵ}.$$
(30)
Particles with $`r_b>r`$, i.e., $`u=r/r_a<\xi `$, have $`P(u)=0`$ and do not contribute to the integral on the right hand side. The relation (29) allows us to obtain the asymptotic behavior in scheme A and B without a detailed calculation of $`P(u,\xi )`$.

### 4.1 Scheme A

The first step in the derivation of the asymptotic behavior is to relate $`\xi `$ to the initial density contrast. This can easily be done by substituting the expression (11) for $`L`$ in (21), and using (18) and (25). The result is

$$\frac{\xi ^2\xi ^{\gamma +1}}{\xi ^21}=\alpha (\alpha 1)(\gamma 1)\delta _i^2.$$
(31)
Since $`\delta _i1`$ this relation implies that $`\xi 1`$ as well. Retaining only the lowest order in $`\xi `$, the relation gives $`\xi \delta _i`$ for $`\gamma >1`$, and $`\xi \delta _i^{2/(\gamma +1)}`$ for $`\gamma <1`$. Since $`r_ar_{}`$, we have $`r_b=\xi r_a\xi r_{}=\xi r_i\delta ^1<r_i`$, i.e., the pericenter of a shell is smaller than its initial radius. Except near $`r_b`$, where a shell spends only a tiny fraction of its orbital time222 Consider the orbit of a particle in the plane $`(r,\psi )`$ where $`\psi `$ is the angular position. The time spent near $`r_b`$ can be estimated as $`\mathrm{\Delta }T=_{\psi _1}^{\psi _2}d\psi r^2/L<2\pi r_m^2/L`$, where $`r_m`$ is a few times $`r_b`$. By (21) we find that $`r_b^2/L\xi ^nT_c`$ where $`n>0`$ and $`T_c`$ is the period of a circular orbit at the apocenter $`r_a`$. So the time fraction spent near $`r_b`$ is negligible for $`\xi 1`$., the motion is almost radial. The initial radius, $`r_{ti}`$, of the shell currently at turnaround sets an upper limit on $`r_b`$. For $`rr_{ti}`$ all shells move on almost purely radial orbits and the analysis of FG84 for collapse on purely radial motions apply. Therefore, according to FG84, we have

$$\gamma =p=\frac{3}{1+3ϵ},s=q=0,\mathrm{f}orϵ\frac{2}{3},$$
(32)
and,

$$\gamma =1,q=\frac{3ϵ2}{9ϵ},s=2q,\mathrm{f}orϵ<\frac{2}{3}.$$
(33)
However $`r_{ti}`$ tends to zero as $`t_i0`$. So for cosmological initial conditions in which $`t_i0`$ the the collapse is identical to that of particles on purely radial motions.

### 4.2  Scheme B

The angular momentum in this scheme is given by (13). So the relation (21) gives

$$\frac{\xi ^2\xi ^{\gamma +1}}{\xi ^21}=^2(\gamma 1),$$
(34)
implying that $`\xi `$ is the same for all particles. The radius $`r_t\xi `$ is an upper limit on the pericenters of all shells. So only particles which have passed their turnaround radii early enough can contribute to the density at $`rr_t\xi `$. Particles at $`rr_t\xi ^2`$ have their apocenters inside $`r_t\xi `$ as well, and so the radius $`r_t\xi ^2`$ marks the boundary of the inner region where we expect the same power law behavior for the mass distribution. In this region only particles with $`r_b=\xi r_ar`$, i.e., $`u\xi `$ contribute to the integral in (29), so

$$\left(\frac{r}{r_t}\right)^{\gamma p}=\frac{1}{p}_\xi ^{\mathrm{}}duu^{(1+p)}P(u,\xi ).$$
(35)
The integral is independent of $`r`$, therefore we must have $`\gamma =p`$ for all $`ϵ>0`$. Using the expression (30), and the relations (18) and (25), yields

$$\gamma =p=\frac{3}{1+3ϵ},s=q=0,\mathrm{f}or\mathrm{a}llϵ>0.$$
(36)
In the collapse with $`L=0`$ (FG84) and in scheme A this result is correct only for $`ϵ2/3`$.

If $`\xi 1`$ then in the region $`r_trr_t\xi `$ particles move on almost radial orbits and the analysis of FG84 for collapse with no angular momentum is applicable. Applying the analysis of FG84 in that region gives $`\gamma =1`$ for $`ϵ2/3`$ and $`\gamma =p`$, as in (36), for $`ϵ>2/3`$. Therefore, as $`r`$ is increased, the mass profile index varies from $`\gamma =p`$ to $`\gamma =1`$ for $`ϵ2/3`$, but remains $`\gamma =p`$ for $`ϵ>2/3`$. In figure 1 we show the mass profiles obtained from spherical N-body simulations of 12000 shells for $`ϵ=0.3`$ with angular momentum introduced according to scheme B. The curves correspond, respectively, to simulations with $`=0.9`$ (solid), 0.3 (long dashed), 0.1 (short dashed), and 0.001 (dotted). To save CPU, a shell participates in the simulation only after it reaches its turnaround radius. A fourth order Runge-Kutta time integration is used to integrate the equations of motion. The mass profiles in the figure behave as predicted from our analysis. A particle with $`=0.9`$ moves in an almost circular, so that it effectively circles the halo at a radius slightly smaller than its turnaround radius. Therefore the corresponding mass profile (solid curve) would have an index close to $`\gamma =3/(1+3ϵ)=1.57`$ almost for all $`r<r_t`$. The mass profile in the inner regions for $`=0.3`$ and $`0.1`$ (long and short dashed, respectively) have $`\gamma 1.57`$, while in the outer regions it becomes $`1`$. For $`=0.001`$, particles move on almost radial orbits over the whole distance range shown in the plot. Hence $`\gamma 1`$.

## 5 discussion

We have studied the effect of non-radial motions on the mass profile of spherical self-similar collapse. Assigning angular momentum at the initial time does not affect the evolved mass profile. This is not surprising as the initial kinetic energy is dominated by the radial Hubble expansion velocity $`r/t_i`$. Any additional finite velocity is negligible compared to the Hubble expansion velocity as $`t_i0`$. This has also been shown explicitly for collapse with purely radial motions where a finite radial velocity component is superimposed on the Hubble expansion at the initial time (e.g., Peebles 1980, Padmanabhan 1993). We have also considered the consequences of a different scheme (scheme B) for assigning angular momenta to particles. In scheme B, a particle acquires its angular momentum at maximum expansion. The acquired angular momentum is chosen in such a way that the ratio of the pericenter to apocenter is the same for all particles and no additional physical scale is introduced. In this case we have confirmed the conjecture by White & Zaritsky (1991) that the scaling $`Mr^{3/(1+3ϵ)}`$ is maintained for all $`ϵ`$ as $`r0`$. But we have also shown that if the added angular momentum is small, then the behavior $`Mr`$ is restored far enough from the center for a perturbation with $`ϵ2/3`$.

What is the justification for scheme B in a generic collapse configuration? Consider a halo forming in a gravitational collapse of a high density region in an initial scale free random Gaussian. The infalling matter is distributed in clumps (satellites). If the halo is approximately spherical then the collapse process cannot introduce any additional scale and angular momentum must be proportional to the product of the mass and the radius at the epoch of maximum expansion. But scheme B also assumes that the angular momentum of a particle is negligible before it reaches maximum expansion. This is a reasonable assumption if angular momentum is mainly generated by tidal interaction between the infalling satellites themselves. Because the halo is assumed spherical, the halo-satellite tidal interaction is not expected to generate any angular momentum with respect to the halo center. The assumptions underlying scheme B are consistent with results of N-body simulations (e.g., Barnes & Efstathiou 1987). Adiabatic invariance was also applied to non scale free initial density perturbations (Hoffman & Shaham 1985, Ryden & Gunn 1987, Lokas 2000, Lokas & Hoffman 2000). N-body results imply that in these situations as well, it is reasonable to introduce angular momentum according to scheme B.

## 6 Acknowledgment

The author thanks L. Chuzhoy and Y. Hoffman for stimulating conversations.