text:
# Novel finite-differencing techniques for numerical relativity: application to black hole excision

## Abstract

We use rigorous techniques from numerical analysis of hyperbolic equations in bounded domains to construct stable finite-difference schemes for Numerical Relativity, in particular for their use in black hole excision. As an application, we present 3D simulations of a scalar field propagating in a Schwarzschild black hole background.

The numerical implementation of Einstein’s equations represents a daunting task. The involved nature of the equations themselves, and a number of technical difficulties (related to the necessarily finite computational domain, the limitations in relative computational power and the need to deal with singularities, constraints and gauge freedom) imply a significant challenge.

Numerical solutions of Einstein’s equations involve solving a nonlinear set of partial differential equations on a bounded domain, and formally constitute an initial-boundary-value problem (IBVP). Constructing stable and long term well behaved numerical approximations for such systems with boundaries is highly non-trivial. Here the term numerical stability is used in the sense that is equivalent, through Lax’s theorem, to convergence (that is, the property that the numerical solution will approach the continuum one when resolution is increased). Delicate complications arise due to corner and edges at outer (and possible inner) boundaries, all of which introduce subtleties for a stable implementation.

Recently, however, several sophisticated tools of rigorous numerical analysis have been developed for systematically constructing stable numerical schemes for IBVP’s. They employ a discrete form of well-posedness and thus are stable by construction, at least for linear systems. At this time, these tools have essentially not been used by the numerical relativity community. The purpose of this paper is to outline their use, in particular for black hole excision.

An IBVP consists of three ingredients: a partial differential equation, initial and boundary data. It is well-posed if a solution exists, is unique, and depends continuously on the initial and boundary data. It is well known that stable finite difference schemes approximating an IBVP can only be constructed for well-posed systems. While problems in general relativity are typically nonlinear we consider here the linear IBVP, as stability in the linearized case is a necessary condition for stability in the full nonlinear system and these methods may also be applied to nonlinear equations. Consider the linear IBVP on a domain $`[0,\mathrm{})\times \mathrm{\Omega }`$

$`_tu`$ $`=`$ $`A(t,\stackrel{}{x})^i_iu+B(t,\stackrel{}{x})u,`$ (1)
$`u(0,\stackrel{}{x})`$ $`=`$ $`f(\stackrel{}{x}),`$ (2)
$`w_+(t,\stackrel{}{x})`$ $`=`$ $`Sw_{}(t,\stackrel{}{x})+g(t,\stackrel{}{x}),\stackrel{}{x}\mathrm{\Omega },`$ (3)
where $`u`$ is a vector-valued function, $`w_+`$ and $`w_{}`$ are incoming and outgoing modes, and $`S`$ is sufficiently small (maximally dissipative boundary conditions maxdissbc ). The system is assumed to be symmetric hyperbolic, i.e., there exists a symmetric positive definite matrix $`H(t,\stackrel{}{x})`$, the symmetrizer, satisfying $`HA^i=(HA^i)^T`$. The usual proof of well-posedness proceeds by defining an “energy”, $`=_\mathrm{\Omega }u^THud^3x`$, and by showing that $``$ can be bounded as a function of the initial and boundary data. Analogously, a way to construct stable numerical schemes is to design them such that energy estimates hold at the discrete level gko . This is the procedure we shall follow in this paper.

Our construction involves four steps. (1) We construct discrete derivative operators and a scalar product so that a semi-discrete energy estimate holds. The spatial derivatives are then discretized using these operators, resulting in a semi-discrete system of ordinary differential equations. (2) We impose boundary conditions in a way that does not spoil the previous semi-discrete energy estimate. In particular, we account for effects of boundary edges and vertices. Steps (1) and (2) guarantee numerical stability of the semi-discrete system. (3) We may add artificial dissipation and/or arrange the spatial discretization to achieve optimal energy bounds. (4) Finally, an appropriate time integrator is chosen so that fully discrete stability holds.

Step (1): The equations are discretized on a domain $`\mathrm{\Omega }`$ with an inner boundary to accommodate for black hole excision. While several grid geometries are possible, we concentrate on the simplest case of a cubical domain from which an inner, smaller cube has been removed (both cubes aligned with the grid). We introduce the grid points $`\stackrel{}{r}_{ijk}=(i\mathrm{\Delta }x,j\mathrm{\Delta }y,k\mathrm{\Delta }z)\mathrm{\Omega }`$ and assume that some of these points lie on the boundary $`\mathrm{\Omega }`$. A scalar product, $`\mathrm{\Sigma }`$, between any two real vector-valued grid functions $`u`$ and $`v`$, is defined as

$$(u,v)_\mathrm{\Sigma }=\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z\underset{i,j,k=0}{\overset{N}{}}\sigma _{ijk}u_{ijk}^Tv_{ijk},$$
(4)
where $`\sigma _{ijk}`$ are defined below and in the sum $`\stackrel{}{r}_{ijk}\mathrm{\Omega }`$. A semi-discrete energy is then defined by $`E=(u,Hu)_\mathrm{\Sigma }`$.

A key ingredient for deriving the continuum energy estimate is integration by parts. Similarly, in order to get a semi-discrete energy estimate, one must construct the difference operators, $`D`$, and $`\mathrm{\Sigma }`$ in such a way that the discrete version of integration by parts, called summation by parts (SBP), holds. In one dimension, SBP is expressed by $`(u,Dv)_\mathrm{\Sigma }+(Du,v)_\mathrm{\Sigma }=u_Nv_Nu_0v_0`$, and this definition can be generalized to higher dimensions and more complicated domains.

A detailed calculation shows that SBP for our chosen $`\mathrm{\Omega }`$ holds by defining (similar definitions for $`y,z`$ directions)

$$D^{(x)}=\{\begin{array}{cc}D_\pm ^{(x)}u_{ijk}\hfill & \text{at }x=\text{const.}\text{ faces,}\hfill \\ & x=\text{const.}\text{ OB edges}\hfill \\ & \text{and OB vertices}\hfill \\ \left(\frac{1}{3}D_\pm ^{(x)}+\frac{2}{3}D_0^{(x)}\right)u_{ijk}\hfill & \text{at }x=\text{const.}\text{ IB edges}\hfill \\ \left(\frac{1}{7}D_\pm ^{(x)}+\frac{6}{7}D_0^{(x)}\right)u_{ijk}\hfill & \text{at IB vertices}\hfill \\ D_0^{(x)}u_{ijk}\hfill & \text{everywhere else}\hfill \end{array}$$

where IB and OB stand for inner and outer boundary respectively, $`D_+^{(x)}u_{ijk}=(u_{i+1jk}u_{ijk})/\mathrm{\Delta }x`$, $`D_{}^{(x)}u_{ijk}=(u_{ijk}u_{i1jk})/\mathrm{\Delta }x`$ and $`D_0^{(x)}u_{ijk}=(u_{i+1jk}u_{i1jk})/(2\mathrm{\Delta }x)`$. The $`\pm `$ signs indicate whether non-excised points are to the right or left, respectively, of the boundary point. The weights $`\sigma _{ijk}`$ in (4) are defined to be $`1`$ in the interior, $`1/2`$ at the faces of IB and OB, $`1/4`$ at the edges of OB, $`3/4`$ at the edges of IB, $`1/8`$ at the vertices of OB, $`7/8`$ at the vertices of IB, and zero in the excised region. This difference operator is second-order accurate at the interior and first-order at boundaries, yielding in principle, a scheme with overall second-order accuracygus .

Step (2): In contrast to the continuum problem, an examination of boundary terms left after SBP in the energy estimate indicates that edge and vertex boundary points make a finite contribution to $`E`$ at fixed resolution. These contributions show precisely how to impose boundary conditions at these points. They have to be imposed on incoming characteristic modes defined by some “effective” normal vectors $`n`$. For instance, on a uniform grid ($`\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }z`$), at an edge one has $`n=(n_a+n_b)/\sqrt{2}`$, and at a vertex $`n=(n_a+n_b+n_c)/\sqrt{3}`$ (where $`n_a`$, $`n_b`$ and $`n_c`$ are unit vectors of the intersecting faces). Boundary conditions along these effective directions have to be imposed without destroying the semi-discrete energy estimate. Following Olsson olsson , we do so by projecting the right hand side (RHS) of the evolution equations at boundary points onto the subspace of the grid functions that satisfy the discretized boundary conditions. The projection is by construction self-adjoint with respect to $`\mathrm{\Sigma }`$, which ensures that the solution of the projected semi-discrete system will satisfy the boundary conditions without compromising numerical stability.

Step (3): At a fixed resolution, even convergent codes may generate significant error growth as time progresses, and it is often desirable to minimize this growth. This may be achieved by adding artificial dissipation, by rearranging the semi-discrete equations in a specific manner to obtain optimal estimates, or by doing both. We briefly summarize these two procedures.

In numerical simulations a dissipative term, $`(Q_d^{(x)}+Q_d^{(y)}+Q_d^{(z)})u`$, is sometimes added to the RHS of Eq. (1), which damps high frequency modes but does not affect the accuracy of the scheme. Requiring $`(u,Q_d^{(x)}u)_\mathrm{\Sigma }0`$ ensures that the semi-discrete energy estimate still holds. We have modified the Kreiss–Oliger dissipation operator for our black hole excision geometry and $`\mathrm{\Sigma }`$. This modified operator has the standard form on the grid interior (shown here for the $`x`$ direction)

$`Q_d^{(x)}u_{ijk}`$ $`=`$ $`ϵ\mathrm{\Delta }x^3(D_+^{(x)}D_{}^{(x)})^2u_{ijk},`$
(where $`ϵ`$ is a parameter satisfying $`ϵ0`$.) However, it is modified near boundaries to satisfy $`(u,Q_d^{(x)}u)_\mathrm{\Sigma }0`$. Let $`(i_0,j,k)`$ label a point on a $`x=\text{constant}`$ boundary with neighbors $`(i_01,j,k)`$ and $`(i_0+1,j,k)`$. Near this boundary point $`Q_d^{(x)}`$ becomes

$`Q_d^{(x)}u_{i_01}=ϵ\mathrm{\Delta }x(D_{}^{(x)2}2D_+^{(x)}D_{}^{(x)})u_{i_01},`$
$`Q_d^{(x)}u_{i_0}={\displaystyle \frac{ϵ\mathrm{\Delta }x}{\sigma _{i_0}}}(\sigma _{i_01}D_{}^{(x)2}+\sigma _{i_0+1}D_+^{(x)2})u_{i_0},`$
$`Q_d^{(x)}u_{i_0+1}=ϵ\mathrm{\Delta }x(D_+^{(x)2}2D_+^{(x)}D_{}^{(x)})u_{i_0+1},`$
where the indices $`jk`$ are suppressed for clarity, and $`Q_d^{(x)}`$ is set to zero for points outside the domain.

A second method for controlling unnecessary growth is by rearranging the discretized equations so the optimal estimates for the continuum hold at the semi-discrete level. For instance, if at the continuum the system exhibits energy conservation, it is important to preserve it at the semi-discrete level. For example, assuming that Eq. (1) has time-independent coefficients and that $`_i(HA^i)=HB+(HB)^T`$, $``$ is conserved if appropriate boundary conditions are given. The solution of $`_tu=A^i(x)D_iu+B(x)u`$ may yield growth in $`E`$, due to the lack of a Leibniz rule at the discrete level. However, one can show that rearranging the semi-discrete RHS of (1) as

$$\frac{1}{2}A^iD_iu+\frac{1}{2}H^1D_i(HA^iu)+\left(B\frac{1}{2}H^1_i(HA^i)\right)u,$$

gives a non-increasing $`E`$, and therefore the discrete solution will not grow. This idea is closely related to the concept of strict stability olsson .

Step (4): Finally, stability of the fully discrete system follows if it is integrated with a scheme that satisfies the (necessary and sufficient) local stability condition or the (sufficient) preservation of energy estimate’s condition (e.g., third or fourth order Runge–Kutta rk ).

An example: Wave propagation on a black hole spacetime $`g_{\mu \nu }`$ presents some of the same challenges as the full Einstein equations (e.g. constraint preservation, boundary conditions, excision of the singularity). Thus it is an ideal test bed to demonstrate the techniques described in this paper. We write the wave equation in first-order form,

$$_\mu \mathrm{\Phi }=V_\mu ,^\mu V_\mu =0,_\mu V_\nu =_\nu V_\mu ,$$
(5)
with $`_\mu `$ the covariant derivative. To split these equations into evolution equations and constraints, we specify a future-directed time-like vector $`u^\mu `$, and contract the first and the last equation with it. The evolution equations are

$$\mathrm{\pounds }_u\mathrm{\Phi }=u^\mu V_\mu \mathrm{\Pi },^\mu V_\mu =0,\mathrm{\pounds }_uV_\mu =_\mu \mathrm{\Pi },$$
(6)
where $`\mathrm{\pounds }_u`$ is the Lie derivative with respect to $`u^\mu `$. This system is symmetric hyperbolic. Due to the evolution equations, the constraints $`C_\mu V_\mu _\mu \mathrm{\Phi }=0`$, $`C_{\mu \nu }_\mu V_\nu _\nu V_\mu =0`$ are Lie-dragged by $`u^\mu `$, i.e., $`\mathrm{\pounds }_uC_\mu =0`$, $`\mathrm{\pounds }_uC_{\mu \nu }=0`$.

Although the discussion below applies to any background geometry, for definiteness we specialize to the Schwarzschild background. This spacetime possesses a time-like Killing field, $`k^\mu =(_t)^\mu `$, and the future directed unit normal $`n^\mu `$ to the $`t=\text{const.}`$ slices is given by $`(k^\mu \beta ^\mu )/\alpha `$. We denote with $`\alpha ,\beta ^i`$ and $`h_{ij}`$ the lapse, shift and three-metric, respectively. A natural choice for the vector field $`u^\mu `$ is $`n^\mu `$, however this requires care specifying boundary data compatible with the constraints. For a particular coupling between the in- and out-going modes, we have found maximally dissipative boundary conditions compatible with the constraints. Unfortunately, they imply reflective boundary conditions in the sense that (for homogeneous boundary data and in the absence of inner boundaries) the physical energy of the scalar field $`\mathrm{\Phi }`$ is exactly conserved. Choosing a different formulation with $`u^\mu =k^\mu `$ allows us to impose radiative boundary conditions, in the sense that the physical energy decreases when the wave reaches the boundary. In this formulation the constraint variables $`C_\mu `$ and $`C_{\mu \nu }`$ propagate tangentially to the boundaries, and thus the constraints are automatically satisfied when satisfied initially. Moreover, this formulation has the advantage that its symmetrizer agrees with the physical energy. Unfortunately, as the Killing field $`k^\mu `$ becomes space-like inside the black hole, the symmetrizer is not positive definite in this region.

We combine the advantages of both formulations by choosing $`u^\mu =(k^\mu b^\mu )/\alpha `$, where $`b^\mu `$ is a smooth vector field which is tangential to the $`t=\text{const.}`$ slices, vanishes in a neighborhood of the outer boundary, and agrees with the shift vector $`\beta ^i`$ as one approaches the event horizon from the outside region. This new interpolating formulation is manifestly symmetric hyperbolic also inside the black hole; and at the outer boundary radiative boundary conditions can be given. Equations (6) yield

$`_t\mathrm{\Pi }`$ $`=`$ $`b^i_i\mathrm{\Pi }+{\displaystyle \frac{\mathrm{\Delta }^i}{\alpha }}_i(\alpha \mathrm{\Pi })+{\displaystyle \frac{1}{\sqrt{h}}}_i(\sqrt{h}\mathrm{\Delta }^i\mathrm{\Pi })`$ (7)
$`+`$ $`{\displaystyle \frac{1}{\sqrt{h}}}_i\left(\alpha \sqrt{h}H^{ij}V_j\right){\displaystyle \frac{1}{\sqrt{h}}}_j\left(\sqrt{h}{\displaystyle \frac{b^j}{\alpha }}\right)\mathrm{\Delta }^iV_i`$
$`+`$ $`{\displaystyle \frac{1}{\alpha }}\left(\beta ^j_jb^ib^j_j\beta ^i\right)V_i+{\displaystyle \frac{1}{\sqrt{h}}}_i(\sqrt{h}b^i)\mathrm{\Pi }`$
$`_tV_i`$ $`=`$ $`_i(\alpha \mathrm{\Pi })+b^j_jV_i+V_j_ib^j`$ (8)
where $`\mathrm{\Delta }^i=\beta ^ib^i`$, $`H^{ij}=h^{ij}\mathrm{\Delta }^i\mathrm{\Delta }^j/\alpha ^2`$, and $`h^{ij}`$ denotes the inverse three metric. This system is symmetric hyperbolic with respect to the energy $`=\frac{1}{2}_\mathrm{\Omega }\alpha \left(\mathrm{\Pi }^2+H^{ij}V_iV_j\right)\sqrt{h}d^3x`$. The change of $``$ under the flow generated by (7,8) can be bounded as

$$\frac{d}{dt}\underset{OB}{}\frac{1}{4}\left(w_+^2w_{}^2\right)\sqrt{h}d^2x+_\mathrm{\Omega }Fd^3x,$$
(9)
where $`w_\pm =\mu _\pm \mathrm{\Pi }+\alpha H^{ij}n_iV_j`$ are the characteristic variables which have nonzero speeds, $`\mu _\pm =\pm \alpha +\beta ^in_i`$, and $`n_i`$ is the normal to the outer boundary, normalized such that $`h^{ij}n_in_j=1`$. The expression $`F`$ in Eq. (9) vanishes if $`b^i=0`$. Then setting $`w_+`$ to zero guarantees that $``$ will decrease when the wave reaches the boundary. In the previous estimate we have assumed that the excised cube is sufficiently small (side length smaller than $`4\sqrt{3}M/9`$) so that the inner boundary is purely outflow.

Equations (7)–(8) have been written such that the replacement $`_i`$ by $`D_i`$ (only on the dynamical variables) automatically leads to a semi-discrete system with energy conservation in the case $`b^i=0`$. The semidiscrete solution will not grow, even at fixed resolution.

The figures show the results of the numerical implementation of the techniques for the scalar wave example we discussed; we have chosen Kerr-Schild coordinates (we have checked that the case of Painlevé-Gullstrand coordinates behaves similarly and so does the Maxwell field on these backgrounds) . Figure 1 shows evolutions with the computational domain outside the event horizon; fig. 2 shows results of a run with singularity excision. Movies can be see in movie .

In this article we have shown that the use of rigorous numerical analysis techniques can be used in practice to obtain stable schemes for numerical relativity. This maps out a precise road for robust implementations, reducing the need for expensive trial and error standard methods. The successful application of these techniques to the simulation of fields propagating on black hole backgrounds, which requires dealing with singularity excision, corners and edges, illustrates the potential of the approach. This opens the door for applying the same techniques to the full equations of general relativity.

Acknowledgments: We thank E. Tadmor, P. Olsson, and H. O. Kreiss for helpful discussions. This work was supported in part by NSF grant PHY9800973, the Horace Hearne Jr. Institute for Theoretical Physics, and Fundación Antorchas. Computations were done at LSU’s Center for Applied Information Technology and Learning and parallelized with the CACTUS toolkit cactus . GC, LL, OR, OS and MT thank the Caltech Visitors Program in numerical relativity for hospitality.


format: text

text:
# Optimized Neural Networks to Search for Higgs Boson Production at the Tevatron.

## 1 The basic idea

In High Energy physics a discrimination between a signal and its corresponding backgrounds by Neural Networks (NN) is especially remarkable when the data statistics are limited. In this case it is important to optimize all steps of the analysis. One of the main questions which arises in the use of NNs is which, and how many variables should be chosen for network training in order to extract a signal from the backgrounds in an optimal way. The general problem is rather complicated and finding a solution depends on having a concrete process for making the choice, because usually it takes a lot of time to compare results from different sets of variables.

One observation which helps in making the best choice of the most sensitive variables is to study the singularities in Feynman diagrams of the processes. Let us call those kinematic variables in which singularities occur as ”singular variables”. What is important to stress here is that most of the rates for both the signal and for the backgrounds come from the integration over the phase space region close to these singularities. One can compare the lists of singular variables and the positions of the corresponding singularities in Feynman diagrams for the signal process and for the backgrounds. It is obvious that if some of the singular variables are different or the positions of the singularities are different for the same variable for the signal and for the backgrounds the corresponding distributions will differ most strongly. Therefore, if one uses all such singular variables in the analysis, then the largest part of the phase space where the signal and backgrounds differ most will be taken into account. One might think that it is not a simple task to list all the singular variables when the phase space is very complex, for instance, for reactions with many particles involved. However, in general, all singular variables can be of only two types, either s-channel: $`M_{f1,f2}^2=(p_{f1}+p_{f2})^2`$, where $`p_{f1}`$ and $`p_{f2}`$ are the four momenta of the final particles $`f1`$ and $`f2`$ or t-channel: $`\widehat{t}_{i,f}=(p_fp_i)^2`$, where $`p_f`$ and $`p_i`$ are the momenta of the final particle (or cluster) and the initial parton. For the $`\widehat{t}_{i,f}`$ all the needed variables can be easily found in massless case: $`\widehat{t}_{i,f}=\sqrt{\widehat{s}}e^Yp_T^fe^{|y_f|}`$, where $`\widehat{s}`$ is the total invariant mass of the produced system, and Y is the rapidity of the total system (rapidity of the center mass of the colliding partons), $`p_T^f`$ and $`y_f`$ are transverse momenta and pseudorapidity of the final particle f. The idea of using singular variables as the most discriminative ones is described in  and the corresponding method was demonstrated in practice in .

Singular variables correspond to the structure of the denominators of Feynman diagrams. Another type of interesting variables corresponds to the numerators of Feynman diagrams and reflects the spin effects and the corresponding difference in angular distributions of the final particles. In order to discriminate between a signal and the backgrounds, one should choose in addition to singular variables mentioned above those angular variables whose distributions are different for the signal and backgrounds. The set of these singular and angular variables will be the most efficient set for a NN analysis.

The third type of useful variables which we called ”Threshold” variables are related to the fact that various signal and background processes may have very different thresholds. Therefore the distributions over such kind of variables also could be very different keeping in mind that effective parton luminosities depend strongly on $`\widehat{s}`$. The variable $`\widehat{s}`$ would be a very efficient variable of that kind. However, the problem is that in case of neutrinos in the final state one can not measure $`\widehat{s}`$ and should use the effective $`\widehat{s}`$ which is reconstructed by solving t-,W-mass equations for the neutrino longitudinal momenta. That is why we propose to use not only the effective variable $`\widehat{s}`$ but the variable $`H_T^{jets}`$ as well.

To apply the method it is important to use a proper Monte-Carlo model of signal and background events which includes all needed spin correlations between production and decays. For the following analysis we have calculated the complete tree level matrix elements for the background processes with all decays and correlations by means of the CompHEP program . The corresponding events are available at the FNAL Monte-Carlo events database .

## 2 Applying the method

The present estimation of the expected sensitivities for the light Higgs boson search at the Tevatron by means of NNs is given in . Based on the method described above we improve the efficiency of the NN technique. In the analysis we choose the Higgs boson mass to be $`M_H=115`$ GeV. We model the detector smearing by the SHW package .

First of all we exclude ineffective variables from the old set , like $`P_T^e`$ from the $`W`$-boson (shown at the left plot in Fig. 1). After the corresponding analysis of Feynman diagrams and comparison of kinematical distributions we added the new variables for NN training. The example distribution for the new variable ($`cos(z_{axis},e)`$) is shown in the right plot of Fig. 1.

At the next step we constructed the set of NNs for pairs of the signal (WH) and each of the background from the complete set of principle backgrounds ($`Wb\overline{b},WZ,t\overline{t},tb(j)`$).

The standard steps of NN training were used for the NNs with the old set of input variables and with the new one. Efficiencies of networks with different sets have been compared based on the criteria that for the better net the “Error function” $`E=\frac{1}{N}_{i=1}^N(d_io_i)^2`$, where $`d_i`$ and $`o_i`$ are the desired and real outputs of the net and $`N`$ is the number of test events, is smaller. Two examples of distributions you can see in Fig 2 for the $`WHt\overline{t}`$ network (left plot) and $`WHWZ`$ network (right plot). One can see a significant improvement for the networks with new input sets in comparison with old sets of variables, since the corresponding curves of the error function are significantly lower.

## 3 Results

Based on the described method we have constructed the new NNs to search for a light Higgs boson at the Tevatron. After checking the improvement in efficiency of new networks we recommend the new sets of input variables for NNs, which are shown below:

* $`Wb\overline{b}`$$`WH`$
NN: $`M_{b\overline{b}},P_T^{b1},P_T^{b2},P_T^{bb},\widehat{s},H_T^{jets},`$
$`cos(b1,b2)|_{lab},cos(b1,b1b2)|_{b1b2}`$
* $`WZ`$$`WH`$
NN: $`M_{b\overline{b}},P_T^{b1},P_T^{b2},H_T^{jets},`$
$`cos(b1,b2)|_{lab},Q\times cos(z,b1)|_{lab},cos(W,e)|_W`$
* $`t\overline{t}`$$`WH`$
NN: $`M_{b\overline{b}},M_{Wb},\widehat{s},M_{Wjetsb},H_T^{jets},`$
$`Q\times cos(\psi _{axis},e)|_{top},cos(b1,b1b2)|_{b1b2}`$
* $`tbj,tb`$$`WH`$
NN: $`M_{b\overline{b}},M_{Wb},P_T^{b2},\widehat{s},`$
$`M_{Wjetsb},P_T^{top},H_T^{jets},cos(z,e)|_{lab},`$
$`Q\times cos(z,b1)|_{top},cos(e,j)|_{top}`$

where there are three types of variables:

* “Singular” variables (denominator of Feynman diagrams):
$`M_{12}`$ is the invariant mass of two particles and/or jets (1 and 2) and corresponds to s-channel singularities;
$`P_T^f`$ (the transverse momenta of f);
$`M_{Wjetsb}`$ is the invariant mass of the $`W`$ and all jets except the $`b`$-jet for which the $`M_t=(p_W+p_b)`$ is closest to the top quark mass;
* “Angular” variables (numerator of Feynman diagrams, spin effects): $`cos(b1,b1b2)|_{b1b2}`$ means the cosine of the angle between highest $`P_T`$ b-quark and vector sum of the two highest $`P_T`$ b-quarks in the rest frame of these two b-quarks. Scalar (Higgs) and vector (gluon, Z-boson) particle decays lead to significantly different distributions on this variable, this is also very much different for the case when b-quarks come from the decay of top and anti-top quarks;
$`cos(b1,b2)|_{lab}`$ characterizes how much two b-quarks are collinear;
$`cos(z,b1)|_{lab}`$ and $`cos(W,e)|_W`$ reflect the difference in t-channel Z-boson and s-channel Higgs-boson production topologies where lab means the laboratory rest frame and $`z`$ means the z-axis;
$`cos(\psi _{axis},e)|_{top}`$  and $`cos(e,j)|_{top}`$  are the top quark spin correlation variables used in the analysis of the top quark pair and single production, the lepton charge $`Q`$ is added here to take uniformly into account the electron and the positron contributions from the $`W`$-boson decays.
* “Threshold” variables. As explained above the $`\widehat{s}`$ and $`H_T^{jets}`$ variables are used in our analysis.

As one can see from the Fig.2 using the new NN variables allows to improve the NN efficiency by a factor of 1.5-2 depending on the background process. It will lead to corresponding improvement in prospects to find a light Higgs at the Tevatron. However, one needs to take into account the ZH production channel as well as a number of detector efficiencies in order to predict a realistic discovery limit.

## Acknowledgements

The work was partly supported by the INTAS 00-0679, CERN-INTAS 99-377, YSF 02/239 and Universities of Russia UR.02.03.002 grants. E. Boos thanks the Humboldt Foundation for the Bessel Research Award.


text:
# Observations of Fourteen Pulsar Glitches

## 1 Introduction

Two kinds of irregularities are observed in the rotation rates of pulsars, timing noise and glitches. Timing noise is a continuous wandering of the rotation rate, while glitches are characterized by a sudden increase in the rate, often followed by a period of relaxation. They are often revealed by the sudden onset of continually decreasing arrival time residuals. Glitches were first observed in the Crab and Vela pulsars \[Boynton et al. 1969, Radhakrishnan & Manchester 1969\] and it was soon realized that they can be important diagnostic tools for studying neutron star interiors \[Ruderman 1969, Baym et al. 1969\]. It is widely believed that these events are caused by sudden and irregular transfer of angular momentum from a faster rotating interior superfluid to the solid crust of the neutron star. The result is a sudden fractional increase in the rotational frequency $`\nu _0`$ of the pulsar with a magnitude in the range $`10^{10}<\mathrm{\Delta }\nu _0/\nu _0<5\times 10^6`$. A characteristic feature of many glitches is a relaxation after the frequency jump, which may occur over a period of days, months or years. However, as we present in this paper, small glitches seem to show little significant relaxation after the glitch.

Although glitches are rather rare phenomena, an increasing number of these events have been reported recently (McKenna & Lyne 1990; Shemar & Lyne 1996; Wang et al. 2000; this paper), permitting more comprehensive statistical studies (Alpar & Baykal 1994; McKenna & Lyne 1990; Lyne, Shemar & Smith 2000). In this paper, we extend the analysis of the Jodrell Bank database by 2.5 year beyond that of Shemar & Lyne (1996) and seek to lower the detection threshold significantly.

## 2 Observations and analysis

Observations were carried out at Jodrell Bank, mostly using the 76-m Lovell radio telescope, but also occasionally using the 30-m Mark II telescope. Measurements were made at intervals of between one and three months, while some, more interesting, pulsars were observed more often. The total list of pulsars that have been observed regularly in this programme at Jodrell Bank is presented in Shemar & Lyne (1996, their Table 1). Those authors analysed the data up to about MJD 49500. This work represents a more detailed study of the same data and also extends the analysis on their list of pulsars to about MJD 50500. The B1950 and J2000 names of these pulsars, their periods and characteristic ages, and the dates spanned by the observations are listed in Table 1. The main improvement in the analysis is a more careful calibration of the systematic effects which arose from the use of a variety of filterbanks and dedispersion procedures over the typically 16-year span of the observations. Greater care was also exercised in removing data which might have been affected by impulsive radio-frequency interference.

Both telescopes were equipped with dual-channel cryogenic receivers at observing frequencies centered close to 408, 610 or 1400 MHz. Each channel was sensitive to one hand of circular polarization. The data were dedispersed using filterbanks and folded synchronously with the nominal topocentric rotation period of the pulsar for sub-integration periods of between one and three minutes. An observation consisted typically of six such integrations which were stored on disk for subsequent processing.

Total intensity profiles were obtained by adding the six sub-integrations. These were then cross-correlated with a standard template to give pulse topocentric times of arrival which were then corrected to the Solar system barycentre using the JPL ephemeris DE200 \[Standish 1982\]. Assessment of arrival time residuals, which are the differences between actual pulse arrival times and times calculated from a simple rotational model, provides information about the behaviour of the pulsar rotation. The fitting procedure used a simple slow-down model involving rotational frequency and its first derivative. The analysis for each pulsar involves such a fit to a period of time which is devoid of any glitch activity. The timing residuals for the whole data set are then inspected visually for the presence of glitches. Pulsar positions assumed in this analysis are given in Table 2. Several of the positions were obtained from the same timing data described here, using data located well away from glitches (Martin 2001).

The parameters of each glitch were obtained from comparison of parameters before and after the glitch. Epochs of glitches were determined by requiring a continuity of phase across the glitch. Pre-glitch parameters were obtained by fitting a simple pre-glitch slow-down model of the form $`\nu (t)=\nu _0+\dot{\nu }_0t`$ to the data. The observed post-glitch frequency residuals are described as a function of the time $`t`$ elapsed since the epoch of the glitch, relative to the pre-glitch ephemeris:

$$\mathrm{\Delta }\nu (t)=\mathrm{\Delta }\nu _p+\mathrm{\Delta }\dot{\nu }_pt+\ddot{\nu }_pt^2/2+\mathrm{\Delta }\nu _1e^{t/\tau _1},$$
(1)
where $`\mathrm{\Delta }\nu _p=\nu _p\nu _0`$ and $`\mathrm{\Delta }\dot{\nu }_p=\dot{\nu }_p\dot{\nu }_0`$ are differences between post-glitch and pre-glitch parameters. The last term in equation (1) describes an exponentially decaying transient component of post-glitch behaviour, while the penultimate term represents the large, approximately constant value of second derivative often seen following large glitches after any short-term transient has decayed. Note that $`\mathrm{\Delta }\nu _p`$ and $`\mathrm{\Delta }\dot{\nu }_p`$ usually differ from the instantaneous changes in $`\mathrm{\Delta }\nu _o`$ and its derivative because of the decaying components. Thus,

$$\mathrm{\Delta }\nu _p=\mathrm{\Delta }\nu _o\mathrm{\Sigma }\mathrm{\Delta }\nu _1\mathrm{and}\mathrm{\Delta }\dot{\nu }_p=\mathrm{\Delta }\dot{\nu }_o+\mathrm{\Sigma }\mathrm{\Delta }\nu _1/\tau _1$$
(2)
A more detailed description of the observation system and analysis can be found in the paper by Shemar & Lyne (1996).

## 3 Results

The parameters of 14 new glitches found in 9 pulsars are given in Table 3, which shows the epoch of the glitch as a Modified Julian Date (MJD), pre-glitch frequency $`\nu _0`$ and its first derivative $`\dot{\nu }_0`$ at that epoch, glitch fractional parameters, and the post-glitch frequency $`\nu _p`$ and its first derivative $`\dot{\nu }_p`$. The glitch epoch is estimated in the following manner. First, two solutions for the pulse phase across the glitch were obtained by extrapolating the pulse ephemeris before and after the glitch respectively. Then, the epoch was estimated from these phases by requiring that the the pulse phase be continuous across the glitch. The errors quoted in brackets are twice the standard deviations obtained from the formal fits. However, this procedure could not be used for the glitch in PSR B1930+22 due to lack of sufficient number of measurements near the glitch as well as the large magnitude of the glitch. The relaxation parameters, $`\dot{\nu }_p`$, $`\mathrm{\Delta }\nu _1`$ and $`\tau _1`$ are usually insignificant in these glitches and are mentioned in the text where appropriate.

The frequency residuals for previously unpublished glitches are presented in the lower panels of Figs. 1, 3, 5$``$9 and Fig. 11, and in the top panel of Fig. 10. They were obtained by performing local fits over about 50 days to the arrival time data, and presented relative to a simple slow-down model. Since most of the glitches presented in this paper are small, we also usually show their timing residuals in the upper panels of Figs. 1, 3, 5$``$9 and Fig. 11 for clarity of presentation. The number of available timing measurements between two glitches was small in the case of PSR B1737$``$30 making it difficult to obtain the frequency residuals in the manner described above except in case of two glitches. The timing residuals in all these pulsars show the familiar negative change in the gradient after the glitch, corresponding to a spin-up.

Below we describe detailed results of the search for glitches in the improved Jodrell Bank pulsar timing data base. Following Shemar & Lyne (1996), we give for each pulsar the conventional B-name, J-name \[Taylor, Manchester & Lyne 1993\] and the characteristic age $`\tau =\nu /2\dot{\nu }`$.

### 3.1 PSR B0154+61 (J0157+6212, $`\tau =200`$ kyr)

Although data on PSR B0154+61 have been collected at Jodrell Bank for more than 10 years since MJD 46866, this pulsar suffered its first observed glitch near MJD 48504. The glitch was quite small with the size of the frequency jump equal to $`2.46\times 10^9`$. The frequency and timing residuals are shown in lower and upper panels of Fig. 1, respectively. This glitch was not reported by Shemar & Lyne (1996), as it was below their completeness level of $`5\times 10^9`$ in $`\mathrm{\Delta }\nu _0/\nu _0`$.

### 3.2 PSR B1737$``$30 (J1740$``$3015, $`\tau =20.63`$ kyr)

This pulsar exhibits frequent glitches and nine glitches have been reported in the past (McKenna & Lyne 1990; Shemar and Lyne 1996). Our analysis included observations carried out up to May 1999 for this pulsar (around 51300 MJD). We present four more glitches detected in these data. The timing residuals for the two smaller glitches are shown in Figure 2. The other two glitches were large with a $`\mathrm{\Delta }\nu _0/\nu _0`$ exceeding $`1.0\times 10^7`$. The timing and frequency residuals for these glitches are shown in Figs. 3a and 3b and the glitch parameters for all the four glitches are presented in Table 3.

The cumulative fractional change in the rotation rate, $`\mathrm{\Delta }\nu _0/\nu _0`$, for the pulsar is shown in Fig. 4. The dashed - dot line indicates the average rate of fractional spin-up due to glitches. The mean spin-up rate due to glitches, $`\dot{\nu }_{glitch}`$ (See Lyne et al. 2000), for this pulsar is about 1.4 percent of its spin-down rate. Thus, a fixed fraction 0.014 of the pulsar’s slowdown is reversed by glitch activity and this is consistent with statistical estimate in other pulsars (Lyne et al. 2000). Shemar and Lyne (1996) noted that there are two types of glitches and this is evident from this figure. The larger glitches occur typically 800 days apart whereas the typical separation for all the glitches is of the order of 300 days.

### 3.3 PSR B1758$``$03 (J1801$``$0357, $`\tau =4,400`$ kyr)

The rotational frequency and timing residuals for this pulsar are shown in Fig. 5. It suffered a glitch after 3 years of regular monitoring at Jodrell Bank and was not reported by Shemar & Lyne (1996), being below their sensitivity threshold. The size of the frequency jump is rather small, with the fractional increase $`\mathrm{\Delta }\nu _0/\nu _03\times 10^9`$. The only older pulsar which has been observed to glitch is PSR B1859+07 ($`\tau =4.5`$ Myr). It is difficult to determine the exact date of the glitch because of the large gap of 140 days between observations. An estimate of the glitch epoch by assuming the continuity of pulse phase across the glitch indicates that the glitch occurred sometime near the beginning of 1990 (around MJD 48016). It was probably followed by a small relaxation, visible in the lower panel of Fig. 5.

### 3.4  PSR B1758$``$23 (J1801$``$2306, $`\tau =60`$ kyr)

The frequency and timing residuals for this pulsar are shown in Fig. 6. A total of 6 glitches have been detected in PSR B1758$``$23 since MJD 46697. The first three glitches were reported by Kaspi et al. (1993), the next one by Shemar and Lyne (1996) and we show two more recent glitches here. Four of these glitches were also reported by Wang et al. (2000). The first occurred around MJD 50055 and the second about a year later. Both newly-reported glitches are of moderate size, having magnitudes of $`\mathrm{\Delta }\nu _0/\nu _0=22.6\times 10^9`$ and $`80.6\times 10^9`$, respectively. These values are consistent with those reported by Wang et al. (2000). The event reported in Shemar & Lyne (1996) is of similar size to the events presented here, while the glitches recorded by Kaspi et al. (1993) are one order of magnitude larger, but still substantially less than those seen in Vela and other youthful pulsars. Thus all glitches in this pulsar are rather small. The post-glitch data analysis does not show any significant relaxation in this pulsar.

### 3.5 PSR B1800$``$21 (J1803$``$2137, $`\tau =16`$ kyr)

Shemar & Lyne (1996) reported a giant glitch with the third largest magnitude of all at the end of 1990, while we have found a small glitch with $`\mathrm{\Delta }\nu _0/\nu _0=5.3\times 10^9`$ that occurred more than five years later around MJD 50269. The relatively dense data coverage around the event permitted the determination of the time of the event to within about one day. The timing and frequency residuals are presented in Fig. 7 (relative to data from about 1000 days between glitches). These plots show that the glitch was accompanied by a significant increase in the rate of slow-down.

### 3.6 PSR B1907$``$03 (J1910$``$0309, $`\tau =3,600`$ kyr)

This is one of the oldest pulsars known to glitch and it suffered two rather small glitches. One of them occurred around MJD 48241 and the other about 1000 days later. In Fig. 8a and 8b, one can see the timing and frequency residuals for this pulsar corresponding to the two glitches. The fractional frequency increases are $`0.6\times 10^9`$ and $`1.84\times 10^9`$, respectively. The first glitch is the smallest glitch known. Again, timing noise is probably the reason for irregular behaviour both before and after glitches. Both glitches reported here were below the threshold of the search of Shemar & Lyne (1996).

### 3.7 PSR B1917+00 (J1919+0021, $`\tau =2,600`$ kyr)

Recently, around MJD 50174, the pulsar PSR B1917+00 suffered a small glitch, with fractional frequency increase $`\mathrm{\Delta }\nu _0/\nu _0=1.29\times 10^9`$. There appears to be a small short-term relaxation after the glitch.

### 3.8 PSR B1930+22 (J1932+2220, $`\tau =40`$ kyr)

A giant glitch occurred in this pulsar between MJD 50244 and 50284. The fractional frequency increase, $`\mathrm{\Delta }\nu _0/\nu _0=4457\times 10^9`$, observed in this glitch is the second largest ever reported, marginally greater than the previous largest glitch in PSR B0355+54 \[Lyne 1987, Shemar & Lyne 1996\] and somewhat smaller than the recently reported glitch in PSR J1614$``$5047 (Wang et al. 2000). The frequency residuals and the frequency derivative for this glitch are shown in Figure 10. Removing the mean value of frequency from both pre-glitch and post-glitch data and expanding the frequency scale by a factor of 100 (Fig. 10b) reveals an almost linear relaxation of the rotation rate following a short term quasi-exponential relaxation. The exponential recovery can also be seen in the frequency derivative.

### 3.9 PSR B2255+58 (J2257+5909, $`\tau =1,000`$ kyr)

PSR B2255+58 exhibited a small glitch at MJD 49463. Timing and frequency residuals corresponding to this glitch are presented in Fig. 11. This glitch is one of the smallest known with $`\mathrm{\Delta }\nu _0/\nu _0=0.92\times 10^9`$.

## 4 Discussion

Shemar & Lyne (1996) found 25 glitches in 10 pulsars after analysing about 2500 years of pulsar rotation in Jodrell Bank pulsar timing data base. Continuing this work we have found another 14 glitches in 9 pulsars in a further 1000 years of pulsar rotation in the improved and larger data base covering over about 7 years up to April 1997 (that is around 50500 MJD). We found 6 pulsars glitching for the first time, which increases the number of known glitching pulsars to 31. We also report new glitches in PSRs B1737$``$30, B1758$``$23 and B1800$``$21, which were already known as glitching pulsars. We report the smallest glitch ever observed, in PSR B1907$``$03 with $`\mathrm{\Delta }\nu _0/\nu _0=0.6\times 10^9`$, and the second largest glitch, in PSR B1930+22 with $`\mathrm{\Delta }\nu _0/\nu _0=4457\times 10^9`$. The fractional frequency increase in most of the glitches described in this paper is of the order of $`10^9`$. Thus, we detected 4 glitches which were not reported by Shemar & Lyne (1996) as they were below their threshold of $`5\times 10^9`$ in $`\mathrm{\Delta }\nu _0/\nu _0`$. We also found 10 more glitches in the epoch interval not covered by them. For those glitches newly found in the Shemar & Lyne interval, the average, the maximum and the minimum values of $`\mathrm{\Delta }\nu _0/\nu _0`$ were 1.9, 2.9 and 0.6, while for those detected since then, the values were 646.2, 4457 and 0.92, respectively (all in units of $`10^9`$). This can be compared with corresponding values of the Shemar & Lyne (1996) search: 1072, 4368 and 1.2, respectively. They detected several giant glitches while our survey resulted mostly in the detection of rather small glitches.

The detection threshold of the present search for glitches varies significantly from one pulsar to another, depending upon the accuracy of the times-of-arrival, the density of the observations and the presence of timing noise intrinsic to the pulsar. We have conducted a number of simulations on a few dozen representative pulsars which had not glitched in our data set. The effects of a number of glitches of different magnitude in turn were introduced into the arrival time data of a pulsar. These data were then inspected using the glitch detection procedure described in section 2. As a result of these tests, we estimate that at least 90% of glitches with magnitude $`2\times 10^9`$ have been detected. This is a factor of about 2.5 smaller than the threshold of Shemar & Lyne (1996). The detection of several glitches below the threshold of those authors indicates that the frequency of occurrence of small glitches does not decrease for smaller sizes of glitch (e.g. Lyne et al. 2000).

We confirm the observation of Shemar & Lyne (1996) and Lyne et al. (2000) that the dominant effect of glitches, particularly the smaller ones, is a sudden increase in rotational frequency with very little or no recovery. The age range of the glitching pulsars is broad, from 16 thousand to 4.4 million years, including PSR B1758$``$03, which is now the second oldest glitching pulsar.

Acknowledgements This work is supported in part by the KBN Grant 2 P03D 015 12 of the Polish State Committee for Scientific Research.


text:
# Non glassy ground-state in a long-range antiferromagnetic frustrated model in the hypercubic cell

## I Introduction

The success in understanding the collective behavior of conventional magnets lead in the past decade to an increased interest in systems that exhibit a novel type of ordering as a consequence of “frustration”. In particular, geometrically frustrated systems (having as a basic building block triangles of antiferromagnetic bonds) have been the object of extensive theoretical studies (see  and references therein). It was found that frustration gives rise to spectacular and often unexpected behaviors at low temperatures due to the high degeneracy of the ground state. Recently, more and more experimental realizations of such geometrically frustrated systems have been achieved. 2D and 3D materials have been assembled from triangles and tetrahedra such that neighbouring geometrical units share either a common edge or a common corner . Such new materials have stirred a new debate, whether frustration alone is strong enough to destroy magnetic long range order. Spin glass behavior is usually associated to the presence of both disorder and frustration. It has been argued that frustration alone can lead to large ground-states degeneracies, accompanied by extensive ground-state entropies, but cannot produce a sufficiently “rough” free energy landscape necessary for the development of a thermodynamical ”glassy” state.

In contrast with the previously mentioned systems, non-geometrical frustration may arise as an effect of antiferromagnetic long-ranged interactions (e.g., dipolar, Coulomb, etc.). It is well known that the competition between short-range ferromagnetic exchange interactions and long-range antiferromagnetic ones in hypercubic Bravais lattices give rise to low temperature long range order associated with lamellar structures. Although those systems presents a very rich behaviour concerning both its equilibrium- and non-equilibrium properties, in all these cases the ground state presents a finite degeneracy. Even in the case of a 3D system with Coulomb interactions, where some evidence of glassy behaviour has been obtained, it is of purely dynamical origin, the ground state still having a finite degeneracy.

In this work we analize the zero and low temperature equilibrium properties of a long range antiferromagnetic spin model defined on a $`D`$-dimensional hypercubic cell. The associated Hamiltonian is related to a recently proposed complexity measure of Boolean functions in the context of neural networks learning processes and it is described in section II. Although the original motivation for studying the statistical mechanics of this model is related to the the field of Boolean functions complexity and neural networks , we found that it presents some very peculiar properties that provide new insights about the relationship between frustration and order-disorder low temperature properties. Spin models defined on hypercubic cells have been repeatedly used in the past to study both equilibrium and dynamical properties of spin glasses-. In that context, the main interest in hypercubic cell models is to analyze how their properties change as the dimension (and therefore the connectivity of the lattice) increases, thus resembling fully connected models like the Sherrington-Kirkpatrick one, whose general properties are rather well understood. One of the main questions today concerning spin glasses is whether the physical scenario that emerges in long-range fully connected models holds in systems with short range interactions or not. One step further towards the answer of that question is to consider systems with long range interactions that decay with the distance between spins. In this work we analize the effects of frustrating long range antiferromagnetic interactions in an hypercuibc cell without disorder. The interactions in our model decay with the Hamming distance between spins, a natural metric in hypercubic cells, instead of the Euclidean distance. Notice that the Hamming distance in the hypercubic cell is equivalent to the chemical distance (i.e., the number of steps required to go from one node to another along the shortest path). Hence, the present study may also provide some insigths about the possible effects of frustrating interactions we expect in models defined on lattices with complex network topologies, such as small world networks (see Ref. and references therein).

We show that, for some particular values of the dimension $`D`$ the frustrating long range interaction gives rise to a ground state that is neither disordered (i.e., non-glassy) nor completely ordered but, it is composed of two equal size anti-aligned antiferromagnetic domains. The ground state is thus partially ordered. Moreover, it has an infinite degeneracy in the thermodynamic limit (although this degeneracy grows logarithmically with the number of spins, and hence its entropy per spin goes to zero).

In section III we analyze the frustrated ground state of the system for the particular set of values of $`D`$ such that it exhibit global frustration. In section IV we analize the finite temperature behaviour of the model using Monte Carlo simulations. In section V we present Monte Carlo results for the ferromagnetic version of the model, which are of interest in the general context of long range interacting systems. Finally, we present some general conclusions in section VI.

## II The model

In Ref. a measure for the complexity of Boolean functions of $`D`$ Boolean variables was introduced in the context of feed-forward neural networks learning processes. The $`2^D`$ possible inputs of a Boolean function of this type are located at the corners of a $`D`$-dimensional hypercubic cell. If we associate a spin variable $`S_i=\pm 1`$ to each corner $`i`$ of the cell every possible Boolean function corresponds to a particular spin configuration of this system. The complexity measure is based on the concept that complex Boolean functions assign different outputs to similar inputs, where the similarity is measured by the Hamming distance between them. Thus, in the magnetic analogy, configurations where neighbour spins tend to be anti-aligned correspond to functions of increasing complexity. The complexity measure can be then mapped into minus the long range antiferromagnetic Hamiltonian

$$H=\underset{i,j}{}J(r_{ij})S_iS_j$$
(1)
where $`r_{ij}`$ is the Hamming distance between spins $`i`$ and $`j`$ in the hypercubic cell and

$$J(r)=\{\begin{array}{cc}\frac{1}{M(r)}\hfill & \text{if }rD/2\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(2)
$`M(r)`$ being the number of nearest neighbors at distance $`r`$ of a given spin and equals the binomial coefficient

$$M(r)=C(D,r)=\frac{D!}{(Dr)!r!}.$$
(3)
The distribution of Boolean functions with high complexity is then associated with the low energy spectrum of Hamiltonian (1) which can be characterized by the low temperature statistical mechanics behaviour of the associated spin model. The Hamiltonian (1) can be rewritten as a sum of terms, each of them taking into account the interactions between pair of spins at a different Hamming distance:

$$H=\underset{r=1}{\overset{r=D/2}{}}J(r)\underset{[i,j/Hamming(i,j)=r]}{}S_iS_j,$$
(4)
Notice that the number of interaction terms in Hamiltonian (4) depends on $`D`$; as $`D`$ increases a new term is included every time $`D`$ takes an even value. For instance, in the cases $`D=2`$ and $`3`$ only nearest neighbor spins interact; in the cases $`D=4`$ and $`5`$ an interaction between next nearest neighbor spins is added, etc.. Hence, the antiferromagnetic nature of the interactions introduce local frustration for $`D4`$, in the sense that not all the couplings inside groups of spins closely located can be simultaneously minimized. However, this does not necessarily generate global frustration, in the sense of a non-ordered ground state with high degeneracy, due to the particular form Eq.(2) of the interaction coupling $`J(r)`$. Global frustration appears only when the number of interaction terms included in the Hamiltonian (4) is even, that is, when the dimension $`D`$ is a multiple of four $`D=4l`$, $`l=1,2,3,\mathrm{}`$. In all the other cases the ground state is the antiferromagnetic one (i.e., the Nèel state with every pair of nearest neighbor spins anti-aligned), with degeneracy $`2`$ (corresponding to the inversion of all the spins). To see this let us calculate the energy of the Nèel state. Every term in Eq.(4) corresponding to a fixed Hamming distance $`r`$ equals in this case $`(1)^rNM(r)`$, $`N=2^D`$ being the total number of spins (notice that the interaction between every pair of spin in Hamiltonian (1) is counted twice for simplicity; this corresponds to a renormalization of the energy units by a factor two, and does not affect the results). Hence, the energy of the Nèel state is given by

$$E_N=N\underset{r=1}{\overset{r=D/2}{}}(1)^r$$
(5)
which equals zero when $`D/2`$ is even and $`E_N/N=1`$ when $`D/2`$ is odd. Moreover, in the last case is easy to see that a single spin flip increases the energy in a quantity $`\mathrm{\Delta }E=2`$. Therefore, the Nèel state is stable against small perturbations. Monte Carlo simulations presented in section IV provide further evidence that the ground state is the Néel one. Moreover, the finite temperature simulations show that there is a finite critical temperature below which the system presents long range antiferromagnetic order when $`D/2`$ is odd.

In the next section we analyze the frustrated ground state for the case $`D=4l`$, $`l=1,2,3,\mathrm{}`$.

## III The frustrated ground state

From exhaustive enumeration of the energy values for the case $`D=4`$, where the Hamiltonian contains two terms, it was found that the ground state of the system consist of two equal size antiferromagnetic domains with an energy value equal to $`0.5`$. Monte Carlo simulations in higher dimensions ($`D=4,8,12,16`$) have also shown that the energy value of the ground state of the system is equal to $`0.5`$ and has the same structure found for the case for $`D=4`$. (Long relaxation times and very low cooling were used in different simulations, all of them confirming the structure and energy value of the ground state). Energy calculations for states composed by $`4`$ and $`8`$ anti-aligned antiferromagnetic domains($`D=8,16`$) also yield values higher than $`0.5`$. Moreover, as we will demonstrate next, the energy per spin of the state composed by two antiferromagnetic domains equals $`0.5`$ for arbitrary values of $`D`$. Since the antiferromagnetic state has always zero energy, all these result suggests that the two antiferromagnetic domains structure has the minimal energy for arbitrary dimension $`D=4l`$, $`l=1,2,3,\mathrm{}`$. The total magnetization and staggered magnetization is zero for the ground state. In figure 1 a representation of the ground state for $`D=4`$ is depicted, where the state of the spins is indicated by filled or empty circles. The antiferromagnetic domains correspond to two hypercubes of dimension $`D1`$ each (see Fig. 1). Each spin in this state has $`D1`$ nearest neighbor spins anti-aligned, which belong to the same domain, and one nearest neighbor spin aligned, which belong to the other domain. Once we fixed two nearest neighbor spins to be aligned, the rest of the structure is completely determined. Since every spin has $`D`$ nearest neighbors, the degeneracy of the state is $`2D=2\mathrm{log}_2N`$, where the factor $`2`$ results from the inversion operation of all the spins.

We now demonstrate that the energy per spin $`E_d`$ of the state composed of two antiferromagnetic domains equals $`0.5`$ for any dimension $`D=4l`$, $`l=1,2,3,\mathrm{}`$. $`E_d`$ can be written as

$$E_d=\underset{n=1}{\overset{n=D/2}{}}E_d(n)$$
(6)
$`E_d(n)`$ being the interaction energy between an arbitrary chosen reference spin and all its neighbors located at a Hamming distance $`n`$.

We denote the two antiferromagnetic domains as regions A and B. Consider now an arbitrary spin $`S_1`$ that belongs to the domain A. Spin $`S_1`$ has $`D1`$ anti-aligned nearest neighbors (NN) spins in the domain A and only one first nearest neighbor spin $`S_2`$ in the domain B. Spin $`S_2`$ has the same orientation than the spin $`S_1`$. As the total number of NN spins in dimension $`D`$ is exactly $`D`$, the first term of Eq.(6) equals

$$E_d(1)=\frac{(D1)}{D}+\frac{1}{D}$$
(7)
which can be expressed as

$$E_d(1)=2\frac{C(D,0)}{C(D,1)}1$$
(8)
Let us now consider the next nearest neighbors (NNN) term $`n=2`$ in Eq.(6). Note that the NNN spins of $`S_1`$ in the domain B (anti-aligned to $`S_1`$) are the $`D1`$ spins NN of $`S_2`$ in B; the remaining spins NNN to $`S_1`$ in A have the same orientation of it. Since the total number of NNN spins is $`C(D,2)`$ we have that

$`E_d(2)`$ $`=`$ $`{\displaystyle \frac{(D1)}{C(D,2)}}+{\displaystyle \frac{C(D,2)(D1)}{C(D,2)}}`$ (9)
$`=`$ $`2{\displaystyle \frac{C(D,1)C(D,0)}{C(D,2)}}+1`$ (10)
Following the same analysis the $`3rd`$ term of Eq.(6) equals

$$E_d(3)=2\frac{C(D,2)C(D,1)+C(D,0)}{C(D,3)}1$$
(11)
and the $`nth`$ order term can be written as

$$E_d(n)=2\frac{_{k=0}^{n1}(1)^kC(D,k)}{C(D,n)}+(1)^n.$$
(12)
Equation (12) can be rewritten as

$$E_d(n)=(1)^{n+1}\left(2a(n)1\right)$$
(13)
where

$$a(n)(1)^{n+1}\frac{_{k=0}^{n1}(1)^kC(D,k)}{C(D,n)}.$$
(14)
It is s easy to verify that $`a(n)`$ satisfy the recurrence equation

$$a(n+1)=\frac{C(D,n)(1a(n))}{C(D,n+1)}=\frac{n+1}{Dn}\left(2a(n)1\right)$$
(15)
from which it can be proved by induction that $`a(n)=n/D`$. Replacing into Eqs.(6) and (13) we finally obtain

$`E_d`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{n=D/2}{}}}(1)^{n+1}\left({\displaystyle \frac{2n}{D}}1\right)`$ (16)
$`=`$ $`{\displaystyle \frac{2}{D}}{\displaystyle \underset{n=1}{\overset{n=D/2}{}}}(1)^nn={\displaystyle \frac{1}{2}}`$ (17)
where in the last step we have used that

$$\underset{n=1}{\overset{n=D/2}{}}(1)^nn=\frac{D}{4}$$
(18)
## IV Finite temperature equilibrium properties

We now present some numerical simulation results at finite temperature. Monte Carlo simulations were performed using heat bath dynamics. We first analyze the frustrated case $`D=4l`$, $`l=1,2,\mathrm{}`$. In Fig. 2 we show the results for the internal energy per spin for system sizes $`D=4,8,12`$. The different curves were obtained by starting at high temperature with random or ferromagnetic initial conditions (there were no difference between both cases in the final results) and slowly cooling, that is, the initial configuration at every subsequent temperature was taken as the last configuration of the previous temperature. At every temperature we left the system thermalize during a varying transient period, which was higher at low temperatures.

For low values of $`D`$ we see that the energy converges to the value $`0.5`$. Direct inspection of the final configurations obtained at $`T=0`$ shows always the conjectured state composed by two antiferromagnetic domains. Notice from Fig.2 that there seems to be a change of behaviour around $`T0.2`$. An analysis of the energy distribution at low temperatures showed that this effect is related to a gap between the ground state and the first excited states that decreases as $`D`$ increases and hence, it is a finite size effect. Both the magnetization and the staggered magnetization are always $`𝒪(N^{1/2})`$ at any temperature and also the specific heat does not show any type of anomaly. So the system does not (globally) order at any temperature.

We next performed Monte Carlo simulations for the antiferromagnetic system when $`D/2`$ is odd.

In figure 3 we show the absolute value of the staggered magnetization, the energy per spin and the staggered susceptibility vs temperature for $`D=6,10,14`$. These results not only confirms that the ground state is antiferromagnetic, but also show that the system undergoes a second order phase transition at a finite critical temperature $`T_c1`$. The comparison of the magnetization curves with the solution of the Curie-Weiss equation $`m=tanh(\beta m)`$ suggests that mean field behaviour is exact in this case at any temperature.

## V The ferromagnetic model

Finally we consider the ferromagnetic version of the model, that is, the model defined by the Hamiltonian (1) with $`J(r)=\frac{1}{N(r)}`$. This system is non-extensive, in the sense that the quantity

$$\varphi =\underset{ji}{}J(r_{ij})=D/2=\frac{\mathrm{log}_2(N)}{2}$$
(19)
diverges for $`N\mathrm{}`$, and therefore the thermodynamic limit is not defined. In this situation the magnetization per spin $`m(T)`$ and the internal energy per spin $`u(T)`$ are expected to scale as

$`m(T)`$ $``$ $`m^{}(T/\varphi )`$ (20)
$`u(T)`$ $``$ $`\varphi u^{}(T/\varphi )`$ (21)
when $`N1`$, $`m^{}(x)`$ and $`u^{}(x)`$ being scaling functions independent of $`N`$. In the case of $`d`$-dimensional hypercubic lattices with interactions that decay with the distance $`r`$ between spins as $`1/r^\alpha `$ the quantity $`\varphi `$ is proportional to

$$N^{}=\frac{1}{1\alpha /d}\left(N^{1\alpha /d}1\right)$$
(22)
This scaling factor generalizes the traditional scaling $`JJ/N`$ of the couplings in the Curie-Weiss model, i.e., a fully connected Ising spin system with distance-independent interactions, which corresponds to the $`\alpha =0`$ case. For $`\alpha >d`$ the scaling factor (22) becomes independent of $`N`$ in the thermodynamic limit and the system is extensive, while for $`\alpha d`$ (this corresponds, for instance, to the case of three dimensional systems wtih dipolar interactions that decay as $`1/r^3`$ for $`r1`$) we have that $`N^{}\mathrm{ln}N`$. Thus, the behavior of the present hypercell model is expected to be analog to that of the borderline case between short and long range interactions $`\alpha =d`$ in hypercubic Bravais lattices.

The data collapse of the rescaled quantities for different values of $`D`$ shown in Fig.(4) verify the scaling (21). Moreover, the comparison of the magnetization curves with the solution of the Curie-Weiss equation $`m=tanh(m/T^{})`$ (where $`T^{}T/\varphi `$) shows that mean field behavior is asymptotically exact for a system of this type when the thermodynamic variables are properly rescaled, as was previously conjectured for general systems with long-range ferromagnetic interactions. Also the specific heat calculation shown in Fig.(5) displays a convergence to a discontinuous curve at the rescaled critical temperature when $`N\mathrm{}`$, consistently with the mean field behaviour.

## VI Discussion

We analyzed the statistical mechanics of a long-range model defined on a $`D`$-dimensional hypercube, which is related to the complexity of Boolean functions. Besides its interest in Boolean functions complexity and neural networks fields, the model exhibits several features that contribute to the understanding of general long-range systems. Thus, the statistical mechanics of the model is of interest on its own.

In the antiferromagnetic (AF) version of the model, the particular form of the coupling terms in Eq.(2) actually define two different models, according to the value of $`D`$. When $`D/2`$ is an odd number the model presents a low temperature AF phase and undergoes a mean-field like second order phase transition at a finite critical temperature $`T_c=1`$. The most interesting situation appears when $`D/2`$ is and even number. In this case long range AF order is suppressed at any finite temperature, while the ground state is composed by two equal size anti-aligned AF domains; this state becomes infinitely degenerated in the thermodynamic limit, but its entropy per spin goes to zero. This result shows that the frustration produced by the interplay between long range interactions and lattice topology can produce large ground state degeneracies without the generation of a “rough” free energy landscape. This result is of importance in spin glasses studies, where models defined on an hypercube are used to analize the influence of large dimensionalities. We have shown that this type of lattice can introduce effects of large ground state degeneracies which are independent of the presence of disorder, at least in the case of long range interactions. Although the contribution of such effects to the entropy per spin goes to zero in the thermodynamic limit, its importance could not be negligible for finite systems, as in the case of numerical simulation studies. It is worth to stress that the antiferromagnetic model do not have a well defined thermodynamic limit for general dimension of the hypercell $`D`$. The system exhibit a well defined thermodynamic limit only when the dimension $`D`$ is restricted to a particular set of values, that is, to be a multiple of four or not. This peculiar behavior illustrates the type of subtleties that may arise when working with long range frustrating interactions in hypercubic cells.

Finally, in the ferromagnetic version of the model we verified a previous conjecture for non-extensive (ferromagnetic) systems, that is, all the thermodynamic behaviour of the model can be exactly described at any temperature by mean field theory if the thermodynamic variables are properly rescaled with the system size.

This work was partially supported by grants from Consejo Nacional de Investigaciones Científicas y Técnicas CONICET (Argentina), Agencia Córdoba Ciencia (Córdoba, Argentina) and Secretaría de Ciencia y Tecnología de la Universidad Nacional de Córdoba (Argentina). L.F. acknowledges support from MRC IRC (Oxford, UK). We acknowledge fruitful suggestions of one anonymous referee that helped to improve the quality of the manuscript.


text:
# On Friedmann integrals coincidence

##

Expansion of the two-component universe with arbitrary spatial curvature has been considered. It has been shown that the Friedmann integrals of the almost flat universe do not coincide.

As it is known, luminosity distance–redshift relationship for the type Ia supernovae reveals the presence of an acceleration of the cosmological expansion . This is possible in existence of a cosmological repulsion, which is described usually by $`\mathrm{\Lambda }`$-term in the Einstein equations or equivalently by the presence of a vacuum-like medium with a negative pressure. The expansion becomes accelerated when the doubled vacuum density exceeds the decreasing matter density. These two components are described by two evolution constants being named Friedmann integrals. In recent works  it is claimed that the known cosmic coincidences  have to be supplemented by a Friedmann integrals coincidence, and various reasons and consequences of this coincidence are investigated. In the present work it is shown that the Friedmann integrals of our almost flat Universe do not and cannot coincide.

The two-component universe dynamics is described by the Einstein field equation for the expansion factor $`a`$,

$$\dot{a}^2=\frac{A_m}{a}+\frac{a^2}{A_v^2}k,$$
(1)
where $`k=1,\mathrm{\hspace{0.17em}0},1`$ for positive, zero and negative spatial curvature in comoving with the cosmological expansion frame. In the case $`k=\pm 1`$ factor $`a`$ is the space curvature radius, and for the flat universe $`k=0`$ factor $`a`$ is defined up to arbitrary scale transformation. Here and below convention $`c=1`$ is used.

The Friedmann integrals $`A_m`$ and $`A_v`$

$$A_m=\frac{8\pi }{3}G\rho _ma^3,A_v^2=\frac{8\pi }{3}G\rho _v=\frac{\mathrm{\Lambda }}{3}$$
(2)
determine the universe evolution. Here $`\rho _m`$ is the dust-like matter energy density including a dark matter, $`\rho _v`$ is the vacuum energy density, $`\mathrm{\Lambda }`$ is the cosmological constant, $`G`$ is the gravitational constant. Define

$$U_m=\frac{A_m}{a},U_v=\frac{a^2}{A_v^2},U=U_m+U_v.$$

Graphs of the functions $`U_m(a)`$, $`U_v(a)`$, $`U(a)`$ are plotted in Fig. 1.

When $`k=1`$, the Friedmann integral $`A_m`$ is a maximum value of the curvature radius in the standard model $`\mathrm{\Lambda }=0`$, the integral $`A_v`$ equals to an initial value $`a`$ in the de Sitter universe $`\rho _m=0`$. At the point $`a=\stackrel{~}{a}`$,

$$\stackrel{~}{a}^3=\frac{1}{2}A_mA_v^2$$
(3)
the maximum $`\stackrel{~}{U}`$ of the potential barrier $`U(a)`$ is disposed,

$$\stackrel{~}{U}=U(\stackrel{~}{a})=\left(\frac{\alpha A_m}{A_v}\right)^{2/3},$$
(4)
where $`\alpha =3\sqrt{3}/2`$. At this point the balance of matter gravitation and vacuum antigravitation is achieved, $`\rho _m=2\rho _v`$, and the deceleration parameter $`q`$ changes the sign:

$$q=\frac{\ddot{a}a}{\dot{a}^2}=\frac{\stackrel{~}{a}^3a^3}{a^3aA_v^2+2\stackrel{~}{a}^3}.$$
(5)
Density parameters $`\mathrm{\Omega }_m=\rho _m/\rho _c`$, $`\mathrm{\Omega }_v=\rho _v/\rho _c`$, where $`\rho _c`$ is the critical density, $`\rho _c=3H^2/8\pi G`$, $`H=\dot{a}/a`$ is the Hubble constant. The Einstein field equation (1) connects the density parameters and the Friedmann integrals by the correlation

$$\left(\mathrm{\Omega }1\right)^3=\left(\frac{A_v}{A_m}\right)^2\mathrm{\Omega }_m^2\mathrm{\Omega }_vk^3,$$
(6)
where $`\mathrm{\Omega }`$ is the total density parameter, $`\mathrm{\Omega }=1`$ for the flat universe; in the model (1) $`\mathrm{\Omega }=\mathrm{\Omega }_m+\mathrm{\Omega }_v`$.

As it is clear from Fig. 1, the Friedmann integrals coincide if the curves $`U_m(a)`$ and $`U_v(a)`$ intersection point is on the level $`U=1`$. If the level line $`U=1`$ passes through the potential barrier summit then $`A_v=\alpha A_m`$. In the case $`|\stackrel{~}{U}|>1`$ i.e. $`A_v>\alpha A_m`$ the closed universe cannot overcome the potential barrier and do not reach the accelerated expansion state .

According to data adduced in , $`A_v`$ does not differ from $`A_m`$ more than one order of magnitude, $`A_v`$ being larger than $`A_m`$, and the expansion factor now $`a_0`$ coincides with $`A_v`$. For our almost flat Universe $`\mathrm{\Omega }_v0.7`$, $`\mathrm{\Omega }_m0.3`$. As it is clear from the correlation (6), irrespective of the curvature sign the Friedmann integrals ratio $`A_v/A_m`$ is connected with the degree of the universe spatial flatness: $`\mathrm{\Omega }1A_v/A_m0`$, and the barrier summit $`\stackrel{~}{U}`$ (4) goes down from the level $`U=1`$. Because of that for the almost flat universe $`a_0`$ and $`A_v`$ do not coincide either:

$$\left(\frac{a_0}{A_v}\right)^3=\frac{\mathrm{\Omega }_v}{\mathrm{\Omega }_m}\frac{A_m}{A_v}1,$$
(7)
and $`a/A_v\mathrm{}`$ when $`\mathrm{\Omega }1`$. For all this

$$\left(\frac{a_0}{A_m}\right)^3=\frac{\mathrm{\Omega }_v}{\mathrm{\Omega }_m}\frac{A_v}{A_m}1.$$
(8)
Last data on anisotropy of the cosmic microwave background restrict spatial curvature of the Universe  $`|\mathrm{\Omega }1|0.03`$. According to (6) it corresponds $`A_v/A_m210^2`$. Emphasize that the foregoing statements on the value of the Friedmann integrals ratio are correct for $`k=\pm 1`$. A flat space is the limit of the curved one when $`A_v/A_m0`$. The latter is also clear from the fact that $`A_v`$ does not depend on the curvature radius and $`A_ma^3`$ (2). For the formally flat universe, when $`k=0`$, the equation (6) turns into an identity and in equation (1) the parameter $`a`$ becomes an arbitrary scale factor. Redefining the latter one can make any value of $`A_m`$ and the ratio of the Friedmann integrals loses the sense.


text:
# Optimal Size of a Complex Network

## Abstract

We investigate the response behavior of an Ising system, driven by an oscillating field, on a small-world network, with particular attention to the effects of the system size. The responses of the magnetization to the driving field are probed by means of Monte Carlo dynamic simulations with the rewiring probability varied. It is found that at low and high temperatures the occupancy ratio, measuring how many spins follow the driving field, behaves monotonically with the system size. At intermediate temperatures, on the other hand, the occupancy ratio first grows and then reduces as the size is increased, displaying a resonance-like peak at a finite value of the system size. In all cases, further increase of the size eventually leads to saturation to finite values; the size at which saturation emerges is observed to depend on the temperature, similarly to the correlation length of the system.

It has been known that under appropriate circumstances, the presence of noise in a driven system can enhance rather than suppress the output of the system. Such attractive phenomena, called stochastic resonance (SR), have been widely investigated in various systems of many practical applications ref:SR:general ; ref:SR:practical . Recently, those SR phenomena have also been examined in the context of the system size resonance: Stochastic flips of the mean field are observed to depend on the system size, leading the linear response of the system to reach a maximum at a certain system size ref:SSR . Those works considered two systems of fully coupled noisy oscillators and one two-dimensional Ising system with nearest neighbor interactions. Namely, the underlying connection topology of dynamic variables was assumed to be regular, with either local or global connections. Meanwhile, recent studies of neuronal networks, computer networks, biochemical networks, and even social networks, have revealed that various real systems in nature possess quite complex structures, which can be described neither by regular networks nor by completely random networks ref:network . Accordingly, it is desirable to study effects of the system size on collective responses in the systems with the connection topology of complex networks, which can be made more realistic. In particular, the interplay between the system size and noise may be relevant in various biological systems such as neural networks and other cell networks, which consist of finite numbers of elements. For example, in the study of stochastic resonance in biological systems, optimal sizes of calcium ion channel clusters have been examined. Observed is that the clustering of the release channels in small clusters increases the sensitivity of the calcium response ref:Jung . This suggests a possible realization of the system, providing motivation for the investigation of the size resonance in the complex network structure.

In this paper we consider an Ising model on a complex network, specifically, on the Watts and Strogatz (WS) type small-world network ref:WS . It is well known that the WS network is characterized by a small characteristic path length $`\mathrm{}\mathrm{ln}N`$, where $`N`$ is the number of nodes constituting the network, and a large clustering coefficient. Both are commonly observed properties of real networks in nature. The WS network in this paper is constructed following Ref. ref:WS : A regular one-dimensional (1D) network of $`N`$ nodes is first constructed with local connections of range $`k`$, under periodic boundary conditions. At this stage each node on the network has $`2k`$ nearest neighbors. Next, each local link is visited once and, with the rewiring probability $`P`$, removed and reconnected to a randomly chosen node. After a whole sweep of the entire network, the average number of shortcuts in the network of size $`N`$ is given by $`NPk`$. Accordingly, the rewiring probability $`P`$ may be regarded as the fraction of the average number of shortcuts over the total number of connections $`Nk`$. In this paper, the local interaction range $`k`$ is set equal to two for convenience; longer ranges ($`k>2`$) are not expected to lead to any qualitative difference. After the WS network is built as above, an Ising spin is put on every node, and an edge (or a link) connecting two nodes is regarded as the coupling between the two spins at the two nodes. Finally, we apply an oscillatory field, driving the Ising spins, and the corresponding responses of the average spin, i.e., the magnetization are probed via Monte Carlo dynamic simulations, with attention to the effects of the system size.

The Hamiltonian for the field-driven Ising model on the WS network, which is constructed as described above, reads

$$H=\frac{1}{2}\underset{i,j}{}J_{ij}\sigma _i\sigma _jh(t)\underset{i}{}\sigma _i,$$
(1)
where the ferromagnetic spin-spin interaction strength $`J_{ij}`$ is given by

$$J_{ij}=J_{ji}\{\begin{array}{cc}J\hfill & \text{for }j\mathrm{\Lambda }_i\text{,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(2)
The neighborhood $`\mathrm{\Lambda }_i`$ of node $`i`$ stands for the set of nodes connected to $`i`$ (via either local edges or shortcuts), and $`\sigma _i(=\pm 1)`$ represents the Ising spin at node $`i`$. The sinusoidally oscillating magnetic field $`h(t)=h_0\mathrm{cos}\mathrm{\Omega }t`$ is applied with the driving amplitude $`h_0`$ and frequency $`\mathrm{\Omega }`$, while the system is assumed to be in contact with an isothermal heat bath at temperature $`T`$. We probe the dynamics of the system described by Eq. (1) by means of Monte Carlo (MC) dynamic simulations, employing the heat bath single spin-flip algorithm ref:Newmanbook  and measuring the time $`t`$ in units of the MC time step. For thermalization, we start from sufficiently high temperatures and lower the temperature $`T`$ slowly with the decrement $`\mathrm{\Delta }T=0.01`$ (in units of $`J/k_B`$ with the Boltzmann constant $`k_B`$). The driving amplitude and frequency are taken to be $`h_0=0.1`$ and $`\mathrm{\Omega }=0.001`$. We have also considered different frequencies, for example, $`\mathrm{\Omega }=0.01`$ and $`0.1`$, and found that the resonance-like peak indicating the system size resonance behavior tends to diminish at higher frequencies (see below). While simulations are performed at a given temperature, the data from the first $`10^5`$ MC steps are discarded, which turns out to be sufficient for stationarity, and measurements are made for next $`10^5`$ MC steps. Networks of various sizes are constructed as described above, and averages are taken over 100 different network realizations for each size.

To investigate the collective response of the system, we measure the occupancy ratio $`R`$ which is defined to be the average fraction of the spins in the direction of the external field ref:ORfirst ; ref:BJKim :

$$R\frac{\text{number of spins in the direction of}h(t)}{\text{total number of spins}},$$
(3)
where $`\mathrm{}`$ denotes the time average. In other words, $`R`$ measures how many spins follow the oscillating magnetic field. It is easy to understand that $`R`$ approaches the value 1/2 in both low- and high-temperature limits (see, e.g., Ref. ref:SR\_WS\_1 ; ref:BJKim ) and becomes increased near the stochastic resonance temperature, reflecting that more spins follow the external driving. Such SR phenomena have been observed in the system of given size, and it has been demonstrated that the matching condition of two time scales, the relaxation time of the system and the inverse frequency of the driving field, yields the optimal noise strength $`T_{SR}`$ at which the system displays maximum responses ref:SR\_WS\_1 ; ref:BJKim . Here we consider the system studied in Ref. ref:SR\_WS\_1  from a different point of view and examine the behavior of the occupancy ratio with the system size at various temperatures, probing the size resonance.

In the absence of long-range shortcuts ($`P=0`$), the network structure reduces to that of the 1D regular network with only local couplings of range $`k`$. Accordingly, when $`P=0`$, the field-driven Ising model described by Eq. (1) as well as the undriven model \[$`h(t)=0`$\] does not exhibit long-range order at finite temperatures, yielding $`T_c=0`$. Note, however, that even such a 1D system displays SR behavior at finite temperatures ref:Brey . In the presence of long-range shortcuts ($`P0`$), on the other hand, it has been found that both driven and undriven Ising model display ferromagnetic order at finite temperatures ref:SR\_WS\_1 ; ref:smallIsing ; furthermore, double SR peaks, which originate from matching of two time scales have been observed.

We first consider the case without long-range shortcuts ($`P=0`$), i.e., the purely 1D system, and show in Fig. 1 the behavior of the occupancy ratio $`R`$ versus the system size $`N`$ at various temperatures. It is observed that the occupancy ratio $`R`$ first increases and eventually saturates to a finite value as the system size $`N`$ is increased. The saturation size $`N_s`$, beyond which the occupancy ratio $`R`$ does not increase any more, reduces as the temperature $`T`$ is raised. Figure 2 displays such temperature dependence of the saturation size $`N_s`$, which has been taken as the size giving the occupancy ratio with the difference from the stationary value less than $`5\times 10^4`$. One can observe the exponential behavior: $`N_se^{c/T}`$ with $`c=5.6\pm 0.7`$, which is reminiscent of the behavior of the correlation length. In the 1D Ising model the correlation length $`\xi `$ diverges in the low-temperature limit as $`\xi e^{k(k+1)/T}`$ ref:smallIsing . Here the local interaction range $`k`$ has been chosen to be two $`(k=2)`$, leading to the behavior $`\xi e^{6/T}`$, essentially the same as that of the saturation size. It is thus concluded that the saturation behavior emerges when the system size reaches the correlation length of the system.

Meanwhile, in the presence of long-range shortcuts $`(P0)`$, substantially different behavior has been obtained for the occupancy ratio. In Fig. 3, the occupancy ratio $`R`$ of the system with the rewiring probability $`P=0.5`$ is displayed as the system size $`N`$ is varied. At low temperatures $`(T1.5)`$, $`R`$ first decreases monotonically with the size $`N`$ and then saturates to the value $`0.5`$; at high temperatures ($`T2.9`$), on the contrary, $`R`$ increases monotonically to the saturation value larger than $`0.5`$, depending on the temperature. In contrast to these monotonic behaviors, at intermediate temperatures $`(1.5T2.9)`$, the occupancy ratio $`R`$ behaves non-monotonically, exhibiting a maximum at a finite value of the system size $`N`$. The height of such a resonance-like peak tends to increase as the temperature is raised. We have also considered different values of the rewiring probability $`P`$ as well as of the driving frequency $`\mathrm{\Omega }`$. It is found that as $`P`$ is increased, the range of the temperature, in which the size-resonance behavior is displayed, becomes wider and that the saturation temperature beyond which $`R`$ shows saturation behavior increases. On the other hand, as the driving frequency is increased, the position of the resonance peak shifts toward smaller values of the system size, thus tending to yield monotonic decrease of $`R`$.

Note that this resonance behavior manifests two kinds of length scale in the system with long range shortcuts: the saturation size $`N_s`$ and the resonance size $`N_m`$ at which $`R`$ reaches the maximum. To understand the possible relation with the correlation length even in the presence of long range shortcuts, we examine the behaviors of $`N_s`$ and $`N_m`$, which are displayed in Figs. 46.

Figure 4 exhibits the saturation size $`N_s`$ versus the temperature $`T`$ in the system with $`P=0.5`$. For comparison with the correlation volume described by $`\xi _V|TT_c|^{\overline{\nu }}`$ ref:XY , the data points are plotted in the logarithmic scale, thus fitted to a linear relation between $`\mathrm{ln}N_s`$ and $`\mathrm{ln}|TT_c|`$ with the proportionality constant (slope) $`\overline{\nu }`$. From this fitting, we obtain $`T_c2.7`$, with which the slope is estimated to be $`\overline{\nu }=2.0\pm 0.1`$. Noting that in $`d`$ dimensions, the correlation volume relates with the correlation length $`\xi `$ via $`\xi _V\xi ^d`$ and the behavior $`\xi |TT_c|^\nu `$, we thus have $`\overline{\nu }=d\nu `$ in a $`d`$-dimensional system. Here it is known that the (effective) dimension $`d`$ of a mean-field system should be taken as the upper-critical dimension $`d_u`$ ref:Botet , leading to $`\overline{\nu }=d_u\nu `$. With $`d_u=4`$ and $`\nu =\nu _{MF}=1/2`$ for a mean-field system, we conclude that the value $`\overline{\nu }=2.0`$ indicates a transition of the mean-field nature ref:XY . We also examine the other length scale $`N_m`$ and show its temperature dependence in Fig. 5, where for convenience $`N_s`$, shown already in Fig. 4, is also plotted. It is shown that both the two length scales behave similarly with the temperature, yielding essentially the same value of the exponent $`\overline{\nu }=2.0\pm 0.1`$ (see the slopes of the two fitted lines). Accordingly, both the two length scales $`N_s`$ and $`N_m`$ apparently measure the correlation length of the system. On the other hand, the fitting parameter $`T_c`$ for $`N_m`$ turns out to be $`3.1`$, which is somewhat higher than the value $`2.7`$ obtained from fitting of $`N_s`$. At this stage it is difficult to discern unambiguously the difference; more extensive simulations and careful analysis are necessary for confirming the origin as well as the presence of this discrepancy.

Finally, we consider systems with different rewiring probabilities, and examine how the rewiring probability $`P`$ affects the resonance behavior. Figure 6 shows the behavior of the resonance size $`N_m`$ with the temperature at various values of $`P`$. It is observed that $`N_m`$ first increases slowly with the temperature $`T`$ then very fast as $`T`$ approaches $`T_c`$, which again reminds us of the behavior of the correlation length. Behavior of $`N_m`$ with the rewiring probability $`P`$ at given temperature $`T=2.2`$ is displayed in the inset of Fig. 6. Note the rather fast decrease of $`N_m`$ for small rewiring probabilities ($`P0.5`$) and the saturation behavior for large rewiring probabilities ($`P0.5`$). Such saturation behavior has also been reported in the synchronization of the system of coupled oscillators on a WS network ref:Hong . It is also noteworthy that $`N_m`$ decreases as shortcuts are introduced to the system; the shortcuts tend to decrease the optimal system size which corresponds to the maximum collective response at given temperature and rewiring probability.

In conclusion, we have investigated the effects of the system size on the collective response, measured by the occupancy ratio, in the oscillatory field-driven Ising model on WS networks. In the purely one-dimensional system without long-range interactions, which does not undergo a phase transition at any finite temperature, the occupancy ratio has been found to display monotonic behavior, not exhibiting a resonance peak. As long-range interactions come into the system, on the other hand, system-size resonance, characterized by nonmonotonic behavior, has been observed to emerge, thus suggesting a possible relation between the size resonance and a finite-temperature phase transition. The resonance size at which the occupancy ratio reaches the peak may be regarded as the optimal size of the network, in view of the maximum response. It is noteworthy that the optimal size as well as the saturation size displays temperature-dependent behavior, which is essentially the same as that of the correlation length of the system. This suggests the interesting possibility of estimating the correlation length from the stochastic resonance behavior at various sizes. Namely, the size resonance phenomena may be used as a tool to measure the correlation length. Note also that at given temperature both length scales, the optimal size and the saturation size, tend to decrease as the amount of long-range interactions is increased.

As a possible application of the system size resonance, we suggest biological systems such as the assembly of beta cells which reside in the islets of Langerhans in a pancreas ref:Lange . It is known that the beta cells form clusters, each with a finite number of cells, rather than gathering together as one unit. Thus speculated is the possibility that the function of the beta cells may be optimized via the mechanism of the system-size resonance. In addition, the system size resonance behavior may also be useful for understanding the formation of the public opinion ref:opinion .

###### Acknowledgements.

H.H. thanks J. Lee for providing the privilege of using the computing facility Gene. B.J.K. was supported in part by the Korea Science and Engineering Foundation through Grant No. R14-2002-062-01000-0. M.Y.C. thanks the Korea Institute for Advanced Study for hospitality during his visit, where this work was performed, and acknowledges the partial support from the Swiss-Korean Outstanding Research Efforts Award Program.