text:
# Beyond storage capacity in a single model neuron: Continuous replica symmetry breaking

## I Introduction

In his seminal paper Hopfield  reformulated Little’s model  of auto-associative memory in terms of an energy function. By this act, the field of the statistical mechanics of neural networks was plowed and sown in, and proved itself since then remarkably fertile. The network is an interconnected set of McCulloch-Pitts neurons , the latter being perhaps the biologically least realistic model of a nerve cell. In its simplest version the model neuron can be in one of only two states, “firing” or “quiescent”, it is nonlinear, and admits a large number of connections from other units. Despite the oversimplification in its node cells, the network can exhibit complex behavior and functions as a reasonable model of content addressable memory. From the practical viewpoint, however, artificial neural networks are generically not superior to other methods in the task of storage and associative retrieval .

An important ingredient of this type of models of neural networks is what is called in statistical mechanics quenched disorder. For example, in the Little-Hopfield model the synaptic coupling strengths, made up of random numbers by the Hebb rule , can be held fixed for the duration of the neural dynamical process. This is the basis for the analogy between spin glass models, the archetypical example of a many-body problem with quenched disorder, and neural networks . Methods borrowed from the theory of spin glasses, in particular from the infinite range interaction Sherrington-Kirkpatrick model, yielded a harvest of results on a variety of neural model systems, as well as on other problems whose motivation came from outside of physics but could be formulated as disordered statistical mechanical systems .

It is safe to say that the technique inherited from spin glass theory and the most widely used in the statistical mechanical approach to neural networks in equilibrium is the replica method . It was first applied thoroughly for infinite range interaction spin glasses, thus it is especially suited for networks where neurons have a large number of connections. In its simplest version, with so called replica symmetry, it can be straightforwardly adopted to many a neural problem .

Another family of artificial neural networks consists of layered feed-forward networks , which accept a number of inputs and for a given set of couplings produce the output as a function of the inputs. Such networks were introduced for the purpose of generalization, i. e., rule extraction, from examples of input-output pairs. Then the a priori unknown target rule is to be reconstructed by adaptive changes in the couplings, called training. Since the introduction of the backpropagation algorithm  for that purpose, feed-forward networks found wide usage .

The statistical physics of neural modeling gained an impetus of lasting effect from the work by Gardner and Derrida  on the storage problem of a single neuron, called also simple perceptron in that context. This is the simplest version  of feedforward networks. While in the Little-Hopfield statistical mechanical system the quenched variables are the couplings and a microstate is a configuration of neural states, the roles in feed-forward networks are reversed. Now the space of synaptic couplings is considered as the configuration space within Boltzmannian thermodynamics and the examples appear as quenched parameters. The error for one example measures the difference between the actual output of the network and the required output. The errors on all training examples add up to form the Hamiltonian, i. e., the cost function, of the statistical mechanical model. Within the canonical statistical mechanical approach the temperature plays its usual roles. On the one hand, it is the Lagrange multiplier associated with a preset value of the error, on the other hand, it is the amplitude of noise if a gradient-descent-like dynamics of the couplings is used so as to reach optimal configuration. The thermodynamical limit is achieved by admitting a large number of adjustable couplings, but for that it is not necessary to have many neurons. The approach of Gardner and Derrida was successful in what is called equilibrium learning, when a Hamiltonian can be associated with the problem. However, statistical physical methods are proving themselves useful also in studying off-equilibrium learning algorithms .

A central quantity of a feed-forward network is its storage capacity, i. e., how many random input-output examples the network can reproduce without error. In terms of the statistical mechanical approach this is in its original formulation a zero temperature problem. The subspace of couplings that reproduce a given set of patterns is called version space, its volume, related to the ground state ($`T=0`$) entropy, vanishes beyond capacity. Since the statistical mechanical solution of the region below capacity of a single neuron by Gardner and Derrida  a number of results have been obtained about storage properties of feed-forward networks, see for example Refs. . Nevertheless, if the task is to minimize the number of incorrectly stored examples, beyond capacity the problem has not been solved. Technically this is because below capacity for completely random examples replica symmetry holds, while beyond it no finite replica symmetry breaking scheme yields thermodynamically stable solution .

In the present paper we reconsider the problem of storage of random patterns, technically generalize Parisi’s solution of the Sherrington-Kirkpatrick model, and obtain beyond capacity a phase reminiscent to the frustrated ground state of the Sherrington-Kirkpatrick model. That phase continues for $`T>0`$ into the analog of the low temperature, or Parisi, phase of the spin glass.

Our work is motivated not only by the problem of storing random patterns. Generalization has also been successfully analyzed by statistical mechanical methods, see Ref.  for reviews. The storage problem below capacity is analogous to equilibrium learning of a learnable task, where the network is compatible with all possible examples, there is no frustration in either systems. For instance, equilibrium generalization properties of the perceptron when the examples are generated by another one, can be understood for arbitrary number of examples within replica symmetry . However, a given feed-forward network may be unable to reproduce a complicated target function on all possible examples . If the target is unlearnable then the network is presumed to get into a frustrated phase if a sufficiently large number of examples are used. Thus the properties of a network beyond capacity are of foremost interest from the viewpoint of rule extraction as well.

The single neuron may be considered as the hydrogen atom of neural problems and studied for its own interest. It is the unit of the Little-Hopfield network, where the symmetry in the couplings has been given up and the couplings of different neurons are considered as independent. As to a feed-forward network, even if the whole network operates without error, its units may still be strained beyond their individual capacities . Thus the description of a single neuron beyond its storage capacity is of importance also from the viewpoint of networked neurons. Furthermore, a close analogy exists between the behavior of model neurons beyond capacity and the glassy, frustrated, phase of disordered spin systems . Therefore, the understanding of the way a single neuron works may have ramifications beyond the field of artificial neural networks.

We firstly review in Sec. II the statistical mechanics of storage, recall the basic thermodynamical quantities and the formula for the replica free energy of a single neuron. In Sec. III we summarize the main ingredients of Parisi’s approach and obtain the free energy functional that we propose describes the equilibrium problem. This way some background is given to our previous communication , wherein we identified a spin glass phase of Parisi type in the high temperature limit. The recipe for the calculation of expectation values by means of Green functions is explained in Sec. IV, producing among other the formulas for the local stability distribution and the stationarity conditions. Sec. V is devoted to the scaling for low temperatures, enabling us to put the extremization of the free energy functional on a computer, and in the end a simulation is discussed.

## II The storage problem and its replica free energy

The model neuron, or perceptron, under consideration is

$`\xi `$ $`=`$ $`\text{sign}(h),`$ (2)
$`h`$ $`=`$ $`N^{1/2}{\displaystyle _{k=1}^N}J_kS_k,`$ (3)
where $`𝐉`$ is the synaptic coupling vector, $`𝐒`$ the input and $`\xi `$ the output. Patterns are given as input-output data,

$$\{𝐒^\mu ,\xi ^\mu \}_{\mu =1}^M.$$
(4)
In the simplest setup the $`S_k^\mu `$-s are independently drawn from any distribution with unit variance and zero average, and $`\xi ^\mu =\pm 1`$ with equal probability. We introduce the local stability parameter as

$`\mathrm{\Delta }^\mu =h^\mu \xi ^\mu ,`$ (5)
where $`h^\mu `$ is given by (3) with $`S_k^\mu `$. If the neuron generates $`\xi ^\mu `$ in response to $`𝐒^\mu `$, i. e.$`\mathrm{\Delta }^\mu >0`$, then we say that the $`\mu `$-th pattern is stored. If $`\mathrm{\Delta }^\mu `$ is a large positive number then high stability of storage against changes in either the couplings or the inputs can be assumed. Large stability is associated with large basin of attraction in memory networks . Given the ensemble of patterns, the local stability parameter obeys some distribution $`\rho (\mathrm{\Delta })`$ . If the number of patterns $`M`$ is of order $`N`$ then it is useful to introduce the relative number of examples

$$\alpha =M/N,$$
(6)
called also load parameter. Since an overall positive factor of $`𝐉`$ does not change the output, we set the norm of $`𝐉`$ to $`\sqrt{N}`$, expressed by the prior distribution

$$w(𝐉)=C_N\delta \left(N\left|𝐉\right|^2\right).$$
(7)
This is called spherical constraint. The factor $`C_N`$ normalizes $`w(𝐉)`$ to unity, it has no thermodynamical significance besides setting the zero point of the entropy scale. Given the distribution of patterns and the length of $`𝐉`$ it can be easily seen that the normalization in (3) results in $`h`$ values of typically $`O(1)`$.

Storage with minimal error can be formulated as an optimization task by our introducing an error measure. If we treat all patterns in the same way we obtain what is called the equilibrium problem. The associated Hamiltonian, or cost function, is

$$=\underset{\mu =1}{\overset{M}{}}V(\mathrm{\Delta }^\mu ),$$
(8)
where the potential $`V(\mathrm{\Delta }^\mu )`$ gives the error on a single pattern $`𝐒^\mu ,\xi ^\mu `$. Obviously, $`V(y)=0`$ for $`y>0`$ in the original storage problem. One can also impose a bound $`\kappa `$ for the local stability, i. e.$`V(y)`$ is set to be zero for $`y>\kappa `$, if $`\kappa >0`$ this is obviously stricter than the original storage criterion. Generically, $`V(y)`$ should be monotonically decreasing for $`y<\kappa `$. In this paper specific results will be presented on the error counting, or Gardner-Derrida , potential

$$V(y)=\theta (\kappa y),$$
(9)
where $`\theta (y)`$ is the Heaviside function. Given the potential (9) the aim is to minimize the number of patterns whose stability is below the bound $`\kappa `$. In the theoretical framework we shall keep the general form $`V(y)`$.

The partition function of the optimization task is

$$Z=d^NJw(𝐉)\mathrm{exp}\left(\beta \underset{\mu =1}{\overset{M}{}}V(\mathrm{\Delta }^\mu )\right),$$
(10)
with $`\beta =1/T`$. Quenched average, $`\mathrm{}_{\text{qu}}`$, is defined as the mean over the patterns. We deal with the idealized equilibrium of the system, when for large $`N`$ the free energy, the energy, and the entropy are assumed to approach their quenched average. This property of self-averaging was proved rigorously only in special cases, see for example , but it is the widely used basis in studies of the equilibrium thermodynamics of disordered systems .

The replica method  consists in writing the mean free energy per coupling as

$$f=\underset{N\mathrm{}}{lim}\frac{\mathrm{ln}Z_{\text{qu}}}{N\beta }=\underset{N\mathrm{}}{lim}\underset{n0}{lim}\frac{1Z^n_{\text{qu}}}{nN\beta }.$$
(11)
Denoting the thermal average with the Boltzmann weight in (10) by $`\mathrm{}_{\text{th}}`$ the mean error per pattern can be written as

$$\epsilon =V(\mathrm{\Delta })_{\text{th}}_{\text{qu}}=\frac{1}{\alpha }\frac{\beta f}{\beta }.$$
(12)
The entropy is

$$s=\beta (\alpha \epsilon f).$$
(13)
For $`T=0`$ and $`\epsilon =0`$ the volume of version space, i. e., the space of couplings that perfectly reproduce the examples is obtained as $`\mathrm{\Omega }=\mathrm{exp}(Ns)`$. In general, $`Ns`$ has the usual meaning of the logarithm of the volume with given error $`\epsilon `$.

Introducing the overlap matrix $`𝖰`$ of synaptic vectors

$$\left[𝖰\right]_{ab}q_{ab}=\frac{1}{N}\underset{k=1}{\overset{N}{}}J_{ak}J_{bk},$$
(14)
where the replica indices $`a,b`$ go from $`1`$ to $`n`$, we can express the free energy as the result of the minimization of the replica free energy

$`f`$ $`=`$ $`\underset{n0}{lim}{\displaystyle \frac{1}{n}}\underset{𝖰}{\mathrm{min}}f(𝖰)`$ (16)
$`f(𝖰)`$ $`=`$ $`f_s(𝖰)+\alpha f_e(𝖰),`$ (17)
$`f_s(𝖰)`$ $`=`$ $`(2\beta )^1\mathrm{ln}\text{det}𝖰,`$ (18)
$`f_e(𝖰)`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}\mathrm{ln}{\displaystyle \frac{d^nxd^ny}{(2\pi )^n}\mathrm{exp}\left(\beta _{a=1}^nV(y_a)+i\mathrm{𝐱𝐲}\frac{1}{2}𝐱𝖰𝐱\right)}.`$ (19)
The subscripts $`s`$ and $`e`$ stand for entropic and energy-like terms. The entropic contribution (18) arises because of the spherical constraint and the definition (14), it is indeed independent of the potential, while (19) depends on it.

A central role is played by the probability density for the local stabilities

$$\rho (\mathrm{\Delta })=\delta \left(\mathrm{\Delta }h^1\xi ^1\right)_{\text{th}}_{\text{qu}},$$
(20)
where (3) is understood. Due to the symmetry of the Hamiltonian with respect to the permutation of patterns we could choose $`\mu =1`$ for convenience. It will turn out to be useful to interpret the integrand in (19) as an effective Boltzmann weight and denote the average over this measure as

$$\mathrm{},$$
(21)
where the $`n0`$ limit is implied. A straightforward replica calculation shows (see Ref.  for a pedagogic presentation) that the local stability distribution can be rewritten as

$$\rho (\mathrm{\Delta })=\delta \left(\mathrm{\Delta }y_1\right).$$
(22)
Here the subscript could be any replica index, for convenience we chose $`1`$. Comparison with (20) allows an intuitive interpretation for the replica average $`\mathrm{}`$, namely, this corresponds to the combined average over thermal and quenched fluctuations. From (20,22) it follows that the combined thermal and quenched average in (12) boils down to

$$\epsilon =V\left(y_1\right)=𝑑y\rho (y)V(y).$$
(23)
The free energy (14) was calculated within the replica symmetric ansatz for the error counting potential (9) and the capacity, i. e., the maximal $`\alpha `$ with $`\epsilon =0`$ at $`T=0`$ was determined in . It has been shown that beyond capacity the replica symmetric solution is thermodynamically unstable . One  and two  step replica symmetry breaking solutions were presented, while Ref.  proved that no finite step symmetry breaking ansatz can possibly be thermodynamically stable. We presented in  a variational free energy functional without derivation that incorporated continuous replica symmetry breaking, but gave concrete results only in the high temperature, large $`\alpha `$ limit. In what follows we provide some background to the general theory and in the end properties of the ground state ($`T=0`$) beyond capacity will also be described.

## III The Parisi scheme

In this section we show how to evaluate the replica free energy $`f_e(𝖰)`$ of (19) within Parisi’s ansatz. The $`R`$ step replica symmetry breaking form is

$$𝖰=\underset{r=0}{\overset{R+1}{}}\left(q_rq_{r1}\right)𝖴_{m_r}𝖨_{n/m_r},$$
(24)
where $`k`$ subscripts $`k`$-dimensional matrices, $`𝖨_k`$ is the identity operator, all elements of $`𝖴_k`$ equal $`1`$, $``$ marks the direct product, and

$`q_1`$ $`=`$ $`0q_0q_1\mathrm{}q_Rq_{R+1}=1,`$ (26)
$`m_{R+1}`$ $`=`$ $`1m_Rm_{R1}\mathrm{}m_1m_0=n,`$ (27)
where the integer $`m_r`$ is a divisor of $`m_{r1}`$. The $`n0`$ limit can be performed smoothly if instead of $`m_r`$ we use $`x_r=(nm_r)/(n1)`$ for parametrization. Thus for arbitrary $`n>0`$ we have the ordering

$$x_{R+1}=1x_Rx_{R1}\mathrm{}x_1x_0=0.$$
(28)
We consider the $`x_r`$-s fixed along the $`n0`$ limiting process, whence follows the formal $`n`$-dependence of the $`m_r(n)`$-s, and for $`n=0`$ we get $`x_r=m_r(0)`$. The inspection of the first few $`R=0,1,2`$ cases  allows, in the spirit of Parisi’s , the generalization of the energy term (19) in the replica free energy to arbitrary $`R`$ as

$`f_e(𝐪,𝐱)`$ $`=`$ $`\underset{n0}{lim}{\displaystyle \frac{1}{n}}f_e(𝖰)`$ (29)
$`=`$ $`{\displaystyle \frac{1}{\beta x_1}}{\displaystyle Dz_0\mathrm{ln}Dz_1}`$ (31)
$`\times \left[{\displaystyle Dz_2\mathrm{}\left[Dz_{R+1}\mathrm{exp}\left\{\beta V\left(\underset{r=0}{\overset{R+1}{}}z_r\sqrt{q_rq_{r1}}\right)\right\}\right]^{\frac{x_R}{x_{R+1}}}\mathrm{}}\right]^{\frac{x_1}{x_2}},`$
where

$$Dz=\frac{dze^{\frac{1}{2}z^2}}{\sqrt{2\pi }}.$$
(32)
This is the analog of Parisi’s formula for the Sherrington-Kirkpatrick model, Eq. (11) in ; a comprehensive derivation will be presented elsewhere . The energy term (19) has become a function of the parameters in (26,28). The evaluation of (29) can be done by iteration,

$`\psi _{r1}(y)`$ $`=`$ $`{\displaystyle Dz\psi _r\left(y+z\sqrt{q_rq_{r1}}\right)^{\frac{x_r}{x_{r+1}}}},`$ (34)
$`\psi _{R+1}(y)`$ $`=`$ $`e^{\beta V(y)},`$ (35)
where $`x_{R+2}=1`$ is understood. Then the sought free energy term is obtained as

$$f_e(𝐪,𝐱)=\frac{1}{\beta x_1}Dz\mathrm{ln}\psi _0\left(z\sqrt{q_0}\right).$$
(36)
where $`𝐪=(q_0,\mathrm{},q_R)`$ and $`𝐱=(x_1,\mathrm{},x_R)`$.

The above iteration can be redressed as a partial differential equation. Parisi’s order parameter function $`x(q)`$ is a concatenation of the $`𝐪`$ and $`𝐱`$ as

$$x(q)=\underset{r=0}{\overset{R}{}}(x_{r+1}x_r)\theta (qq_r),$$
(37)
where $`x_1=0`$. Next we introduce $`\psi (q,y)`$ such that at the discontinuities

$`\psi (q_r^{+0},y)`$ $`=`$ $`\psi _r(y),`$ (39)
$`\psi (q_r^0,y)`$ $`=`$ $`\psi (q_r^{+0},y)^{\frac{x(q_r^0)}{x(q_r^{+0})}},`$ (40)
that is, $`\psi (q,y)^{1/x(q)}`$ is continuous in $`q`$. Furthermore, along the plateau in the open interval $`(q_{r1},q_r)`$

$$\psi (q,y)=Dz\psi (q_r^0,y+z\sqrt{q_rq}).$$
(41)
It is easy to show that the field so defined satisfies the partial differential equation (PDE)

$`_q\psi (q,y)`$ $`=`$ $`{\displaystyle \frac{1}{2}}_y^2\psi (q,y)+{\displaystyle \frac{\dot{x}(q)}{x(q)}}\psi (q,y)\mathrm{ln}\psi (q,y),`$ (43)
$`\psi (1,y)`$ $`=`$ $`e^{\beta V(y)},`$ (44)
which evolves from $`q=1`$ to $`q=0`$. Indeed, along the plateaus $`\dot{x}(q)=0`$ when only the first term on the r. h. s. of (43) remains, thus producing (41). Near jumps of $`x(q)`$ the second term dominates, and at a fixed $`y`$ the resulting ordinary differential equation in the variable $`q`$ is separable. Hence it follows that $`\psi (q,y)^{1/x(q)}`$ is continuous in $`q`$, thus (40) is recovered. An equivalent field can be defined by

$$f(q,y)=\frac{\mathrm{ln}\psi (q,y)}{\beta x(q)},$$
(45)
satisfying for a continuous potential $`V(y)`$ the PDE

$`_qf(q,y)`$ $`=`$ $`{\displaystyle \frac{1}{2}}_y^2f(q,y)+{\displaystyle \frac{1}{2}}\beta x(q)\left(_yf(q,y)\right)^2,`$ (47)
$`f(1,y)`$ $`=`$ $`V(y).`$ (48)
The fact that $`f`$ denotes the free energy (16) as well as the field $`f(q,y)`$ should not cause misunderstandings. Equation (47) with initial condition $`\mathrm{ln}2\mathrm{cosh}\beta y`$ has been discovered by Parisi  while studying the Sherrington-Kirkpatrick model.

If the potential $`V(y)`$ is not continuous, the PDE (45) holds only from any $`q^{}<1`$ onward where $`\mathrm{ln}\psi (q^{},y)`$ is continuous in $`y`$. In the generic case such is $`q_R`$, so the evolution along the first plateau from $`1`$ to $`q_R`$, where $`x(q)1`$, is to be done explicitly as

$$f(q_R,y)=\frac{\mathrm{ln}\psi (q_R,y)}{\beta x_{R+1}}=\beta ^1\mathrm{ln}Dze^{\beta V\left(y+z\sqrt{1q_R}\right)},$$
(49)
and for $`q<q_R`$ the PDE (47) with the initial condition (49) can be used. Alternatively, one can introduce the field $`m(q,y)`$ as

$$\psi (q,y)m(q,y)=\frac{_y\psi (q,y)}{\beta x(q)},$$
(50)
which equals $`_yf(q,y)`$ when continuity in $`y`$ holds. Then the PDE (47) should be replaced by

$$_qf(q,y)=\frac{1}{2}_ym(q,y)+\frac{1}{2}\beta xm^2(q,y),$$
(51)
while the initial condition (48) can be kept.

Finally we obtain the sought free energy term (19) as a functional of the order parameter function

$$f_e[x(q)]=\underset{n0}{lim}\frac{1}{n}f_e(𝖰)=f(0,0).$$
(52)
It should be emphasized that the above PDE-s do not require infinite refining of the partition by $`q_r`$-s of the interval $`(0,1)`$. They are valid for discrete as well as continuous replica symmetry breaking schemes, i. e., they admit $`x(q)`$ with steps and plateaus, as well as strictly monotonically increasing continuous segments.

The entropic term (18) can also be cast in the form of the energy term (19) with the substitution

$$e^{\beta V(y)}=\sqrt{2\pi }\delta (y).$$
(53)
If we conceive the Dirac delta as a Gaussian with small variance, the initial condition (48) becomes a quadratic function and the PDE (47) can be solved analytically. The analogue of (52) for the entropic term, after going with the variance to zero, can be cast into

$$f_s[x(q)]=\underset{n0}{lim}\frac{1}{n}f_s(𝖰)=\frac{1}{2\beta }_0^1𝑑q\left[\frac{1}{D(q)}\frac{1}{1q}\right],$$
(54)
where

$$D(q)=_q^1𝑑\overline{q}x(\overline{q}).$$
(55)
The free energy of the neuron is then obtained as

$$f=\begin{array}{c}\text{ }\\ \text{max}\\ \text{x(q)}\end{array}f\left[x(q)\right],$$
(56)
where the free energy functional is

$$f[x(q)]=f_s[x(q)]+\alpha f_e[x(q)],$$
(57)
with $`f_s`$ and $`f_e`$ defined in Eqs. (54,52). It is due to the $`n0`$ limit that the maximization in (56) replaces the minimization by the matrix elements of $`𝖰`$ in (16), see, e. g., Ref. .

## IV Linear response theory, stationarity conditions, and expectation values

The least obvious part of the extremization condition (56) is the variation of $`f(0,0)`$ by $`x(q)`$. This can be calculated from linear response theory for the PDE (47). Moreover, linear response theory yields a technique to calculate replica averages as introduced in (21), essential for the evaluation of physical quantities.

The Green function for the PDE (47) can be introduced formally as

$$𝒢(q_1,y_1;q_2,y_2)=\frac{\delta f(q_1,y_1)}{\delta f(q_2,y_2)},$$
(58)
whence $`𝒢(q_1,y_1;q_2,y_2)=0`$ for $`q_1>q_2`$. The Green function for the Parisi solution of the Sherrington-Kirkpatrick model has been studied in Refs. . In the fore and hind variable pairs the Green function $`𝒢(q_1,y_1;q_2,y_2)`$ satisfies the respective PDE-s

$`_{q_1}𝒢={\displaystyle \frac{1}{2}}_{y_1}^2𝒢+\beta x(q_1)m(q_1,y_1)_{y_1}𝒢\delta (q_1q_2)\delta (y_1y_2),`$ (60)
$`_{q_2}𝒢={\displaystyle \frac{1}{2}}_{y_2}^2𝒢+\beta x(q_2)_{y_2}\left[m(q_2,y_2)𝒢\right]+\delta (q_1q_2)\delta (y_1y_2),`$ (61)
where $`m`$ is given in (50). The first equation without the Dirac delta excitation is the linearization of the PDE (47). The minus sign of the Dirac deltas follows from the fact that (60) evolves towards decreasing “time” $`q`$. The homogeneous part of (61) is obtained from the requirement that

$$𝒢(q_1,y_1;q_3,y_3)=𝑑y_2𝒢(q_1,y_1;q_2,y_2)𝒢(q_2,y_2;q_3,y_3)$$
(62)
does not depend on $`q_2`$, and the plus sign of the inhomogeneous term is due to the fact that evolution goes towards increasing $`q`$. The homogeneous parts of two PDE-s (61) and (61) are called adjoint to each other.

For the sake of simplicity we gave formula (58) for the case of continuous potential $`V(y)`$. If the potential is discontinuous then the definition (58) should and can be appropriately modified when a $`q`$-argument is near $`1`$, but the PDE-s (58) for the Green functions hold as they are.

The significance of the Green function is in that it helps to solve the linear PDE with the source term $`h(q,y)`$

$$_q\vartheta (q,y)=\frac{1}{2}_y^2\vartheta (q,y)+\beta x(q)m(q,y)_y\vartheta (q,y)+h(q,y),$$
(63)
as

$$\vartheta (q,y)=𝑑y_1𝒢_\phi (q,y;1,y_1)\vartheta (1,y_1)_q^1𝑑q_1𝑑y_1𝒢_\phi (q,y;q_1,y_1)h(q_1,y_1).$$
(64)
A prominent role will be played by

$$P(q,y)=𝒢(0,0;q,y),$$
(65)
which solves the PDE (61) with $`q_1=y_1=0`$, i. e., with initial condition

$$P(0,y)=\delta (y).$$
(66)
This function first appeared in the context of the Sherrington-Kirkpatrick model in Ref. . Note that the PDE for $`P(q,y)`$ is in fact a Fokker-Planck equation, producing a nonnegative solution and conserving the norm $`𝑑yP(q,y)1`$ for all $`q`$-s. This suggests the intuitive interpretation of $`P(q,y)`$ as probability density of $`y`$.

Now we are in the position to calculate the variation of the free energy functional (57). As to the energy term (52), by varying the functions $`f(q,y)`$ and $`x(q)`$ in the PDE (47) one obtains (63) with $`\vartheta =\delta f`$, $`\vartheta (1,y)=0`$, and $`h=\frac{1}{2}\beta (_yf)^2\delta x`$. Hence (64) gives at $`q=0,y=0`$ the sought $`\delta f(0,0)/\delta x`$. The variation of $`f_s[x(q)]`$ can be calculated straightforwardly, and, with the notation (65), we arrive at

$$F(q,[x(q)])=\frac{2}{\beta }\frac{\delta f[x(q)]}{\delta x(q)}=\left(_0^q\frac{d\overline{q}}{\beta ^2D(\overline{q})^2}\alpha 𝑑yP(q,y)m(q,y)^2\right).$$
(67)
The stationarity condition in case $`x(q)`$ can be freely varied is thus

$$F(q,[x(q)])=0.$$
(68)
If in an interval $`I`$ the $`x(q)`$ is supposed to have a plateau, we differentiate by the plateau value of $`x(q)`$ to get

$$_I𝑑xF(q,[x(q)])=0.$$
(69)
In isolated points $`q_r`$ where plateaus meet (68) should hold pointwise. This summarizes the stationarity conditions for an arbitrary order parameter function, i. e., arbitrary replica symmetry broken scheme of Parisi type.

We have seen in Sec. II instances when the double, thermal and quenched, average could be replaced by the replica average (21). Replica averages can be calculated by the Green function technique as described below. The procedure can be viewed as the generalization of the groundbreaking results from Refs. , where the local magnetization and some of its moments in the Parisi phase of the Sherrington-Kirkpatrick model were evaluated.

A simple case is when $`A(y_a)`$ is to be calculated for an arbitrary function $`A(y)`$. Because of the symmetry with respect to the permutation of single replica indices we have $`A(y_a)=n^1_{a=1}^nA(y_a)`$. This quantity can be easily evaluated if one replaces in (19) $`V(y)`$ by $`V(y)+\lambda A(y)`$, thus obtains $`f_e(𝖰;\lambda )`$, and calculates its initial slope at $`\lambda =0`$. Reversing the limits $`n0`$ and $`\lambda 0`$ then using the first equality of (29) and Eq. (52) we get

$`A(y_a)`$ $`=`$ $`\underset{n0}{lim}{\displaystyle \frac{1}{n}}{\displaystyle \frac{f_e(𝖰;\lambda )}{\lambda }}|_{\lambda =0}={\displaystyle \frac{f(0,0;\lambda )}{\lambda }}|_{\lambda =0}`$ (70)
$`=`$ $`{\displaystyle 𝑑y\frac{\delta f(0,0)}{\delta f(1,y)}A(y)}={\displaystyle 𝑑yP(1,y)A(y)}.`$ (71)
The third equality comes from the fact that $`\lambda `$ is in the initial condition for $`f(q,y)`$ at $`q=1`$, and the last one comes from the definitions for the Green function (58) and for $`P(q,y)`$ (65). Again, for the sake of brevity we gave the derivation for continuous potential $`V(y)`$, however, the result holds also for discontinuous ones. Immediately follows from (22) the formula for the probability density of the local stabilities

$$\rho (y)=P(1,y),$$
(72)
and thus from (23)

$$\epsilon =𝑑yP(1,y)V(y).$$
(73)
From the practical viewpoint interesting are effective averages of products like

$$A_1(y_{a_1})A_2(y_{a_2})\mathrm{}A_k(y_{a_k}).$$
(74)
One can show by using elementary properties of the Fourier transformation that the average of a product of $`x_a`$-s over the effective Boltzmann weight on the r. h. s. of (19) can be expressed as averages over functions of variables $`y_a`$-s. Thus the knowledge how to evaluate (74) also resolves the problem of averages of polynomials in $`x_a`$-s. The latter quantities are of importance because they appear when the replica free energy is differentiated in terms of $`q_{ab}`$-s.

Here we shall only describe the recipe for calculating (74), details can be found in Ref. . If $`k=2`$ then the average depends on $`q=q_{a_1a_2}`$ and is given by the formula

$`C_{12}(q)`$ $`=`$ $`A_1(y_{a_1})A_2(y_{a_2})`$ (75)
$`=`$ $`{\displaystyle 𝑑y𝑑y_2𝑑y_3P(q,y)𝒢(q,y;1,y_1)𝒢(q,y;1,y_2)A_1(y_1)A_2(y_2)}.`$ (76)
Of such type is $`f_e(𝖰)/q_{a_1a_2}`$ where $`A_1(y)=A_2(y)=im(1,y)`$, see (50) for definition, whence we obtain the second term in $`F(q,[x(q)])`$ given in Eq. (67). This is related to the fact that the stationarity condition can also be obtained by first differentiating the replica free energy (17) by $`q_{ab}`$ and then equating the result to zero. For $`k=3`$ suppose without restricting generality that

$$q=q_{a_1a_3}=q_{a_2a_3}<\overline{q}=q_{a_1a_2}.$$
(77)
Then we have

$`C_{123}(q,\overline{q})`$ $`=`$ $`A_1(y_{a_1})A_2(y_{a_2})A_3(y_{a_3})`$ (78)
$`=`$ $`{\displaystyle 𝑑y𝑑\overline{y}𝑑y_1𝑑y_2𝑑y_3P(q,y)𝒢(q,y;1,y_3)𝒢(q,y;\overline{q},\overline{y})𝒢(\overline{q},\overline{y};1,y_1)𝒢(\overline{q},\overline{y};1,y_2)}`$ (80)
$`\times A_1(y_1)A_2(y_2)A_3(y_3).`$
Special versions of the above formulas, for the case of the second and third moments of the magnetization in the Sherrington-Kirkpatrick model, were worked out in Refs. . The integrals in (78) admit a simple graphic representation as shown on Fig. 1.

The cases considered above are the effective average for $`k=1`$, given by Eq. (70), represented by one line, and for $`k=2`$, as calculated in (76), represented by a fork with one handle and two branches. Averages of more than three functions can be analogously constructed, for a given $`k`$ a graph has $`k+1`$ “legs”. Obviously there are two topologically possible graphs for $`k=4`$, depending on the overlaps $`q_{a_ia_j}`$, and more for larger $`k`$-s.

The ability to calculate $`k=4`$ effective averages allows us to study linear stability of the replica free energy (17) at the stationary $`𝖰`$. Using the results of Ref.  on ultrametric matrices we expressed the so called replicon eigenvalues in terms of Green functions. While a general proof of the fact that there are no negative eigenvalues in the case of continuous replica symmetry breaking, i. e., when $`x(q)`$ has a continuously increasing segment, is not available, in the high temperature limit  we confirmed the absence of linear instability against replicons whenever we encountered such a stationary state. For any temperatures we recovered analytically the zero eigenvalues, corresponding to Goldstone modes, as well as the lowest order Ward-Takahashi identities predicted by algebra .

The generalization of (74) to non-factorizable functions is straightforward. For example, such functions would simply replace the products of $`A_k`$-s in (76) and (78).

## V Low temperature results

In Ref.  we have shown that in the high temperature limit, i. e., for $`\alpha ,\beta \mathrm{}`$ with $`\gamma =\alpha \beta ^2`$ finite, the problem simplifies to the extent that if $`x(q)`$ has a continuously increasing segment, this can be given in a closed analytic form. In that limit the problem becomes equivalent to the spherical, multi-$`p`$-spin interaction spin glass , where a similar observation has been made. Four different phases have been found  for the error counting potential (9): for small $`\gamma `$ replica symmetry holds, and for $`|\kappa |<2`$ and large $`\gamma `$ there is a Parisi phase with a single continuously increasing segment of $`x(q)`$ between the trivial plateaus $`x0`$ and $`x1`$. When $`|\kappa |>2`$ there is also a narrow one-step replica symmetry broken regime, and for large enough $`\gamma `$-s equilibrium is characterized by an $`x(q)`$ that is a concatenation of a nontrivial plateau and a continuously increasing segment. While for $`T>0`$ there cannot be error free storage, it is plausible to conceive the replica symmetric regime as the continuation of the $`T=0`$ phase of perfect storage and the symmetry broken phases as the analog of the regime at $`T=0`$ beyond capacity.

In the case of a $`V(y)`$ potential that vanishes for $`y>\kappa `$ the limit of capacity is given for $`\kappa 0`$ at $`T=0`$ by

$$\alpha _c(\kappa )=\left(_{\mathrm{}}^\kappa Dt(\kappa t)^2\right)^1,$$
(81)
as it follows from the replica symmetric solution when $`q1`$. This formula also gives the limit of the de Almeida-Thouless (AT) local stability  in the case of the potential (9). For $`\kappa =0`$ one has $`\alpha _c=2`$, and, for increasing $`\kappa `$, $`\alpha _c(\kappa )`$ understandably decreases. As it has been already mentioned, beyond capacity none of the finite-step replica symmetry breaking schemes gives a locally stable equilibrium state . Thus in this regime the order parameter function is no longer of the step-like form of (37), rather it has a continuously increasing part. This makes it necessary to numerically solve the extremization problem (56).

The ground state (T=0) has its special scaling properties. The PDE (47) stays meaningful if for $`q<1`$ the function $`\beta x(q)`$ does not diverge, implying that $`x(q)`$ goes to zero. Given the meaning of $`x(q)`$ as the probability of a $`q`$ being in the interval $`(0,q)`$ , we can say that at $`T=0`$ the overlap $`q`$ is $`1`$ with probability 1 for all $`\alpha >\alpha _c`$. Nevertheless, a physically meaningful order parameter is obtained after scaling by $`\beta `$. Firstly we introduce for $`T>0`$ the parameters of the classic Parisi shape as

$$\begin{array}{cc}x(q)0\hfill & \text{ if }0q<q_{(0)}\hfill \\ \dot{x}(q)>0\text{ and }x_{(0)}<x(q)<x_{(1)}\hfill & \text{ if }q_{(0)}<q<q_{(1)}\hfill \\ x(q)1\hfill & \text{ if }q_{(1)}<q1.\hfill \end{array}$$
(82)
The scaled quantities

$`q(t)`$ $`=`$ $`q_{(1)}\left(q_{(1)}q_{(0)}\right)\left(1(1+q_{(1)})t+q_{(1)}t^2\right),0t1,`$ (84)
$`\xi (t)`$ $`=`$ $`\beta x(q(t))\dot{q}(t),`$ (85)
$`\eta `$ $`=`$ $`\beta \left(1q_{(1)}\right),`$ (86)
$`\mathrm{\Delta }(t)`$ $`=`$ $`\beta D(q(t))={\displaystyle _t^1}\xi (\overline{t})𝑑\overline{t}+\eta ,`$ (87)
are expectedly regular even in the $`T0`$ limit, when $`q_{(1)}1`$. Note that there is an arbitrariness in the parametrization by $`t`$, the main features being that $`q(0)=q_{(0)}`$, $`q(0)=q_{(0)}`$, $`q(1)=q_{(1)}`$, and $`\dot{q}(1)=0`$. With this parametrization the PDE (45) becomes

$`_tf(t,y)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{q}(t)_y^2f(t,y)+{\displaystyle \frac{1}{2}}\xi (t)\left(_yf(t,y)\right)^2,`$ (89)
$`f(1,y)`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}\mathrm{ln}{\displaystyle Dze^{\beta V\left(y+z\sqrt{1q_{(1)}}\right)}}.`$ (90)
For $`T=0`$ the initial condition becomes

$$f(t=1,y)|_{T=0}=\underset{\overline{y}}{\mathrm{min}}\left(V(\overline{y})+\frac{(y\overline{y})^2}{2\eta }\right)=\{\begin{array}{cc}1\hfill & \text{if }y\kappa \sqrt{2\eta }\hfill \\ \frac{(\kappa y)^2}{2\eta }\hfill & \text{if }\kappa \sqrt{2\eta }y\kappa \hfill \\ 0\hfill & \text{if }y\kappa ,\hfill \end{array}$$
(91)
where (9) was substituted to get the second equality. By Gaussian integration hence the replica symmetric solution can be obtained . Note that $`f(t=1,y)`$ is a continuous function even though $`V(y)`$ is the step function. Although we used the same symbols for functions of $`q`$ and $`t`$, misunderstanding are avoided by our marking which argument we mean.

The PDE for $`P(q,y)`$ follows from the definition (65) and from the evolution equation of the Green function (61). The latter can be properly rescaled according to (82), yielding in principle the solution $`P(t,y)`$. The probability density of local stabilities is obtained by evolving $`P(t=1,y)=P(q=q_{(1)},y)`$ to $`P(q=1,y)`$. At $`T=0`$ with (9) we have

$$\rho (\mathrm{\Delta })=\{\begin{array}{cc}P(1,\mathrm{\Delta })\hfill & \text{if }\mathrm{\Delta }\kappa \sqrt{2\eta }\hfill \\ 0\hfill & \text{if }\kappa \sqrt{2\eta }<\mathrm{\Delta }<\kappa \hfill \\ P(1,\mathrm{\Delta })+\delta (\mathrm{\Delta }\kappa )_{\kappa \sqrt{2\eta }}^\kappa 𝑑\overline{y}P(1,\overline{y})\hfill & \text{if }\mathrm{\Delta }\kappa ,\hfill \end{array}$$
(92)
where $`P(1,\mathrm{\Delta })=P(t=1,\mathrm{\Delta })`$ is understood. For arbitrarily small $`T>0`$ the gap in the support of $`\rho (\mathrm{\Delta })`$ immediately vanishes. For details we again refer to .

The numerical extremization was done with complementing the free energy functional (57) by constraints. The PDE (V) was added giving rise to a Lagrange multiplier field. This field can be shown to be just $`P(t,y)`$, see Refs. . Further technical requirements are $`x(q_{(0)})0`$, $`x(q_{(1)})1`$, and $`\dot{x}(q)0`$ for $`q_{(0)}<q<q_{(1)}`$, which were taken into account by soft constraints. A few results for the Gardner-Derrida potential (9), with $`\kappa =0`$, $`\alpha =3`$, are displayed on Figs. 2 and 3 for various low temperatures. Note that the point $`(\kappa =0,\alpha =3)`$ lies beyond capacity, while beyond about $`\beta ^1=T=0.2`$ the RS solution satisfies the AT stability condition .

On Fig. 2 the scaled order parameter function $`\beta x(q)`$ is shown. For small temperatures the replica symmetric solution is AT unstable, and we indeed obtain the Parisi form (82) for $`x(q)`$. At $`q_{(0)}`$ near $`0.75`$ the functions $`x(q)`$ jump to zero and remains there as $`q`$ further decreases. The upper plateau with $`x(q)1`$ starts at $`q_{(1)}`$. The consistency of the scaling (82) is confirmed by our finding at $`T=0`$ finite values for $`\beta x(q)`$ if $`q<1`$ and for $`\eta 0.26`$. Interestingly, the curved segment of $`\beta x(q)`$ does not change much with increasing temperature, the main effect being the decrease of $`q_{(1)}`$.

FIG. 2. Scaled order parameter function $`x(q)`$ for $`\kappa =0`$, $`\alpha =3`$ at $`T=0`$ (solid), $`T=0.01`$ (dashed), and $`T=0.1`$ (dotted).

The local stability distribution $`\rho (\mathrm{\Delta })`$ of (22) is displayed on Fig. 3. It was numerically obtained for $`T=0`$ from (92) and for $`T>0`$ from the original formula (72), with the same parameter values as in Fig. 3. For the sake of better visibility of the other details, the very high peaks near $`\mathrm{\Delta }=\kappa =0`$ for $`T=0.1`$ and $`T=0.01`$ as well as the corresponding $`\delta `$-peak for $`T=0`$ have been omitted from this plot. For $`|\mathrm{\Delta }|>3`$ the curves approach zero very quickly. A true gap with $`\rho (\mathrm{\Delta })=0`$ to the left of $`\mathrm{\Delta }=0`$ develops only if $`T=0`$, while for the positive $`T`$-values $`\rho (\mathrm{\Delta })`$ is positive albeit small there. While the gap at $`T=0`$ is present in $`R=0,1,2`$ step replica symmetry breaking schemes, see Refs.  respectively, in all these cases a jump appears near the lower edge. This can be associated with the thermodynamic instability of those saddle points . Our present solution gives linearly vanishing $`\rho (\mathrm{\Delta })`$ at the lower edge, signaling the absence of replicon instability .

FIG. 3. Density of local stabilities $`\rho (\mathrm{\Delta })`$ from theory for $`\kappa =0`$, $`\alpha =3`$ at $`T=0`$ (solid), $`T=0.01`$ (dashed), and $`T=0.1`$ (dotted).

In order to compare theory with practice we performed a medium scale simulation. The standard Hebbian algorithm was modified by Wendemuth  to provide convergence for negative stabilities. Since we chose $`\kappa =1`$, the final steps during stabilization of a pattern went on with $`\mathrm{\Delta }>0`$, so in our case the modification was not essential. The algorithm goes as follows. Firstly random patterns (4) are generated uniformly from an interval centered about zero and normalized as $`_{k=1}^N(S_k^\mu )^2=N`$, and all outputs $`\xi ^\mu `$ are taken uniformly $`1`$. This does not restrict generality since $`S_k^\mu `$ have random signs. The initial coupling vector is set to be proportional to

$$J_k(0)\underset{\mu =1}{\overset{M}{}}S_k^\mu ,$$
(93)
such that its Eucledian norm is $`N`$. In the $`t`$-th step of the algorithm one calculates the local stabilities

$$\mathrm{\Delta }^\mu (t)=\frac{𝐉(t)𝐒^\mu }{|𝐉(t)|}$$
(94)
and selects the least unstable pattern, i. e., the one with the largest $`\mathrm{\Delta }^\mu (t)<\kappa `$. Let us denote its index by $`\mu _0`$, whose argument $`t`$ we omit. Next one augments the couplings as follows. If $`\mathrm{\Delta }^{\mu _{0(t)}}(t)>0`$ then

$$𝐉(t+1)=𝐉(t)+\lambda 𝐒^{\mu _0},$$
(95)
and if $`\mathrm{\Delta }^{\mu _0(t)}(t)<0`$ we have following Wendemuth

$$𝐉(t+1)=𝐉(t)+\lambda \left(𝐒^{\mu _0}+𝐉(t)\frac{N/|𝐉(t)|\mathrm{\Delta }^{\mu _0}(t)}{|𝐉(t)|\mathrm{\Delta }^{\mu _0}(t)}\right).$$
(96)
Here $`\lambda `$ is the gain parameter, the overall scale of increments of $`𝐉`$. In Ref.  the gain parameter was $`\lambda =N^{3/2}`$, after experimentation we chose $`\lambda =N^1`$. Such an increase in the gain parameter did not endanger, rather sped up convergence. At time $`t+1`$ we again look for the least unstable pattern, and so on. The algorithm goes on until it gets stuck with one pattern that we are not able to stabilize in a reasonable time. The intuitive idea behind the algorithm is that since an unstable pattern counts as error irrespective of the distance of the stability parameter $`\mathrm{\Delta }^\mu `$ from $`\kappa `$, one assumes that it is the easiest to stabilize the pattern with $`\mathrm{\Delta }^\mu `$ closest to $`\kappa `$. So one may hope that thus the largest possible number of patterns can be stabilized.

The program ran on $`28`$ PC-s in parallel, each having an AMD K6 processor of $`333`$ MHz, during about one day. Fig. 4 shows both the theoretical curve and the results of the simulation for $`\alpha =\kappa =1`$, a point known to fall beyond the capacity curve (81). The full line is the result of numerical extremization of the free energy functional in the way Fig. 3 was obtained. The Dirac delta peak of the theoretical probability density at $`\kappa `$ is not illustrated. The discontinuous lines represent the histograms for the local stabilities from simulation for two sizes, $`M=N=500`$ and $`1000`$, after normalization.

FIG. 4. Density of local stabilities $`\rho (\mathrm{\Delta })`$ at $`\alpha =\kappa =1`$, axes as in Fig. 3. The theoretical prediction is given by the full line. The two empirical densities are normalized histograms, taken with $`M=N=500`$ and $`1000`$.

The closeness of the two histograms demonstrates that size effects were probably not the cause for the systematic difference between theory and the numerical experiment. A possible ground for the discrepancy is that the algorithm may have been halted prematurely. However, the time necessary for the stabilization of patterns was allowed to grow for each subsequent pattern, and the algorithm was ended only when stabilization did not occur even within the multiple of such an extrapolated time. Another possible reason for the deviation may be that the algorithm got stuck in a “local optimum” without being able to globally maximize the number of stable patterns. In this regard several modified initial conditions were tested but the number of stabilized patterns did not grow in the end. A source of concern can be that the built in random number generator of the C compiler was used; we did not test other routines for this purpose. As to the algorithm, despite its intuitive appeal, there is no proof that it would be able to globally minimize the Hamiltonian (8) with error measure (9). Furthermore, it is likely that with the present parameters the learning task is an NP-complete problem , thus explaining imperfect convergence.

We emphasize that the results represent a significant improvement with respect to the earlier simulation in Ref. . The error per example $`\epsilon `$ found in  is about $`0.21`$, while the present data correspond to $`0.15`$ and theory predicts $`0.1358`$. Thus the deviation between simulation and theory has been decreased by 80%. That means that we stabilized more patterns than , although, given the difference between the theoretical and simulation results, we still could not find the global optimum. Furthermore, an important feature of the density $`\rho (\mathrm{\Delta })`$ is that it should continuously vanish, with a nonzero slope, at the lower edge of the gap. This property is reproduced by the simulation data, in a sharper fashion with the larger $`M=N=1000`$ size, while the value of the edge remains slightly overestimated.

###### Acknowledgements.

The authors acknowledge support by OTKA grant No. T017272 (G. Gy.) and from a special grant for young scientists from the University of Augsburg and by the State of Bavaria within the postgraduate scheme Graduiertenkolleg GRK283 “Nonlinear Problems in Analysis, Geometry, and Physics” (P. R.). It is a pleasure to thank F. Csikor and Z. Fodor for offering us for the simulation their PC farm, supported by grants Nos. OTKA-T22929 and FKFP-0128/1997. Thanks are due to F. Pázmándi for his pointing out and discussing with us the problem of a discontinuity in the potential.


format: text

text:
# DIAGONALIZATION OF A BOSONIC QUADRATIC FORM USING CCM: APPLICATION ON A SYSTEM WITH TWO INTERPENETRATING SQUARE LATTICE ANTIFERROMAGNETS

## 1 Introduction – The Model

It is always possible to diagonalize quadratic bosonic forms (which appear frequently in physics) using a Bogoliubov transformation , but it is a tedious analytical task to find one for a complicated form with many different magnons. Therefore we want to show here how the coupled cluster method (CCM), one of the most powerful and universal techniques in quantum many-body theory (s.  and references therein), can be used in a straightforward scheme to find the exact ground state of such a form.

To be concrete we consider the 2D spin $`1/2`$ Heisenberg model

$$H=J_{AA}\underset{iA_1,jA_2}{}𝐒_i𝐒_j+J_{BB}\underset{iB_1,jB_2}{}𝐒_i𝐒_j+J_{AB}\underset{iA,jB}{}𝐒_i𝐒_j,$$
(1)
which is related to the situation in Ba2Cu3O4Cl2 , a layered quantum antiferromagnet showing significant differences to its parent cuprats like La2CuO4 (see e.g.  for recent experiments). In contrast to La2CuO4 we have two different types of Cu-sites in the Cu-O-planes, namely there are additional Cu(B) atoms located in the centre of every second Cu(A)-O2 square. Within the Cu(A) subsystem we have a strong 180 Cu-O-Cu superexchange yielding to strong antiferromagnetic couplings ($`J_{AA}`$) between Cu(A) atoms, whereas the couplings within the Cu(B) subsystem ($`J_{BB}`$) and between the subsystems ($`J_{AB}`$) are weaker. A recent calculation of $`J_{AA}`$, $`J_{BB}`$ , finding $`J_{AA}10J_{BB}`$ (both antiferromagnetic) agrees with the experimental values . There are also some arguments  for a ferromagnetic $`|J_{AB}|J_{BB}`$.

In the classical ground state (1) shows for $`|J_{AB}|2\sqrt{J_{AA}J_{BB}}`$ a Néel like order for the two subsystems A and B, where the energy is degenerated with respect to the angle $`\phi `$ between the spins of these two subsystems.

## 2 The Method

In this paper we study the ground state properties of (1), using a four-magnon linear spin wave approximation  around the classical ground state, i.e. for each of the four sublattices $`A_1,A_2,B_1,B_2`$ of the two coupled bipartite antiferromagnetic square lattices we introduce different bosonic operators. Thus we get for (1)

$$H=\frac{2N}{3}s^2\left(2J_{AA}+J_{BB}\right)+\underset{𝐤}{}H_𝐤,\text{with}$$
(2)
$$\begin{array}{c}H_𝐤=4J_{AA}s\left(a_{1𝐤}^+a_{1𝐤}+a_{2𝐤}^+a_{2𝐤}\gamma _{𝐤AA}\left[a_{1𝐤}^+a_{2𝐤}^++a_{1𝐤}a_{2𝐤}\right]\right)\hfill \\ +2J_{BB}s\left(b_{1𝐤}^+b_{1𝐤}+b_{2𝐤}^+b_{2𝐤}\gamma _{𝐤BB}\left[b_{1𝐤}^+b_{2𝐤}^++b_{1𝐤}b_{2𝐤}\right]\right)\hfill \\ +J_{AB}s(1+\mathrm{cos}\phi )/2\left(b_{1𝐤}^+a_{1𝐤}+b_{1𝐤}a_{1𝐤}^+b_{2𝐤}^+a_{1𝐤}^+b_{2𝐤}a_{1𝐤}\right)\gamma _{𝐤AB}^1\hfill \\ +J_{AB}s\left(1\mathrm{cos}\phi \right)/2\left(b_{2𝐤}^+a_{2𝐤}+b_{2𝐤}a_{2𝐤}^+b_{1𝐤}^+a_{2𝐤}^+b_{1𝐤}a_{2𝐤}\right)\gamma _{𝐤AB}^2\hfill \\ +J_{AB}s\left(1\mathrm{cos}\phi \right)/2\left(b_{2𝐤}^+a_{1𝐤}+b_{2𝐤}a_{1𝐤}^+b_{1𝐤}^+a_{1𝐤}^+b_{1𝐤}a_{1𝐤}\right)\gamma _{𝐤AB}^1\hfill \\ +J_{AB}s\left(1+\mathrm{cos}\phi \right)/2\left(b_{1𝐤}^+a_{2𝐤}+b_{1𝐤}a_{2𝐤}^+b_{2𝐤}^+a_{2𝐤}^+b_{2𝐤}a_{2𝐤}\right)\gamma _{𝐤AB}^2,\hfill \end{array}$$
(3)
using the lattice structure factors $`\gamma _{𝐤AA}=\mathrm{cos}(k_x/2)\mathrm{cos}(k_y/2)`$, $`\gamma _{𝐤BB}=(\mathrm{cos}k_x+\mathrm{cos}k_y)/2`$ and $`\gamma _{𝐤AB}^{1(2)}=\mathrm{cos}(k_{x(y)}/2)`$.

As stated, we use the coupled cluster method (CCM) to find the exact ground state of (2). To do this we notice the following property of $`H`$

$$\underset{𝐤}{}H_𝐤=\underset{𝐤}{}(H_𝐤+H_𝐤)/2\underset{𝐤}{}H_𝐤^{};[H_𝐤^{},H_𝐤^{}^{}]_{}=0𝐤,𝐤^{}.$$
(4)
Hence it is possible to treat each $`H_𝐤^{}`$ seperately within the CCM, since they all commute with each other. So we have to deal with a bosonic system with eight different bosonic operators $`a_{1\pm 𝐤}`$, $`a_{2\pm 𝐤}`$, $`b_{1\pm 𝐤}`$, $`b_{2\pm 𝐤}`$ denoted with $`a_1,\mathrm{},a_8`$.

The ket and bra ground state of such a system (i.e a many-mode bosonic field theory with bosonic operators $`a_i`$, $`a_i^+`$ in the Hamiltonian) in CCM-SUB$`l`$ approximation is given by

$$\begin{array}{cc}|\mathrm{\Psi }=e^S|0,\hfill & S=_{i_1,i_2,\mathrm{},i_l}A_{i_1,i_2,\mathrm{},i_l}a_{i_1}^+a_{i_2}^+\mathrm{}a_{i_l}^+,\hfill \\ \stackrel{~}{\mathrm{\Psi }}|=0|\stackrel{~}{S}e^S,\hfill & \stackrel{~}{S}=1+_{i_1,i_2,\mathrm{},i_l}\stackrel{~}{A}_{i_1,i_2,\mathrm{},i_l}a_{i_1}a_{i_2}\mathrm{}a_{i_l}\hfill \end{array},$$
(5)
where $`|0`$ is the bosonic vacuum state (i.e. $`a_i|0=0`$), and $`A_{i_1\mathrm{}}`$ and $`\stackrel{~}{A}_{i_1\mathrm{}}`$ are the CCM correlation coefficients. These coefficients are calculated by two systems of equations (one of them is a system of nonlinear equations).

$$\frac{\overline{H}}{\stackrel{~}{A}_{i_1\mathrm{}i_l}}=0,\frac{\overline{H}}{A_{i_1\mathrm{}i_l}}=0,\overline{H}=\stackrel{~}{\mathrm{\Psi }}|H|\mathrm{\Psi },$$
(6)
using the expectation value ($`\overline{H}`$) of the Hamiltonian, i.e. the ground state energy.

Note, that the CCM-SUB2 approximation (i.e. having only quadratic terms of bosonic operators in $`S`$ and $`\stackrel{~}{S}`$ (5)) gives the exact ground state of a quadratic bosonic Hamiltonian, since the ground state wave function of such a Hamiltonian has the form $`|\mathrm{\Psi }=\mathrm{exp}[_{ij}f_{ij}a_i^+a_j^+]|0`$, which can easily be shown using a Bogoliubov transformation (see appendix). Therefore the CCM correlation operator $`S`$ (and $`\stackrel{~}{S}`$ respectively) (5) consist of products of two bosonic creation operators only, all other coefficients $`A_{i_1,\mathrm{},i_l}`$ are zero; so we just have to use SUB2.

To calculate the CCM equations (6) easily using computer algebra, we make use of the Bargmann representation

$$a^+z,a\frac{d}{dz},|01,0|f(a,a^+)|0f(\frac{d}{dz},z)|_{z=0},$$
(7)
which maps the original many-mode bosonic field theory into the corresponding (classical) field theory of complex functions in a particular normed space. So instead of bosonic operators we just have to handle with (complex) numbers and differential operators, which is much easier. Once the (partial nonlinear) equations are obtained they can be solved numerically.

## 3 Results and conclusions

We apply the CCM-scheme described above to calculate the exact ground state of (3) and by doing this getting a spin wave approximate ground state of the model (1). We discuss the energy as a function of the angle between spins of the two subsystems $`A`$ and $`B`$ and the correlation between spins of different subsystems as a function of $`J_{AB}`$ (Fig.1). We find as a typical order from disorder effect, that the degeneracy of the ground state with respect to the angle $`\phi `$ is lifted by quantum fluctuations and a collinear ordering ($`\phi =0,\pi `$) is stabilized. This can clearly be seen by the energy vs. $`\phi `$ picture in Fig.1 and by the correlation $`\text{S}_i\text{S}_j_{A,B}`$ vs. $`J_{AB}`$, which is zero in the classical case, independent of the value of $`J_{AB}`$ (for $`|J_{AB}|2\sqrt{J_{AA}J_{BB}}`$). In the quantum case however that correlation does depend on $`J_{AB}`$, showing again an order effect induced by quantum fluctuations.

In addition we find a lowering of the magnetic order within the subsystems $`A`$ and particular $`B`$ by frustrating $`J_{AB}`$ in the quantum case.

## Acknowledments

This work has been supported by the DFG (Project Nr. Ri 615/7-1).

## Appendix A Proof that CCM-SUB2 gives exact ground state

Using the fact, that a Bogoliubov transformation $`\beta _\nu =_\mu (u_{\mu \nu }^{}a_\mu v_{\mu \nu }^{}a_\mu ^+)`$ exactly diagonalize a quadratic bosonic Hamiltonian with the bosonic operators $`a_i`$, $`a_i^+`$, one can easily show that its ground state must have the form $`|\mathrm{\Psi }=\mathrm{exp}[_{ij}f_{ij}a_i^+a_j^+]|0`$, by showing that $`\beta _\nu |\mathrm{\Psi }=0`$ $`\nu `$. We use the Bargmann representation (7) and get

$$\beta _\nu |\mathrm{\Psi }\stackrel{!}{=}0\underset{\mu }{}\left(u_{\mu \nu }^{}\frac{d}{dz_\mu }v_{\mu \nu }^{}z_\mu \right)\mathrm{exp}\left[\underset{ij}{}f_{ij}z_iz_j\right]\stackrel{!}{=}0z_i$$

$$\underset{\mu }{}(u_{\mu \nu }^{}2\underset{i}{}f_{i\mu }z_iv_{\mu \nu }^{}z_\mu )\stackrel{!}{=}0z_i,2\underset{\mu }{}f_{i\mu }u_{\mu \nu }^{}\stackrel{!}{=}v_{i\nu }^{}$$

and this last matrix equation is allways fulfilled for some $`f_{i\mu }`$.


text:
# 1 Introduction

## 1 Introduction

Horizontal air showers have been studied for many years for several different reasons . Much of the relevance of horizontal air showers induced by cosmic rays is in the understanding of the background against which high energy neutrino showers could be detected. Although no neutrino events can be expected in the Haverah Park data set, even for the most optimistic neutrino flux predictions, it has recently been shown that relevant bounds could be obtained with the Auger Observatories .

Horizontal cosmic ray showers are of great interest in their own right for two principle reasons. Firstly the acceptance of an air-shower array could be doubled if events above $`60^{}`$ can be adequately analysed. Secondly because, uniquely, they deal with surviving particles that are created very close to the shower core, they complement the information obtained in the study of near vertical cosmic ray showers. Moreover understanding the azimuthal asymmetries at large zenith angles \[6,7,8,9,10 (hereafter Paper I)\] can lead to significant improvement of ultra high energy shower analysis at moderate zenith angles ($`30^{}60^{}`$) .

The Haverah Park array, being made of 1.2 m deep water-Čerenkov tanks , is quite possibly the array detector so far constructed which is best suited on geometrical considerations for the analysis of very large inclined showers. Moreover it can be considered as an early prototype of the Auger Observatories which will employ water-Čerenkov tanks of identical depth. The quantitative aspects of our results are very specific to the water-Čerenkov technique.

During the 14 years of operation of the Haverah Park array  which are considered here nearly 10000 air showers were recorded for $`\theta >60^{}`$. The analysis of these data was difficult because of the complex geomagnetic field effects which distort the circular symmetry of air showers, the relatively small size of the array (12 km2) and the limited computing power then available, both for data analysis and simulation. A single remarkable event, perhaps of 1020 eV, at $`85^{}`$ from the zenith which triggered 20 of the array detectors was extensively discussed , but no systematic study was made with those data for which the initial analysis gave a zenith angle $`\theta >60^{}`$.

Here we present the results of the first systematic analysis of horizontal showers from primaries of energy greater than $`10^{17}`$ eV at Haverah Park with zenith angles exceeding $`60^{}`$, comparing them to those expected from the known cosmic ray spectrum. Preliminary results of this work were reported in . For the calculation of the expected rates we make use of a parameterization for the muon number densities at ground level described elsewhere (Paper I). We take into account two possible compositions (proton and iron nuclei) and use alternative models for the high energy interactions namely QGSJET  and SIBYLL 1.6 .

The article is organized as follows: In section 2 we discuss the relative signals and fluctuations expected from the electromagnetic (electrons and photons) and muon components in horizontal air showers (HAS) induced by hadrons and in section 3 we address the muon component of HAS. In section 4 we discuss the detector simulation paying particular attention to the origin of the different parts of the signal from the individual tanks, stressing the differences between vertical and horizontal particle signals. In section 5 we describe in detail the procedure implemented to obtain a prediction for the horizontal shower rate and in section 6 we give the results obtained in the different cases. Our conclusions are presented in section 7.

## 2 Muon and Electromagnetic Components in Horizontal Air Showers

An air shower induced by a proton or a nucleus can be qualitatively understood as a pion shower that is continuously feeding an electromagnetic component (photons, electrons, and positrons) through $`\pi ^0`$ decay, and both a muonic and a neutrino component through charged pion decay. At ground level the typical number density ratio of photons to electrons in a vertical shower is stable. For a low energy threshold of 90 keV this ratio has a value typically between $`20:1`$ and $`30:1`$ for proton or iron as the distance to shower axis $`r`$ changes from 100 m to 1 km. This value is characteristic of electromagnetic showers. The ratio of electrons to muons on the other hand depends strongly on $`r`$. For example for a $`10^{19}`$eV proton shower it drops from $`100:1`$ to $`3:1`$ as $`r`$ rises in the same interval.

In spite of outnumbering the muon component, the average electron and photon energies are typically in the few MeV range, compared to GeV and above for muons. As a result the relative contributions to the signal in a water-Čerenkov tank are not so different because the water-Čerenkov technique is on average much more efficient at converting single muons into signals than electrons or photons.

Direct measurements with Fly’s Eye  have shown that above $`10^{17}`$ eV the shower maximum is usually between 600 g cm-2 and sea level ($`1000`$ g cm$`{}_{}{}^{2})`$. As the zenith angle varies from the vertical, $`\theta =0^{}`$, to the horizontal, $`\theta =90^{}`$, direction, the slant matter depth rises to $`36000`$ g cm-2 and showers are observed well past maximum. The behavior of the electromagnetic and muon components beyond shower maximum is shown in Fig. 1A. While the electromagnetic component of an air shower becomes exponentially attenuated with depth, the muons which do not decay propagate practically unattenuated to the ground, except for energy loss and geomagnetic deviations. As a result the ratio of the electromagnetic to the muon component of an air shower drops as the zenith angle increases up to $`60^{}`$. Above this angle the ratio levels out because the muons themselves produce electromagnetic particles. The remaining electromagnetic component is mainly due to muon decay and to a smaller extent to hadronic interactions, pair production, and bremsstrahlung.

We have studied the effects of muon decay through Monte Carlo simulation using AIRES . Indeed even at $`60^{}`$ the electromagnetic component of a $`10^{19}`$ eV proton shower which can be directly associated to $`\pi ^0`$ decay is already low and confined within a relatively small region of about 200 m around shower axis. At higher zenith angles this component can be neglected and the muon decay contribution becomes stable, roughly at a proportion of 0.8 electrons for each muon, more than three times below the number obtained for vertical showers 1 km away from the shower axis. Unlike the electromagnetic component from pion decay in vertical showers, the lateral distribution follows that of the muons rather closely. This is not difficult to understand as the decay muons give rise to relatively small electromagnetic sub-showers that preserve the muon spatial distribution. The energy distribution of the electromagnetic component is essentially the same as that of an electromagnetic air shower (i.e. for vertical showers) for the same reason.

The contribution of muon interactions to the electromagnetic component can be simply estimated by considering the muon energy spectrum of a single shower and folding it analytically with the bremsstrahlung cross section and the Greisen parameterization, see . For an $`80^{}`$ zenith and $`10^{19}`$eV proton shower the total number of electrons and positrons ($`N_e`$) obtained is about $`2.510^4`$. These are mostly due to the muons in the energy range between 30 GeV and 500 GeV. As this component arises also from electromagnetic sub-showers its energy distribution should also reflect that of electromagnetic cascades. If we multiply, conservatively, the total number of electrons by a factor three to account for the two other muon interactions, it is still a factor of $``$50 below the number of muons in the shower (see figure 1A).

The contribution to the signal from the electromagnetic component at high zenith angles is in the end small and dominated by muon decay. It has been simulated using WTANK  using the muon and electromagnetic signals from different zenith angle showers generated with AIRES  as an input. In Fig. 1B we show the simulated ratio of electromagnetic to muon signal in a Čerenkov tank of 1.2 m depth (as used in Haverah Park and planned for the Auger Observatory) as a function of distance to the shower axis for a vertical shower compared to two inclined showers. The shower particles have been fed through the tank simulation as if they were coming from the vertical direction to eliminate geometric tank effects. The results illustrate the genuine decrease in the electromagnetic to muon signal ratio due only to the shower composition varying with zenith angle. Even at distances from the shower axis of order $`1.5`$ km the ratio is about a factor 3 smaller than for vertical showers.

As will be shown in section 4, in the Haverah Park detectors the relative electromagnetic signal drops even further as the zenith angle increases because of geometric effects. At large zenith angles the muon track length is enhanced compared with that in the vertical direction and the output signal increases accordingly. In the analysis that follows we will consider a track dependent, but otherwise constant correction due to the electromagnetic component from muon decay, in agreement with , and neglect that from all muon interactions in the atmosphere.

## 3 The muon component of horizontal showers

The behavior of the muonic component at high zenith angles becomes extremely complex because of magnetic field effects. It is described in detail in Paper I where an accurate model is presented to account for the average muon number densities at ground level. The inputs for this model are the lateral distribution function (LDF), the average muon energy as a function of radius ($`(r)`$) and the mean distance to the muon production point, all evaluated in the absence of magnetic effects. The model is based on an anticorrelation between the average muon energy and distance of the muon from the core. We have generated these inputs, using AIRES, for two primaries proton and iron, and two hadronic models Quark Gluon String Model (QGSJET) and SIBYLL and for primary energies in the range $`10^{16}`$eV to $`10^{20}`$ eV. The muon energy distribution at a fixed distance to the shower axis is assumed to have a log-normal distribution of width 0.4 (Paper I).

Both the $`(r)`$ and the shape of the LDF are essentially invariant over this large energy interval and only mildly dependent on zenith angle as shown in Fig. 2. As a result the dependence on shower energy of the muon number density distributions with the magnetic field can be parameterized with an absolute normalization to a high level of precision. The dependence of the normalization factor with energy can be obtained by monitoring the total muon number in the showers. The results are plotted in Fig. 3 for four different zenith angles. The normalization or energy scaling can be taken into account by a relation of the following form:

$$N_\mu =N_0E^\beta $$
(1)
where $`N_0`$ and $`\beta `$ are constants for a given model and mass composition as shown in Table 1.

Fluctuations can enhance trigger rates for air showers produced by lower energy primaries because of the steep cosmic ray spectrum . The fluctuations to high number of particles allow the more numerous low energy showers to trigger. We have studied muon number fluctuations at ground level and how they depend on shower development (mean muon production height) and average muon energy. We have found that the mean energy correlates strongly with production height but that most of the number density fluctuations can be accounted for by fluctuations in the total muon number. This is mainly due to fluctuations in the neutral to charged pion ratios in the first interactions. According to the results shown in Table (1) of Paper I it is sufficient to implement a $`20\%`$ RMS fluctuation in the average total muon number as obtained with Eq. 1.

Although the actual distributions obtained with simulation do show a long tail to low muon numbers (up to a factor of 4 reduction), such showers are not expected to be relevant for triggering the array precisely because they have so few muons. As a result it is sufficient to assume a gaussian distribution.

## 4 Detector simulation

### 4.1 The Haverah Park Array

The trigger rate of an air-shower array at large zenith angles is extremely sensitive to the geometry of the array. Factors such as the shape and relative height of detectors become very important for such showers. Figure 4A shows the layout of the Haverah Park array. The relative heights and orientations of the four A-site detectors, the triggering detectors, are shown in figure 4B. A gradient across the array is apparent and this has a significant effect on the observed azimuthal distribution (see section 6). Figure 4C shows the positions of individual tanks within a detector hut. The signals from 15 of the 16 tanks, each of area 2.25 m2, were summed to provide the signal used in the trigger. One tank in each hut was used to provide a low gain signal. See  for a more detailed description of the array.

Water-Čerenkov densities are expressed in terms of the mean signal from a vertical muon (1 vertical equivalent muon or VEM). It has been shown that this signal was equivalent to approximately 14 photoelectrons (pe) .

The formation of a trigger was conditional on: 1. A density of $`>`$0.3 VEM m-2 in the central detector (A1) 2. At least 2 of the 3 remaining A-site detectors recording a signal of $`>`$0.3 VEM m-2.

The rates of the triggering detectors were monitored daily. Over the life of the experiment, after correction for atmospheric pressure effects, the rates of the detectors were stable to better than 5%.

### 4.2 Čerenkov signals of vertical particles

The calculation of the water-Čerenkov signal from horizontal showers is complex. It is informative to consider first the simpler case of vertical showers. A GEANT simulation of the propagation of vertical electrons, gammas, and muons through Haverah Park tanks has been performed (WTANK ). The wavelength dependences of the physical properties of the tank have been taken from  but normalised to the following peak values which are the best estimates for the Haverah Park tanks:

* reflectivity of the walls - 83%
* absorption length of the water - 15 m
* photomultiplier (PMT) quantum efficiency - 22%

The Thorn/EMI, 9618 PMT wavelength acceptance function is taken from the manufacturers specifications.

The results of this simulation are summarized in figure 5. It can be seen that the values used enable us to reproduce the measured signal from an average vertical muon (assuming a mean muon energy of 1 GeV ).

The signal from a vertical muon is composed primarily of the Čerenkov light emitted from the muon track. However there is a significant contribution from Čerenkov light emitted by $`\delta `$-ray electrons (2 pe of the 14 pe total for an average vertical muon).

By contrast the mean energy of electrons and gamma-rays in vertical showers is below 10 MeV. Convolving the energy spectrum of electrons and gammas with the response given by figure 5, we find mean signals of 2.6 and 0.9 photoelectrons for electrons and gamma-rays respectively.

### 4.3 Signals produced by horizontal showers

Three factors complicate the picture for horizontal showers. Firstly for horizontal showers it is possible for Čerenkov photons to fall directly onto the PMT without reflection from the tank walls (we refer to such photons as “direct light”). Secondly the azimuthal asymmetry of the Haverah Park tanks (and detector huts) becomes increasingly important at large zenith angles. Thirdly the mean muon energy grows with zenith angle. For 87 showers the mean muon energy is 200 GeV. For higher energy muons the probability of interaction in the tank is much greater. The production of secondary electrons via pair-production, bremsstrahlung, nuclear interactions (collectively referred to as PBN interactions), and electron knock-on ($`\delta `$-rays) is therefore enhanced. For example the correction due to $`\delta `$-ray production increases from 2 pe at typical vertical muon energies of 1 GeV to around 3 pe for $`>10`$ GeV.

Signal enhancements due to direct light and muon interactions appear in only a fraction of events and hence do not greatly affect the peak of the photoelectron distribution. They do however add long tails to this distribution which are of great importance in calculation of the array trigger rate. Figure 6 shows the different contributions to the water-Čerenkov signal of a horizontal muon as a function of energy. A single Haverah Park tank is considered with muons propagating parallel to the short (1.24 m) side of the tank.

For showers of $`>60^{}`$ the electromagnetic component is almost entirely a product of muon decay as explained in section 2. Convolving the number and energy spectrum of electrons and gammas with the tank response shown in figure 5 we find that the electromagnetic component contributes 2.8 pe for each muon corresponding to 20$`\pm `$2 % of the vertical equivalent muon signal.

## 5 Implementation

We have simulated the rate of the Haverah Park array as a function of the zenith angle following the steps described below. First of all we simulated 100, 1019 eV, showers for zenith angles from $`60^{}`$ to $`88^{}`$ in steps of $`2^{}`$, without magnetic field, for each of the four combinations of primary mass and interaction model as described above. Use is made of the ideas developed in Paper I.

For each zenith angle the following procedure is applied:

1. The inputs of the analytical model (LDF, $`(r)`$ and $`<d>`$) are parameterized.
2. The analytical model is applied to these parameterizations to generate muon density maps of the showers in the transverse plane taking into account magnetic field effects.
3. The density maps for different primary energies are obtained with the energy scaling relationship given by Eq. 1 and the parameters shown in Table 1.
4. Different azimuthal angles are considered by simple rotation algorithms in the approximation explained in section 3 of Paper I.

For each of these zenith angles we have calculated the triggering probability of the Haverah Park array in 40 energy bins between $`10^{16}`$ eV and $`10^{20}`$eV and 18 azimuth bins as follows:

1. The shower is directed on to the array N times with random core locations on the ground plane up to X km from the centre of the array. N and X are varied according to the primary energy and zenith we are dealing with. Typical values are $`N=10000`$ and $`X=8`$ km.
2. Each time a shower is directed at the array, the total muon number ($`N_\mu `$) is fluctuated with a gaussian distribution of spread 0.2$`N_\mu `$ to take into account shower fluctuations.
3. The density in the ground plane at the location of each of the trigger detectors is read from the muon density maps with an appropriate projection, taking into account the corrections due to the different heights of the detectors and the effects of the magnetic field as explained in section 7 of Paper I.
4. The corresponding signal in each of the trigger detectors is generated (see next subsection).
5. The trigger condition of the Haverah Park array is tested.
6. If an event is deemed to trigger the array then the primary energy and the combination of detectors contributing to the trigger are recorded.

The triggering probability derived in this way has been convolved with a primary energy spectrum , derived from Akeno and Haverah Park data, to obtain the event rate at each zenith and azimuth angle. The final step before comparison with data is to convolve the obtained zenith angle distribution with the appropriate measurement errors. The arrival directions of all the data that had zenith angle greater than $`60^{}`$ have been reanalysed assuming a plane front for the muons and using all available timing information. Previously only times from the four central triggering detectors were used to compute the arrival directions . Here up to 16 times were used in the analysis with the additional data being from the detectors $``$2 km from A1 in the ground plane. The uncertainly in arrival direction for each event has accordingly been significantly reduced. Fig. 7 shows the average uncertainty in zenith angle of the reconstructed events as a function of zenith angle. The uncertainty in azimuth angle is approximately constant and is $`1^{}`$. The events considered here have not been analysed for core position. Accordingly some of the very large zenith angle events may have fallen a long way outside the boundary of the array and consequently have a larger uncertainty in zenith angle than suggested by fig. 7.

We have rejected about $`3\%`$ of the events having an error exceeding $`4.5^{}`$ in the reconstructed zenith angle $`\theta `$. The RMS error in zenith angle shown in fig. 7 is used to *smear* the calculated distribution.

### 5.1 Implementation of the signal in detector

The signal is generated as follows:

1. The projected area of the detector in the shower plane is calculated.
2. Given the local muon density and the projected area, we sample the number of incident muons from a Poisson distribution.
3. The track length of each muon through the detector is sampled from a distribution obtained analytically from the detector geometry (see figure 4C).111This distribution accounts for the possibility that a single muon may traverse several tanks.
4. The contribution of indirect Čerenkov light from the incident muons and from $`\delta `$-ray electrons is calculated from the sampled track lengths (12 pe for each 1.2 m of track, with an additional 3 pe/1.2 m to account for the signal from $`\delta `$-rays as described in the previous section).
5. The signal from direct light on the PMTs is related to the detector geometry in a more complex way and is implemented using WTANK to simulate the passage of muons through the whole detector for a range of zenith and azimuth angles.
6. The probability of PBN interactions in the detector is calculated from the track length and the muon energy (obtained from the $`(r)`$ with a zenith angle dependent correction for the magnetic field).
7. If an interaction is deemed to occur then the energy of the resulting em-cascade is sampled from the appropriate differential cross-section. The signal produced by this cascade is calculated using an analytical approach.
8. The electromagnetic component of the shower is approximated by the addition of 2.8 pe for each muon as discussed in the previous section.

The contributions to the signal depend in different ways on the zenith and azimuth. The mean signal for those contributions proportional to the track are constant with azimuth and zenith, because of the compensating effect of the projected area of the detector. In fact the mean signal is simply proportional to the volume of the detector. For larger zenith angles fewer muons produce the same signal so Poisson fluctuations become more important. The mean signal produced by direct light increases with zenith and has a complex dependence on azimuth. The mean signal from the electromagnetic component is proportional to the projected area, and hence decreases with zenith and has a simpler azimuthal dependence.

We note that the use of distributions (rather than mean values) for different contributions to the signal is essential for accurate calculation of the event rate.

## 6 Results and comparison with data

### 6.1 Energy distribution of showers that trigger the array

In Fig. 8 we show the energy distribution of showers which trigger the array as calculated with the QGSJET and SIBYLL models and for proton and iron primaries at two zenith angles. The mass dependence seen in Fig. 8 is a consequence of the greater number of muons produced by heavier primaries. The choice of interaction model has a some what smaller effect on the energy response.

Fig. 9 shows the mean and width of the energy distribution as a function of zenith angle (for QGSJET with proton primaries). The median energy changes from $`6\times 10^{17}`$ eV at 70 to $`5\times 10^{18}`$ eV at 84.

### 6.2 Azimuthal distribution

One of the most stringent tests for the simulation of horizontal air showers is the prediction of the azimuthal distributions. Minor corrections in the altitude of the detectors, or the orientation of the tanks can lead to large differences in the rates for inclined showers. The total event rate is dominated by events triggering only 3 of the 4 A-site detectors. The azimuthal distribution is best understood by separating these 3-fold events into 3 sets, each set requiring that a particular detector *does not* trigger. The recorded azimuthal distributions for each of these three sets (A2 out, A3 out and A4 out) are shown in Fig. 10. The calculated distributions (using QGSJET with primary protons) are shown normalised to the data for comparison of the shape. The agreement is reasonable with the exception of the set “A2 out”. There are many possible reasons for this disagreement. In particular we will in future consider reflection and absorption in the ground around the detectors, but we are encouraged by the overall result.

### 6.3 Zenith angle distribution

In Fig. 11 we show the contributions to the total rate (from each of the processes considered) as a function of zenith angle. The QGSJET model and proton primaries are assumed. Similar plots for other hadronic models and primaries differ mainly in their normalization. It should be noted that inclusion of only the muon track contribution would lead to underestimating the total rate by a factor $`3`$ ($`4`$) at 65 (85). The effects of delta rays, direct light, electromagnetic corrections, and pair production and bremsstrahlung are all very significant in the calculation of the overall rate. The effect of fluctuations in the total muon number in a shower represents at most a $`20\%`$ rate enhancement because of the large width of the energy response at a given angle.

The left hand plot of Fig. 12 shows the number of events detected at Haverah Park as a function of the zenith angle compared to the expected number of events for two different models and two mass compositions. The error bars on the data points are statistical. The right hand plot illustrates the effect of the uncertainty in the primary energy spectrum on our result. The maximum and minimum predictions from the spectrum described in  (using one sigma errors) are shown together with the result obtained using the primary spectrum parameterized in  (all for proton primaries using the QGSJET model).

We stress that no normalization was done with the simulated events. The shape of the simulated distribution is in good agreement with data between $`60^{}`$ and $`70^{}`$ while at higher zenith angles there is a mild disagreement, at most a $``$30% effect. Within the uncertainties resulting from the zenith angle measurement, the primary mass, the interaction models and the primary energy spectrum, the simulated overall rate of showers above $`60^{}`$ is completely consistent with the data.

For most assumptions of primary spectrum and interaction model an intermediate primary mass in the decade around $`10^{18}`$ eV seems to be required to fit the data. Further work is in progress to select higher energy events by considering detectors (B – H of Fig. 4A) that are not in the trigger. For example the event discussed in  and believed to be of energy $`10^{20}`$eV struck 20 of the detectors.

## 7 Conclusion

We have presented the first comparison of $`>10^{17}`$ eV air-shower data above $`60^{}`$ with simulation. The agreement is reasonable and allows us to state that the detection of HAS induced by hadronic primaries in an array of water-Čerenkov tanks can be explained in terms of the muonic component. The effects of all the considered contributions to the tank signals have been shown to be extremely important in understanding the total rate.

Our final results differ from the preliminary ones reported earlier  in several respects. We have taken into consideration the electromagnetic component from muon decay separately. Previously this contribution was taken together with the correction due to $`\delta `$rays. We have further improved the implementation of the PBN corrections taking into account their increase with the muon energy. We have also simulated the ground number density profiles in a completely different fashion making use of the model described in Paper I.

The general shape agreement between data and simulation both for the zenith angle and azimuthal angle differential rates, together with agreement in the absolute normalization makes us confident that the geometrical considerations have been taken into account correctly, at least at the 20% level. This, together with the general agreement with our preliminary rate simulation, which used Monte Carlo shower simulation rather than the analytical approach described in Paper I, provides strong supporting evidence for this approach as an important tool for HAS studies.

The normalization of the shower rate is sensitive to the model used for high energy hadronic interactions. For proton showers the absolute normalization of the rate differs by up to a factor of two for the two models considered. The absolute rate has also been shown to be very sensitive to the primary composition. There is about a factor 2.2 (3) enhancement, reasonably independent of zenith angle, for the differential rate normalization if the primary particles are assumed to be iron nuclei in the QGSJET (SIBYLL) model. It is interesting to note that current air-shower models are not able to reproduce the observations of high energy particles close to the shower axis above 1016 eV . This may be a possible explanation of the $``$30% disagreement above 80 degrees (where the surviving muons all originate close to the shower core).

These effects illustrate that HAS measurements may provide another tool with which to study both primary composition and hadronic interactions in ultra high energy cosmic rays. Future detectors such as the Auger Observatories should benefit significantly from the study of cosmic ray showers at large zenith angles.

Acknowledgements: We thank Xavier Bertou for helping us with the angular reanalysis of the data described in  and Gonzalo Parente for suggestions after reading the manuscript. This work was partly supported by a joint grant from the British Council and the Spanish Ministry of Education (HB1997-0175), by Xunta de Galicia (XUGA-20604A98), by CICYT (AEN99-0589-C02-02) and by PPARC(GR/L40892).