text:
# Ground state of the Bethe-lattice spin glass and running time of an exact optimization algorithm

## I Introduction

Recent years have seen an increasing interaction between the fields of combinatorial optimization and statistical physicshayes1997 ; nature1999 ; TCS2001 . On one hand, several problems in the statistical physics of disordered systems have been mapped onto combinatorial problems, for which fast combinatorial optimization algorithms are available duxbury ; opt-phys2001 . This has provided valuable insights into questions that are hard to investigate with traditional techniques, such as Monte Carlo simulations. On the other hand, concepts and methods from statistical physics are increasingly applied to combinatorial optimization TCS2001 .

Easy/hard thresholds analogous to phase transitions have been observed in random instances of optimization and decision problems, including satisfiability (SAT) MiSeLe ; MoZe , vertex-cover cover  (VC), number partitioning mertens , and others. There is currently much interest in understanding how phase transitions affect the performance of combinatorial algorithms, following the observationcheeseman  that the averageaveragecase  or typical (i.e. median) running time of some algorithms exhibits a sharp change in the vicinity of a phase transition. For example, in 3SAT, a special case of SAT, and in VC, the typical running time of exact backtracking algorithms changescocco2001 ; cover-time  from a polynomial dependence on the input size in the ”solvable” region, to exponential dependence in the ”unsolvable” region. This provides an insight into the performance of algorithms that goes beyond the worst-case running time traditionally considered in complexity theory. (Note, however, that from the behavior of individual algorithms, strictly speaking, one cannot draw conclusions about the “typical hardness” of a problem itself). Recently, statistical physics techniques have been fruitfully applied to study easy/hard transitions and algorithmic performance TCS2001 .

In this paper, we apply a branch-and-cut algorithm, a technique developed in combinatorial optimization, to find the ground state of the Ising spin glass on random graphs with fixed coordination number (also called Bethe lattices).

Our motivation is twofold. The first goal is algorithmic: we want to characterize the typical running time of our algorithm, notably its behavior across the zero-temperature spin-glass/ferromagnet phase transition that occurs when varying the mean of the random couplings. The interest of this stems from the importance of branch-and-cut as a general technique in combinatorial optimization, and from the fact that finding the ground state of a spin glass is a prominent example of a hard optimization problem arising from statistical physics (in general, it is NP-hardbarahona1982 , see Section III). The performance of branch-and-cut for this application has not been investigated in detail before (see, however, Refs. simone96, , excited2d,  and bulk\_perturb, ), and here we fill this gap. An interesting aspect is that, unlike in SAT, VC and other classical combinatorial problems, here averaging over random instances is physically motivated. To our knowledge, the only other study relating a “physical” phase transition to algorithmic performance is that of Middletonmiddleton2002 , which investigates the typical running time of the matching algorithm for the random-field Ising model, which however is polynomial everywhere in the parameter space.

We find that the median running time of our algorithm varies sharply near the spin-glass/ferromagnet transition, indicating a change from polynomial time deep in the ferromagnetic phase, to slower than polynomial in the spin glass phase. We also observed a similar behavior for spin glasses on regular lattices in two and three dimensions, but will not report it here.

The second motivation for the present work lies in the ground-state properties of the Bethe-lattice spin glass, which recently have attracted a renewed interest mezard2001 ; mezard2002 ; parisi-tria ; boettcher . Using branch-and-cut, we compute the ground state energy and magnetization, and locate the spin-glass/ferromagnet phase transition. This provides a useful test of recently developed analytical methods to treat diluted spin glass models mezard2001 ; mezard2002 ; franz2002 . We solve the model in the Bethe-Peierls (BP) approximation (equivalent to the replica symmetric approximation in the replica formalism) using a stochastic approach proposed in Refs.mezard2002, ; thouless2, . By comparing the branch-and-cut results with the BP results, we estimate the magnitude of the replica symmetry breaking corrections to the ground state energy and magnetization, finding that they are small.

The rest of the paper is organized as follows. In Section II, we introduce the Bethe-lattice spin glass model. In Section III, we describe the branch-and-cut algorithm used to calculate the exact ground states of the model. In Section IV, we describe the Bethe-Peierls approximation and the stochastic procedure used to solve it. In Section V, we present our branch-and-cut and BP results for the ground state energy and the zero temperature phase transition. In Section VI, we show that this transition coincides with a change of the typical running time. Finally, Section VII summarizes our results.

## II Model

The system considered consists of $`N`$ Ising spins $`S_i=\pm 1`$ sitting on the nodes of a graph $`G=(V,E)`$, where $`V=\{1,\mathrm{},N\}`$ is the set of nodes and $`E=\{(i,j)\}V\times V`$ is the set of edges of the graph. The energy of the system is given by

$$=\underset{(i,j)E}{}J_{ij}S_iS_j$$
(1)
where the couplings $`J_{ij}`$ are independent, identically distributed random variables drawn from a Gaussian distribution $`P()`$ with mean $`\mu `$ and unit variance,

$$P(J)=\frac{1}{\sqrt{2\pi }}\mathrm{exp}[(J\mu )^2/2]$$
(2)
We consider the case in which $`G`$ is a random graph with fixed connectivity $`z`$, or $`z`$-regular graph, where each spin interacts with exactly $`z`$ neighbors. This provides a convenient realization of a Bethe lattice, which avoids some complications associated to the usual construction of a Bethe lattice from a Cayley tree mezard2001 . Frustration is induced by large loops, the typical size of a loop being of order $`\mathrm{log}(N)`$. Small loops are rare, giving the graph a local tree-like structure, and therefore the mean field approximation is exact. A related model is the Viana-Bray model viana1985 , in which the connectivity is a Poisson variable with finite mean. Finite-connectivity or “diluted” models provide a better approximation to finite-dimensional spin glasses than the infinitely-connected Sherrington-Kirkpatrick model. Furthermore, they are directly related to optimization problems such as graph partitioning and coloring.

Although it has long been known that replica symmetry is broken in these two models viana1985 ; mottishaw ; thouless , until recently a replica symmetry broken solution could be found only in some limit cases. Mézard and Parisi recently introduced mezard2001 ; mezard2002  a “population dynamics” algorithm which allows a full numerical solution at the level of one step of replica symmetry breaking. Explicit results were derived for the Bethe-lattice spin glass with the symmetric $`\pm J`$ disorder distribution, but not for the Gaussian distribution considered here or for a non-zero mean. Previous numerical studies of this model can be found in Refs. lai, ; boettcher, ; banavar1987, . For a complete discussion of the Bethe-lattice spin glass, see Ref. mezard2001,  and references therein.

## III Branch-and-cut algorithm

The problem of finding a ground state of the Hamiltonian in Eq. (1) is in general computationally demanding. For a generic graph $`G`$, it is NP-hardbarahona1982 ; opt-phys2001 . For NP-hard problems, currently only algorithms are available, for which the running time increases faster than any polynomial in the system size, in the worst case (see Section VI for a brief description of complexity classes). In the special case of a planar system without magnetic field, e.g. a square lattice with periodic boundary conditions in at most one direction, efficient polynomial-time matching algorithmsbieche1980  exist. For the square lattice with periodic boundaries in both directions, polynomial algorithms exist for computing the complete partition function for the $`\pm J`$ distributionsaul  and for the case in which the coupling strenghts are bounded by a polynomial in the system sizeGLV2001 . In practice, both algorithms can only reach relatively small system sizes.

For the Bethe lattice considered here (and for regular lattices in dimension higher than two), no polynomial algorithm is known. Heuristic algorithms recently used include simulated annealing kirkpatrick1983 , “multicanonical” simulationSG-berg92 , genetic algorithmsSG-pal1996b ; py , extremal optimizationboettcher , a hierarchical renormalization-group based approachhoudayer1999 , and the cluster-exact approximation algorithmalex-stiff . Here, however, we are interested in investigating the running time of an exact, deterministic algorithm, since in this case the running time to find the exact ground state is a well defined quantity. We study the branch-and-cut methodsimone95 ; simone96  (see Ref. JN2001,  for a tutorial on optimization problems and techniques, including branch-and-bound and branch-and-cut), which is currently the fastest exact algorithm for computing spin glass ground statesbulk\_perturb , with the exception of the polynomial-time special cases mentioned above. As the branch-and-cut method is basically branch-and-bound with cutting planes, we also did some experiments with a pure branch-and-bound algorithmhartwig84 ; klotz , which however can only deal with much smaller system sizes, finding a qualitatively similar, but less pronounced, variation of the running time across the transition.

In the rest of this Section, we repeat a short description of the branch-and-cut method already given in Ref. bulk\_perturb, , to the benefit of the reader. It is convenient to map the problem of minimizing the Hamiltonian in Eq.(1) into a maximum cut problem. Consider a graph $`G=(V,E)`$, and let assign weights $`\{K_{ij}\}`$ to the edges. Given a partition of the node set $`V`$ into a subset $`WV`$ and its complement $`VW`$, the cut $`\delta (W)`$ associated to $`W`$ is the set of edges with one endpoint in $`W`$ and the other endpoint in $`VW`$, namely $`\delta (W)=\{(ij)EiW,jVW\}`$. The weight of $`\delta (W)`$ is defined as the sum of the weights of the cut edges, $`_{(ij)\delta (W)}K_{ij}`$. The maximum cut is a node partition with maximum weight among all partitions. It can be shown simone95  that minimizing the Hamiltonian in Eq.(1) is equivalent to finding a maximum cut of $`G`$ with the assignment $`K_{ij}=J_{ij}`$.

The branch-and-cut algorithm solves the maximum cut problem through simultaneous lower and upper bound computations. By definition, the weight of any cut gives a lower bound on the optimal cut value. Thus, we can start from any cut and iteratively improve the lower bound using deterministic heuristic rules (local search and other specialized heuristics, see Ref. juen-vlsi1988,  for details). How do we decide when a cut is optimal? This can be done by additionally maintaining upper bounds on the value of the maximum cut. Upon iteration of the algorithm, progressively tighter bounds are found, until optimality is reached.

Since the availability of upper bounds marks the difference between a heuristic and an exact solution, we now summarize how the upper bound is computed (for more details, see Ref. juen-vlsi1988, .) To each edge $`(ij)`$ we associate a real variable $`x_{ij}`$ and to each cut $`\delta (W)`$ an incidence vector $`\chi ^{\delta (W)}^E`$ with components $`\chi _{ij}^{\delta (W)}`$ associated to each edge $`(ij)`$, where $`\chi _{ij}^{\delta (W)}=1`$ if $`(ij)\delta (W)`$ and $`\chi _{ij}^{\delta (W)}=0`$ otherwise. Denoting by $`P_C(G)`$ the convex hull of the incidence vectors, it can be shown that a basic optimum solutionchvatal1983  of the linear program

$$\mathrm{max}\{\underset{(ij)E}{}J_{ij}x_{ij}xP_C(G)\}.$$
(3)
is a maximum cut. In order to solve (3) with linear programming techniques we would have to express $`P_C(G)`$ in the form

$$P_C(G)=\{x^EAxb,0x1\}$$
(4)
for some matrix $`A`$ and some vector $`b`$. Whereas the existence of $`A`$ and $`b`$ are theoretically guaranteed, even subsets of $`Axb`$ known in the literature contain a huge number of inequalities that render a direct solution of (3) impractical.

Instead, the branch-and-cut algorithm proceeds by optimizing over a superset $`P`$ containing $`P_C(G)`$, and by iteratively tightening $`P`$, generating in this way progressively better upper bounds. The supersets $`P`$ are generated by a cutting plane approach. Starting with some $`P`$, we solve the linear program $`\mathrm{max}\{_{(ij)E}J_{ij}x_{ij}xP\}`$ by Dantzig’s simplex algorithmchvatal1983 . Optimality is proven if either of two conditions is satisfied: (i) the optimal value equals the lower bound; (ii) the solution vector $`\overline{x}`$ is the incidence vector of a cut.

If neither is satisfied, we have to tighten $`P`$ by solving the separation problem. This consists in identifying inequalities that are valid for all points in $`P_C(G)`$, yet are violated by $`\overline{x}`$, or reporting that no such inequality exists. The inequalities found in this way are added to the linear programming formulation, obtaining a new tighter partial system $`P^{}P`$ which does not contain $`\overline{x}`$. The procedure is then repeated on $`P^{}`$ and so on.

At some point, it may happen that (i) and (ii) are not satisfied, yet the separation routines do not find any new cutting plane. In this case, we branch on some fractional edge variable $`x_{ij}`$ (i.e. a variable $`x_{ij}\{0,1\}`$), creating two subproblems in which $`x_{ij}`$ is set to 0 and 1, respectively. We then we apply the cutting plane algorithm recursively for both subproblems.

## IV Bethe-Peierls approximation

We recall here the zero-temperature formulation of the BP approximation, loosely following Ref. mezard2002, . We consider the Hamiltonian Eq.(1) on a random graph with fixed connectivity $`z=k+1`$, in which however some spins $`S_i`$ (cavity spins) have only $`k`$ neighbors. The random couplings are drawn from a distribution $`P(J)`$. The BP approximation consists in assuming that the ground state energy of this system is given by $`E=const._ih_iS_i`$, where the sum runs over all cavity spins. The cavity fields $`h_i`$, implicitly defined by this relation, are independent, identically distributed random variables when considered as a function of the random couplings. Their distribution $`P(h)`$ is the central object of interest, and satisfies a recursion relation derived as follows. Suppose we add a new spin $`S_0`$ to the system, which interacts with $`k`$ pre-existing cavity spins $`S_1,\mathrm{},S_k`$ through couplings $`J_1,\mathrm{},J_k`$, and we minimize the energy with respect to $`S_1,\mathrm{},S_k`$. Now $`S_0`$ is a cavity spin, and it is easily shown that its cavity field $`h_0`$ is given by

$$h_0=\underset{i=1}{\overset{k}{}}u(J_i,h_i),$$
(5)
where $`u(J_i,h_i)=\frac{1}{2}\left(|h_i+J_i||h_iJ_i|\right)`$. This provides a recursion relation for $`P(h)`$ as the $`J`$’s fluctuate according to $`P(J)`$.

Given a spin $`S_0`$ interacting with $`k+1`$ neighbors with couplings $`J_1,\mathrm{},J_{k+1}`$, the internal field $`H`$ acting on $`S_0`$ is

$$H=\underset{i=1}{\overset{k+1}{}}u(J_i,h_i).$$
(6)
Therefore, if we know $`P(h)`$ we can determine the probability distribution $`P(H)`$, and the magnetization

$$m_{BP}=\underset{ϵ0^+}{lim}𝑑HP(H)\text{sgn}(H),$$
(7)
where $`\text{sgn}(x)`$ is the sign function and $`ϵ`$ is a small field that breaks the symmetry of Eq.(5) with respect to changing the sign of all cavity fields.

The knowledge of $`P(h)`$ is also sufficient to determine the ground state energy of the system. As shown in Ref. mezard2002, , this can be expressed as

$$e_{BP}=[\mathrm{\Delta }E^{(1)}]\frac{k+1}{2}[\mathrm{\Delta }E^{(2)}],$$
(8)
where $`[\mathrm{}]`$ is the expectation value with respect to $`P(J)`$ and $`P(h)`$, and the quantities $`\mathrm{\Delta }E^{(1)},\mathrm{\Delta }E^{(2)}`$ are given by

$$\mathrm{\Delta }E^{(1)}=\underset{i=1}{\overset{k+1}{}}a(J_i,h_i)|\underset{i=1}{\overset{k+1}{}}u(J_i,h_i)|$$
(9)
and

$$\mathrm{\Delta }E^{(2)}=\underset{S_i,S_j}{\mathrm{max}}(h_iS_i+h_jS_j+J_{ij}S_iS_j),$$
(10)
where $`a(J_i,h_i)=\frac{1}{2}\left(|h_i+J_i|+|h_iJ_i|\right)`$ and $`S_i,S_j`$ in Eq.(10) are two randomly chosen cavity spins which we connect with the coupling $`J_{ij}`$.

The BP recursion, especially at finite temperature, has been studied extensively (see Ref. mezard2001,  and references therein). In particular, Mézard and Parisi mezard2002  have given an analytic expression of $`P(h)`$ for a binary $`P(J)`$. Klein et al.klein  solved the finite temperature BP recursion for Gaussian couplings with $`\mu =0`$ in the vicinity of the spin-glass/paramagnet transition. No analytical solution has been derived for Gaussian couplings at $`T=0`$, to our knowledge, although Klein et al. klein  derived an analytic solution near the spin-glass/ferromagnet transition $`\mu =\mu _c`$ within the mean random-field approximation.

Here, we employ the stochastic iterative procedure proposed by Mézard and Parisimezard2002  for the more general one-step replica symmetry broken case (see also Ref. thouless,  for a previous application of a similar method). We consider a population of $`𝒩`$ sites, to which we associate $`𝒩`$ cavity fields, which are initially assigned at random (with a small positive bias). We then select $`k`$ sites at random, extract $`k`$ couplings from $`P(J)`$, compute $`h_0`$ from Eq.(5), and assign $`h_0`$ as the new cavity field of a randomly chosen site. We iterate this procedure $``$ times per site. After a certain number of iterations, the distribution $`P(h)`$ will satisfy Eq.(5). At each iteration, by merging $`k+1`$ randomly chosen sites we compute the internal field $`H`$ with Eq.(6), and $`\mathrm{\Delta }E^{(1)}`$ with Eq.(9), and by merging two sites we compute $`\mathrm{\Delta }E^{(2)}`$ with Eq.(10). After discarding the first $`/4`$ iterations, by averaging $`\text{sgn}(H)`$, $`\mathrm{\Delta }E^{(1)}`$ and $`\mathrm{\Delta }E^{(2)}`$ over the remaining iterations we compute the estimates of $`m_{BP}`$ and $`e_{BP}`$ from Eqs.(7) and (8), and their statistical error from a binning procedure. We repeated the procedure for many values of $`\mu `$, choosing $`=10^4`$ and $`𝒩`$ between $`10^3`$ and $`10^5`$, the larger population being for $`\mu `$ near the transition point $`\mu _c`$. With $`𝒩=10^5`$, the iteration requires about one hour of computer time.

We note that the BP approximation is known to be wrong, being equivalent to the replica symmetric solution which is unstable. We have not attempted to use the generalization of the above procedure to one step of replica symmetry breakingmezard2001 , since for the Gaussian case considered it would require a significant computing time, and since the BP approximation gives sufficiently accurate results for our purposes.

## V Results

We have studied the Ising spin glass on random graphs with fixed connectivity $`z=4`$ and $`z=6`$. The instance generator first builds a random regular graph with the algorithm described in Ref. SteWor1999, . We then assign the couplings $`J_{ij}`$ according to the distribution $`P(J)`$ in Eq.(2).

Using the branch-and-cut approach we were able to study graph sizes up to $`N=400`$ for $`z=4`$ and $`\mu 0.9`$, and up to $`N=200`$ for $`z=6`$ and $`\mu 0.7`$. For larger values of $`\mu `$, we considered sizes up to $`N=1280`$. The ground states for the smallest systems can be obtained within a second, while the longest computations lasted at most one day on a typical workstation. Incidentally, for Ising spin glasses on a regular grid, specialized heuristics exist that exploit the grid structure, making it possible to consider larger system sizes than for the model reported here. More one timing issues, in relation to the spin-glass/ferromagnet phase transition, is presented in Section VI.

All the results were averaged over many samples (a sample, or instance, is a realization of the random graph with a realization of the couplings). The largest number of samples were considered in the vicinity of the phase transition, where the fluctuations of the magnetization are larger. Near the transition, for sizes $`N240`$ ($`z=4`$) and $`N160`$ ($`z=6`$) we computed around 5000 samples for each value of $`\mu `$; for $`N=400`$ ($`z=4`$) and $`N=200`$ ($`z=6`$), around 500 samples for each value of $`\mu `$. For sizes larger than these, we computed up to 280 samples for each $`\mu `$. In the following analysis of the ground state energy and magnetization, we consider only sizes up to $`N=400`$ $`(z=4)`$ and $`N=200`$ ($`z=6`$), since for larger sizes the statistical error is quite large. In the analysis of running times we will include all sizes.

### V.1 Ground state energy

We start by showing, in Fig. 1, the average ground state energy $`E(\mu ,N)`$, divided by $`zN`$, as a function of $`\mu `$ for $`z=4,6`$ and two different system sizes. For sufficiently large $`\mu `$, the system is completely magnetized, therefore the ground state energy depends linearly on $`\mu `$, $`E(\mu ,N)/Nz\mu `$, as visible in the figure. For small $`\mu `$ the system is frustrated, hence the energy saturates. Note that $`E(0,N)`$ scales as $`\sqrt{z}`$, not as $`z`$, therefore the two curves diverge at small $`\mu `$. The lines in Fig. 1 represent the numerical solution of the BP recursion obtained with a population size $`𝒩=10^3`$ (we verified that with $`𝒩=10^5`$ the results are unchanged) and $`=10^4`$ iterations of the stochastic algorithm. Clearly, the branch-and-cut results agree well with the BP approximation.

We extrapolate the branch-and-cut results to $`N=\mathrm{}`$ by fitting the data with the form $`E/N=e_{\mathrm{}}+bN^{2/3}`$. As shown in Fig. 2, the finite size corrections are well described by a $`N^{2/3}`$ dependence for small $`\mu `$, although an $`N^\omega `$ correction fits reasonably well the data for other values of $`\omega `$ between 0.6 and 1 as well. For large $`\mu `$, the finite size corrections are very small. A $`N^{2/3}`$ correction was also found to fit well the numerical data by Boettcher boettcher , who computed the average ground state energy of the $`\pm J`$ model for $`z`$ up to $`z=26`$ and $`N`$ up to $`N=2048`$ using a heuristic algorithm. In Ref. mezard2001, , the finite-size dependence of the energy at $`T=0.8`$ was studied, for the $`\pm J`$ distribution and $`z=6`$, finding a finite-size exponent $`\omega =0.767(8)`$, not far from 2/3. For the Viana-Bray model with fluctuating connectivity with mean $`z=6`$, the value $`\omega =0.62\pm 0.05`$, compatible with 2/3, was foundpalassini-phd , also using a heuristic algorithm.

Fig. 2 also shows that the extrapolated energy, $`e_{\mathrm{}}`$, is very close to the BP result, $`e_{BP}`$, in the whole range of $`\mu `$. Of course, the agreement is not surprising for large $`\mu `$, where replica symmetry holds. For smaller $`\mu `$, the observed agreement indicates that replica symmetry breaking corrections to the ground state energy are small (less than $`1\%`$). A similar conclusion was reached in Ref. mezard2002,  for the $`\pm J`$ distribution with zero mean.

In particular, for $`\mu =0`$ we obtain $`e_{\mathrm{}}=1.38\pm 0.04`$ ($`z=4`$) and $`e_{\mathrm{}}=1.72\pm 0.02`$ ($`z=6`$), where the errors take into account the uncertainty on the correction exponent $`\omega `$, to be compared with our BP result $`e_{BP}=1.351\pm 0.002`$ ($`z=4`$) and $`e_{BP}=1.737\pm 0.002`$ ($`z=6`$). It is also interesting to compare this with the ground state energy per spin found in twoalex-2d-unpublished  and three dimensionsSG-pal1996b  (which have coordination number $`z=4`$ and $`z=6`$, respectively) with Gaussian couplings and $`\mu =0`$, which is $`e_{\mathrm{}}=1.31453(3)`$ and $`e_{\mathrm{}}=1.7003(1)`$ respectively.

### V.2 Ground state magnetization

In Figs. 3 and 4 the symbols show, for $`z=4`$ and $`6`$ respectively, the average ground state magnetization $`m=[M]_J`$, where $`M=\frac{1}{N}_iS_i`$ and $`[\mathrm{}]_J`$ denotes the sample average, as a function of $`\mu `$ for different system sizes $`N`$. The lines show the $`N=\mathrm{}`$ result in the BP approximation. For small $`\mu `$, the magnetization vanishes as $`1/\sqrt{N}`$. For large $`\mu `$, the finite-$`N`$ data agree with the BP result within the error bars, with negligible finite-size corrections (again, we recall that the BP approximation is exact for sufficiently large $`\mu `$, hence the agreement is expected). From the point at which the BP magnetization vanishes, we estimate the critical coupling strength $`\mu _c=0.742\pm 0.005`$ ($`z=4`$) and $`\mu _c=0.546\pm 0.005`$ ($`z=6`$). Note that recursion relation Eq.(5) admits two symmetric solutions for $`\mu >\mu _c`$. Hence, in the stochastic procedure the magnetization will oscillate between positive and negative values, with an oscillation “time” (number of iterations) $`_0`$ that increases with the population size $`𝒩`$ and with $`\mu `$. Therefore, to compute the magnetization correctly, we need $`_0`$. To do this, we increased the size of the population progressively from $`𝒩=10^3`$ to $`𝒩=10^5`$ as $`\mu `$ approached $`\mu _c`$. (Residual oscillations very close to $`\mu _c`$ introduce a small systematic errorfootsmaller , which is reflected in the errors for $`\mu _c`$ quoted above.)

Another estimate of $`\mu _c`$ can be obtained from the Binder cumulantbinder

$$g(\mu )=\frac{1}{2}\left(3\frac{[M^4]}{[M^2]^2}\right),$$
(11)
where $`[\mathrm{}]`$ is now the “time” average. In the limit $`𝒩\mathrm{}`$, $`g(\mu )=0`$ for $`\mu <\mu _c`$ and $`g(\mu )=1`$ for $`\mu >\mu _c`$, hence $`g(\mu )`$ can be used to locate $`\mu _c`$. As shown in Fig. 5, the variation of the Binder cumulant with $`\mu `$ sharpens as $`𝒩`$ increases, an effect of the sign oscillations of the magnetization, which become less important as $`𝒩`$ increases. From $`𝒩=10^5`$ we estimate

$`\mu _c^{BP}`$ $`=`$ $`0.743\pm 0.005(z=4)`$
$`\mu _c^{BP}`$ $`=`$ $`0.547\pm 0.005(z=6)`$
which agrees with the above estimate from the average magnetization. We also verified that with these values of $`\mu _c`$, the magnetization obeys $`m_{BP}=a(\mu \mu _c)^\beta `$ for $`\mu \mu _c`$, with the mean-field exponent $`\beta =1/2`$ and $`a0.23`$.

Klein et al. klein  solved the BP recursion in the vicinity of $`\mu _c`$ using the mean random field approximation (MRF). Their results $`\mu _c^{MRF}=0.775`$ ($`z=4`$) and $`\mu _c^{MRF}=0.587`$ ($`z=6`$) (obtained after rescaling their value by an appropriate normalization factor $`\sqrt{z}`$) are slightly larger than our result $`\mu _c^{BP}`$.

In order to obtain an estimate of $`\mu _c`$ from the finite-$`N`$ branch-and-cut data, we computed the Binder cumulant $`g(\mu ,N)`$, defined as in Eq.(11) but with the time average replaced by the sample average. According to finite-size scaling, the curves for $`g(\mu ,N)`$ as a function of $`\mu `$ for various $`N`$ must cross at the critical point $`\mu =\mu _c`$. In Fig. 6 we plot the Binder cumulant in the vicinity of the intersection point (note that the horizontal scale is much larger than that of Fig. 5), from which we obtain

$`\mu _c`$ $`=`$ $`0.77\pm 0.02(z=4)`$
$`\mu _c`$ $`=`$ $`0.56\pm 0.02(z=6).`$
This agrees with $`\mu _c^{BP}`$ within the error bars, suggesting that also for the magnetization replica symmetry breaking corrections are small, causing a shift of $`\mu _c`$ of less than $`34\%`$. Replica symmetry breaking corrections are expected to increase with $`z`$. In the Sherrington-Kirkpatrick model (which is the $`z\mathrm{}`$ limit of the present model), corrections shift $`\mu _c`$ from 1.25 to 1. Although our numerical estimate of $`\mu _c`$ is slightly larger than $`\mu _c^{BP}`$ instead, this could be a statistical fluctuation or a finite-size effect. The small size of replica-symmetry-breaking corrections to $`\mu _c`$ suggests that the mixed ferromagnetic spin-glass phase is narrow for these values of $`z`$, as also recently indicated for the three-dimensional Ising spin glassmartin\_mixed .

The Binder cumulant is expected to satisfy the following finite-size scaling relationprivman  for $`\mu \mu _c`$:

$$g(\mu ,N)=\stackrel{~}{g}(N^{1/(d_u\nu )}(\mu \mu _c))$$
(12)
where $`d_u`$ is the upper critical dimension, which for the Ising spin glass is $`d_u=6`$. As usual, by plotting $`g(\mu ,N)`$ against $`N^{1/(d_u\nu )}(\mu \mu _c)`$ with correct parameters $`\mu _c`$ and $`\nu `$, the data points for different system sizes should collapse onto a single curve near $`(\mu \mu _c)=0`$. As shown in Fig. 7, using the estimates of $`\mu _c`$ obtained above and the mean-field exponent $`\nu =1/2`$ we obtain a good data collapse, showing that finite size scaling is well satisfied in our range of sizes.

In Fig. 8 we also show scaling plots for the average magnetization $`m(\mu ,N)=[M]_J`$, whose scaling form is

$$m(\mu ,N)=N^{\beta /(d_u\nu )}\stackrel{~}{m}(N^{1/(d_u\nu )}(\mu \mu _c)),$$
(13)
with the mean field exponent $`\beta =1/2`$. The data show a good scaling collapse for $`\mu \mu _c`$.

## VI Typical running time of our branch and cut algorithm

In this section we study the running time of our program as a function of the mean coupling strength $`\mu `$. In computer science the complexity of a problem is classified in terms of the worst-case running time of its solution algorithmsGaJo ; mertens\_lecture . Central notions here are the complexity classes P and NP. Informally, the class P consists of all decision problems (namely, problems whose solution can only be “yes” or “no”) for which at least one algorithm is known that can generate an answer in polynomial time, even in the “worst case”. The class NP consists of all decision problems for which, if for a given instance the answer is “yes”, then there is a certificate from which the correctness of the answer can be verified in polynomial time. For example, the question “Given a spin glass instance, is there a spin configuration with energy less than or equal to $`E_0`$?” belongs to NP. If for a given instance the answer is positive, then there is a spin configuration with correct energy, and its correctness can be verified in polynomial time. Only the existence of such a certificate (spin configuration, in the above example) is required, not the ability to find it in polynomial time. The class NP contains P, but it might be larger (many believe it is larger, and answering the question whether P = NP is an important open problem). NP-complete problems are the “most difficult” in the class NP, in the sense that no polynomial algorithm is known for solving them, and if a polynomial algorithm could be found for one of them, this would imply that all of them are polynomially solvableACGKMP1999 . The classes P and NP are defined for decision problems, but similar ideas apply to combinatorial optimization problems as well. Informally, an optimization problem is called NP-hard, if it is at least as difficult as every NP-complete problem. In particular, an optimization problem is NP-hard if the associated decision problem is NP-complete. This is true for many optimization problems, e.g. the maximum cut problem or the travelling salesman problem.

In practice, the running time can vary greatly from an instance of the problem to another, and the worst-case running time might very rarely occur. Recent work has therefore focused on the average running time with respect to random instances drawn from some probability distribution. Instead of the average one can also analyze the median, or typical running time, which has the advantage of being less influenced by the occurrence of exponentially rare samples with huge running times.

It should be noted that, unlike the worst-case complexity classification discussed above, which is an algorithm-independent feature of the problem itself, in general the typical running time can be different for different algorithms and implementations that solve the same problem.

Returning to our problem, as mentioned in Section III, finding the ground state of the Bethe-lattice spin glass is an NP-hard problem. For all values of $`\mu `$ the possible realizations of the disorder are the same as for $`\mu =0`$. Hence, the algorithm has an exponential worst-case running time even on instances deep in the ferromagnetic phase. However, for large $`\mu `$ highly frustrated realizations are very unlikely to appear, hence the typical running time will decrease as $`\mu `$ increases. The question we ask here is whether, for large $`N`$, the running time undergoes a sharp transition as a function of $`\mu `$ and, if so, whether the transition coincides with the spin-glass/ferromagnet phase transition.

One may use the CPU time as a measure of the running time. However, the CPU time is machine-dependent, hence it is not suitable when different computers are used. Furthermore, it is hard to separate out the influence of size-dependent hardware effects on the CPU time (for example, small problems can be fully stored in the cache and therefore run faster). To avoid these problems, in the following we use the number of linear problems solved, $`n_{lps}`$, as a measure of the running timebulk\_perturb .

In Fig. 9 we show the median running time so defined as a function of $`\mu `$ for $`z=4,6`$ and different system sizes. Clearly, ground states are calculated quickly in the ferromagnetic region, while in the spin-glass phase the running time increases dramatically (note the logarithmic scale on the vertical axis), and is approximately constant within the entire spin glass phase. The variation becomes more pronounced as $`N`$ increases, suggesting a sharp discontinuity in the $`N\mathrm{}`$ limit around $`\mu 0.8`$ ($`z=4`$) and $`\mu 0.6`$ ($`z=6`$), which is close to the spin-glass/ferromagnet transition point $`\mu _c`$ determined in Section V.

As shown in Fig. 10, deep in the ferromagnetic phase the data is consistent with a polynomial increase of the running times with $`N`$. For smaller values of $`\mu `$, the curves are bending upwards, indicating that the running time increases faster than any polynomial. This is also the case for $`\mu =0.8`$ ($`z=4`$) and for $`\mu =0.6`$ ($`z=6`$, not shown). Hence from this data, it seems that the change in the typical running time occurs at a value of $`\mu `$ larger than $`\mu _c`$, although it is difficult to locate a precise transition point. A mismatch between phase transition and change of the running time has been observed before, e.g. for a simple algorithm solving vertex covercover-time .

We have fitted the data in Fig. 10 with a function of the form $`n_{lps}(N)\mathrm{exp}(bN^c)`$. For $`\mu =0`$, we find $`b=0.026(9)`$, $`c=0.87(5)`$ for $`z=4`$, and $`b=0.007(3)`$ and $`c=1.24(8)`$ for $`z=6`$, but the data exhibits in both cases a considerable scatter around the fitting region, prohibiting to conclude in a definite way that the typical running time is exponential. Nevertheless, the data strongly suggest so.

## VII Conclusions

We have studied the ground state of a diluted mean-field Ising spin glass model with fixed connectivities $`z=4,6`$ and Gaussian distribution of the couplings, with mean $`\mu `$ and unit variance. We have applied a branch-and-cut algorithm, a sophisticated technique originating in combinatorial optimization which guarantees to find exact ground states. Our motivation was to study the spin-glass/ferromagnet transition and relate it to the change in the typical running time of our algorithm.

From the study of the Binder cumulant, we have obtained values for the critical coupling strength, $`\mu _c`$. We have also solved the model in the Bethe-Peierls approximation, using an iterative stochastic procedure. In this approximation we obtain a critical coupling strength, $`\mu _c^{BP}`$, which agrees with the branch-and-cut estimate within the error bars of the latter, indicating that replica symmetry breaking effects are quantitatively small. Finite-size scaling is well satisfied for systems of size larger than $`N30`$.

We have also analyzed the ground state energy, and shown that the branch-and-cut results, extrapolated to the thermodynamic limit, are in very good agreement with the Bethe-Peierls results, again indicating that replica symmetry breaking effects are quantitatively small. In the spin glass region, finite-size corrections are well described by a $`N^{2/3}`$ dependence.

We have investigated the typical running time of our implementation of the branch-and-cut algorithm, which we defined as the median number of linear programs needed to find the ground state, with respect to a uniform distribution over the space of instances. We have shown that finding ground states is “hard” in the spin-glass phase, and “easy” deep in the ferromagnetic region, with a sharp variation at a value of $`\mu `$ slightly larger than $`\mu _c`$. The data indicate that while the worst-case running time is always exponential in the system size, the typical running-time is polynomial in the ferromagnetic phase and super-polynomial in the spin-glass phase.

Our understanding of what makes a problem computationally hard is still very weak. In this paper, we have shown that in a standard hard problem from physics, the Ising spin glass, a “physical” phase transition has a dramatic effect on the performance of a solution algorithm. Although in principle the “typical hardness” is algorithm-dependent, it is reasonable to expect that the phase transition will influence to some extent the running time of many other solution algorithms. Furthermore, we expect that similar phenomena occur in other well-known physical models.

###### Acknowledgements.

We thank A.P. Young for comments on an earlier version of the paper, and the Regional Centre of Computing in Cologne for the allocation of computer time and various support. MP also thanks Alan Bray for a useful correspondence. AKH obtained financial support from the DFG (Deutsche Forschungsgemeinschaft) under grants Ha 3169/1-1 and Zi 209/6-1.


format: text

text:
# A simple plan for the Δ resonance

## I Introduction

The $`\mathrm{\Delta }`$ resonance is present in most of baryonic physics, due to its low energy and its strong coupling to the pion-nucleon systems. PJ90 ; SL96 ; HDT96 ; P01 ; VDL ; Arndt ; Surya  Therefore it is is crucial to have good and simple methods to handle the $`\mathrm{\Delta }`$ resonance and its effects on other hadronic physics. Such method should be based on a local theory where the renormalization is handled properly, and divergences can be re-absorbed consistently in the interactions. Otherwise, the method is nothing more than a particular fit of the data, where the $`\pi N\mathrm{\Delta }`$ interactions have no universal meaning. Especially this interaction should be handled properly, due to strength, which makes the application of perturbation theory dubious. In this paper we propose such a method with local interactions, which can be solved in closed form, without resorting to perturbation theory.

The $`\mathrm{\Delta }`$ resonance is a spin 3/2 state, which is hard to handle in field theory, as it is generally part of a non-renormalizable field theory, due to the high-order derivatives in the coupling and the propagator, which yield high-order divergences in the self-energy corrections. Furthermore, there are at least two other issues that blur our image of the $`\mathrm{\Delta }`$ resonance. First, the propagator has some unknown off-shell dependence, reflected by a free parameter. Second, field transformations allow different representations for the interactions and the propagators which, however, lead to identical observable results. KLF00

In this short paper all the problems raised above are circumvented. The propagator and the interaction are lumped together, so there are no issues concerning field transformations that shift terms around. Furthermore, dimensional analysis, combined with local field theory, implies that the interaction depends, in lowest order, on four parameters only: two $`\pi N\mathrm{\Delta }`$ couplings, one $`\pi \pi NN`$ couplings, and the $`\mathrm{\Delta }`$ mass. The finite renormalization constants are fixed by the $`\mathrm{\Delta }`$ mass and the wave function renormalization.

## II Theory

The self-energy diagram Kor97  of the $`\pi N`$ loop correction to the $`\mathrm{\Delta }`$, which has an imaginary part because of the $`\pi N`$ decay, has only two types of local momentum dependence of the loop integral in the rest frame:

$$\mathrm{\Sigma }(E)=\frac{1}{4\pi }\frac{d^3𝐩}{\omega _N\omega _\pi }\frac{[1;𝐩^2]}{E\omega _\pi \omega _N},$$
(1)
due to the different terms in the traces over polarizations and spins. The attention focusses on the real part. The same diagram occurs in other spin-isospin channels, with, possibly, different masses and coupling constants, however, the momentum and energy dependence will be the same. All the singularities are treated a principal value integrals. The imaginary part, associated with the homogeneous solution of the wave equation will follow from the eigenvalue equation below. Fan61  The two different momentum dependencies are associated with two different vertices (see Fig. 1). They serve as the two $`\pi N\mathrm{\Delta }`$ couplings. Furthermore, there is a third coupling for the $`\pi \pi NN`$ contact interaction, which will yield the phase shift in the absence of the $`\mathrm{\Delta }`$ resonance in other spin-isospin channels. In lowest order, the contact interaction is dominated by the pionic current $`[\pi ^{},_\mu \pi ]`$ which acts on the side of the triangle in Figure 1. The two different time-orderings will be associated with two different terms in the Hamiltonian.

For each channel and partial wave, in the rest frame, we can write down a general Hamiltonian for low-energy $`\pi N`$ scatttering based on the observations above:

$`H`$ $`=`$ $`|\mathrm{\Delta }m_\mathrm{\Delta }\mathrm{\Delta }|+{\displaystyle _{m_\pi +m_N}^{\mathrm{}}}𝑑ϵ|\pi N(ϵ)ϵ\pi N(ϵ)|`$ (2)
$`+`$ $`{\displaystyle \underset{n=0}{\overset{1}{}}}ϵ_{nk}{\displaystyle _{m_\pi +m_N}^{\mathrm{}}}𝑑ϵ^{}\sqrt{V_k(ϵ^{})}|\pi N(ϵ^{}){\displaystyle _{m_\pi +m_N}^{\mathrm{}}}𝑑ϵ\pi N(ϵ)|\sqrt{V_n(ϵ)}`$
$`+`$ $`{\displaystyle \underset{n=0}{\overset{1}{}}}{\displaystyle _{m_\pi +m_N}^{\mathrm{}}}𝑑ϵ|\mathrm{\Delta }\sqrt{W_n(ϵ)}\pi N(ϵ)|+|\pi N(ϵ)\sqrt{W_n(ϵ)}\mathrm{\Delta }|,`$
where $`ϵ_{nk}=1`$ if $`nk`$ and zero otherwise. The momentum dependence of the contact interaction is based on a single $`\pi \pi NN`$ coupling through the pionic current $`[\pi ,_\mu \pi ]`$ with the nucleon current $`\overline{N}\gamma ^5\gamma _\mu N`$, where the single derivative acts either on the incoming or the outgoing pion state, which leads to the asymmetric form.

The coupling functions $`V`$ and $`W`$, which correspond to local interactions, contain both two terms with the same functional form, however, different coupling constants ($`V_n=\frac{\lambda }{g_n}W_n`$):

$$W_n(ϵ)=g_n\frac{k^{2n+2}}{\omega _N\omega _\pi }\frac{k}{ϵ}=\frac{g_n}{ϵ}\left(\frac{\sqrt{ϵ^42ϵ^2(m_\pi ^2+m_N^2)+(m_\pi ^2m_N^2)^2}}{2ϵ}\right)^{2n+1},$$
(3)
The kinematical factors $`(\omega _\pi \omega _N)^1`$ make the interaction relativistic invariant, and the factor $`k^2`$ for $`n=1`$ is due to the local, derivative coupling, in the rest frame. The coupling constants $`g_n`$ have a dimension $`[\mathrm{energy}^{22n}]`$(GeV), and $`\lambda `$ has a dimension of $`[\mathrm{energy}^{12n}]`$(GeV). All the numerical factors from angular integrations, etc. are lumped in with the coupling constants $`g_n`$ and $`\lambda `$. The self-energy diagram, Eq. (1), is equivalent to ($`g_n=1`$):

$$\mathrm{\Sigma }(E)=_{m_\pi +m_N}^{\mathrm{}}𝑑ϵ\frac{[W_0(ϵ);W_1(ϵ)]}{Eϵ},$$
(4)
which is again equivalent with the first order perturbation theory results for the Hamiltonian Eq. (2).

This single channel problem with separable, but local interactions can be solved exactly, apart from the divergencies for which the finite renormalization has to be determined. The coupling functions $`W`$ and $`V`$ have the same functional form, therefore the ansatz wave function for the energy $`\omega `$ yields:

$`|\mathrm{\Psi }(\omega )`$ $`=`$ $`\alpha (\omega )|\mathrm{\Delta }`$
$`+`$ $`{\displaystyle 𝑑ϵ\left(\frac{1}{\omega ϵ}z(\omega )\delta (\omega ϵ)\right)(\sqrt{W_0(ϵ)}\beta _0(\omega )+\sqrt{W_1(ϵ)}\beta _1(\omega ))|\pi N(ϵ)}`$
$``$ $`\alpha (\omega )|\mathrm{\Delta }+\beta _0(\omega )|\omega (\pi N)_0+\beta _1(\omega )|\omega (\pi N)_1.`$ (6)
The function $`z(\omega )`$ models the ratio of the homogeneous and the inhomogeneous part of the $`\pi N`$ scattering state. For this single channel scattering state it is the only observable, which is directly related to the phase shift. When this wave function is inserted into the continuum eigenvalue equation $`0=(\omega H)|\mathrm{\Psi }(\omega )`$ it yields the relations between the spectroscopic densities $`\beta _i`$ and the occupation number $`\alpha `$:

$`0`$ $`=`$ $`\beta _0|\omega (\pi N)_0+\beta _1|\omega (\pi N)_1+(\omega m_\mathrm{\Delta })\alpha |\mathrm{\Delta }`$ (7)
$``$ $`(\stackrel{~}{W}_0(\omega )zW_0(\omega ))\beta _0|\mathrm{\Delta }(\stackrel{~}{W}_1(\omega )zW_1(\omega ))\beta _1|\mathrm{\Delta }`$
$``$ $`\alpha (|\omega (\pi N)_0+|\omega (\pi N)_1)`$
$``$ $`(\stackrel{~}{V}_0(\omega )zV_0(\omega ))\beta _0|\omega (\pi N)_1(\stackrel{~}{V}_1(\omega )zV_1(\omega ))\beta _1|\omega (\pi N)_0,`$
where the tilde refers to the Hilbert transform, without a factor of $`\pi ^1`$:

$$\frac{g_n}{\lambda }\stackrel{~}{V}_n(\omega )=\stackrel{~}{W}_n(\omega )=𝑑ϵ\frac{W_n(ϵ)}{\omega ϵ}.$$
(8)
Projecting on the states $`|\omega (\pi N)_0`$ and $`|\omega (\pi N)_1`$ yields: $`(i,j=0,1)`$

$$\beta _i=ϵ_{ij}\frac{1+\stackrel{~}{V}_j(\omega )zV_j(\omega )}{1[\stackrel{~}{V}_0(\omega )zV_0(\omega )][\stackrel{~}{V}_1(\omega )zV_1(\omega )]}\alpha ,$$
(9)
inserting this back into Eq. (7) and projecting onto $`|\mathrm{\Delta }`$ yields a relation for $`z`$:

$$\omega m_\mathrm{\Delta }=\frac{_{i=0}^1ϵ_{ij}[\stackrel{~}{W}_i(\omega )zW_i(\omega )][1+\stackrel{~}{V}_j(\omega )zV_j(\omega )]}{1[\stackrel{~}{V}_0(\omega )zV_0(\omega )][\stackrel{~}{V}_1(\omega )zV_1(\omega )]},$$
(10)
which yields a straightforward but lengthy quadratic equation for $`z`$. In the case of weak contact interactions $`\lambda 0`$, the solution can be approximated by:

$$z=\frac{1}{W_0+W_1}(\omega m_\mathrm{\Delta }\stackrel{~}{W}_0\stackrel{~}{W}_1),$$
(11)
which can be cast in the form of an analyticity corrected energy-dependent Breit-Wigner resonance.

The general phase shift is given by:

$$\mathrm{tan}\delta =\frac{\pi }{z}.$$
(12)
Since there is only one channel, the scattering process is fully characterized by $`z`$ through the phase shift. Only the Hilbert transforms of $`W_0=\frac{g_0}{\lambda }V_0`$ and $`W_1=\frac{g_1}{\lambda }V_1`$ still need to be determined.

The backward or $`Z`$-diagram contributions are included, to make the results covariant and the expressions simpler. The singularity is treated as a principal value singularity:

$$\stackrel{~}{W}_n(\omega )=g_n\frac{ds}{\sqrt{s}}\left(\frac{\sqrt{s^22s(m_\pi ^2+m_N^2)+(m_\pi ^2m_N^2)^2}}{2\sqrt{s}}\right)^{2n+1}\frac{1}{\omega ^2s}.$$
(13)
The integral is divergent. Subtractions make the integrals finite:

$$\stackrel{~}{W}_n^\mathrm{r}(\omega )=\frac{\omega ^{2n+2}}{s^{n+1}}𝑑sW_n(\sqrt{s})\frac{1}{\omega ^2s},$$
(14)
which yield an analytical results for the Hilbert transform for an arbitrary integer value of $`n`$:

$$\stackrel{~}{W}_n^\mathrm{r}(\omega )=\frac{g_n(1𝒯_{2n+1})}{\omega (2\omega )^{2n+1}}\mathrm{Re}\left[\left(\sqrt{A}\right)^{2n+1}\mathrm{log}\left(\frac{\omega ^2+m_\pi ^2+m_N^2+\sqrt{A}}{m_\pi m_N}\right)\right],$$
(15)
where $`𝒯_m`$ stands for the $`m`$-th order Taylor expansion in $`\omega ^2`$ around $`\omega ^2=0`$, and

$$\sqrt{A}\sqrt{\omega ^42\omega ^2(m_\pi ^2+m_N^2)+(m_\pi ^2m_N^2)^2}.$$
(16)
It should noted that the real part of the function above arises through the real part of the logarithm and the real part of the squareroot for $`\omega ^2>(m_N+m_\pi )^2`$ and $`\omega ^2<(m_Nm_\pi )^2`$ and through the product of the imaginary parts of both for $`(m_Nm_\pi )^2<\omega ^2<(m_N+m_\pi )^2`$.

The subtractions correspond to two divergent constants that should be reabsorbed in the mass and the wave function of the bare $`\mathrm{\Delta }`$ state:

$`\stackrel{~}{W}_0(\omega )`$ $`=`$ $`\stackrel{~}{W}_0^\mathrm{r}(\omega )+c_0^{(0)},`$ (17)
$`\stackrel{~}{W}_1(\omega )`$ $`=`$ $`\stackrel{~}{W}_1^\mathrm{r}(\omega )+c_0^{(1)}+c_1^{(1)}\omega ^2,`$ (18)
where the renormalization constants $`c`$ are fixed through: $`\stackrel{~}{W}_n(m_\mathrm{\Delta })=(\omega ^2m_\mathrm{\Delta }^2)^n+𝒪(n+1)`$. Arguably, $`W_0`$ could have a finite wave function renormalization. However, as one can see from Eq. (11) this renormalization factor $`c_1^{(0)}`$ can be reabsorbed in an overall shift of the coupling constants $`g_ng_n/c`$, for the only observable quantity $`z`$, where $`\omega ^2m_\mathrm{\Delta }^2`$ is approximated by $`2m_\mathrm{\Delta }(\omega m_\mathrm{\Delta })`$. The coupling functions $`W`$ can clearly be identified with the bare $`\mathrm{\Delta }`$ state; the constant part of $`\stackrel{~}{W}`$ is indistinquishable from the mass $`m_\mathrm{\Delta }`$. For the contact interaction there is no such identification, therefore these renormalization constants will be determined from the data.

## III Without bare resonance state

In the case of only contact interactions and no resonance state the expressions simplify somewhat. It corresponds to the Hamiltonian, Eq. (2), with $`W_0=0`$ and $`W_1=0`$, which decouples the resonance from the scattering states. The ratio of the homogeneous to inhomogeneous solution $`\pi /z(\omega )`$ that leads to the phase shift is given by:

$$z=\frac{V_0\stackrel{~}{V}_1+V_1\stackrel{~}{V}_0+\sqrt{(V_0\stackrel{~}{V}_1+V_1\stackrel{~}{V}_0)^2+4V_0V_1(1\stackrel{~}{V}_0\stackrel{~}{V}_1)}}{2V_0V_1}.$$
(19)
Depending on the coupling constant $`\lambda `$, the phase shift can be positive or negative. In the physical range $`m_N+m_\pi <\omega <2.0`$ GeV the sign change in the phase shift occurs for $`\lambda `$ varying between 1.0 and 3.0. In that case the finite renormalization is set to zero, such that $`V_i`$ vanishes at $`\omega =0`$.

It is possible to mimic part of the resonance with a mass $`m_V`$ through a non-trivial renormalization of $`\stackrel{~}{V}_0`$ and $`\stackrel{~}{V}_1`$:

$`\stackrel{~}{V}_0`$ $`=`$ $`\stackrel{~}{V}_0^\mathrm{r}+c_0^{(\lambda ,0)}`$ (20)
$`\stackrel{~}{V}_1`$ $`=`$ $`\stackrel{~}{V}_1^\mathrm{r}+c_0^{(\lambda ,1)}+c_1^{(\lambda ,1)}\omega ^2,`$ (21)
where $`m_V^2=(\stackrel{~}{V}_1^\mathrm{r}(m_V)+c_0^{(\lambda ,1)})/(c_1^{(\lambda ,1)}\omega ^2+_{\omega ^2}\stackrel{~}{V}_1^\mathrm{r}(m_V))`$. This is normally referred to as a dynamically generated resonance.

## IV Results

In order to test the quality of the model a whole set of fits were performed, excluding certain parts of the interaction, and including others. The least-square fit is determined by minimizing the $`\chi `$-square value. However, since there is a definite non-linear dependence of the scattering data on the coupling constants some care is required. As expected already above, fitting with both the contact interaction and the bare resonance is indeterminate, as both yield a similar effect and a flat minimum in parameter space for $`\chi `$-square.

However, as it turns out all the possible fits for the $`\mathrm{\Delta }`$ resonance in the $`P_{33}`$ channel do much worse that a fit with a quasi-local interaction $`n=\frac{1}{2}`$. This interaction has the wrong dimensionality, however, produces the best fit with two parameters: the coupling constant and the renormalized resonance mass. This is also the only way to yield the right threshold behavior, if the fit includes data points beyond the resonance.

The various results are listed below and shown in Fig. 2. In the cases that the bare resonance is presence, i.e., $`g_00`$ or $`g_10`$, its mass is set to zero, i.e., $`m_\mathrm{\Delta }=0`$, the renormalization constants $`c_0`$ will account for the mass, with a finite renormalization. The energy for which $`z=0`$ and the phase shift crosses $`90^{}`$ is for every fit within a few MeV of the $`\mathrm{\Delta }`$ resonance at $`1.232`$ GeV. The number of data points for the fits ranges between 30 and 50 depending on the range.

* For the range 1.1 GeV to 1.6 GeV, fitting only with nonzero $`g_0`$ yields $`\chi ^2/\mathrm{\#}=120`$ and $`g_0=0.082`$ and $`c_0^{(0)}=1.294`$.
* For the range 1.1 GeV to 1.6 GeV, fitting only with nonzero $`g_1`$ yields $`\chi ^2/\mathrm{\#}=208`$ and $`g_1=76000.`$ and $`c_0^{(1)}=10100.`$. The results do not depend significantly on $`c_1^{(1)}`$, which is subsequently set to zero. However, if we set $`c_1^{(1)}=1`$, we find more appropriate numbers for $`g_1=3.78`$ and $`c_0^{(1)}=0.493`$, since for this particular energy $`_\omega W_1^\mathrm{r}(m_\mathrm{\Delta })`$ gives a finite wave-function renormalization near unity, which makes the energy dependence of $`z`$ very weak.
* For the range 1.1 GeV to 1.6 GeV, fitting only with nonzero $`\lambda `$ yields $`\chi ^2/\mathrm{\#}=140`$, $`\lambda =2.31`$, $`c_0^{(\lambda ,1)}=2.26`$, $`c_1^{(\lambda ,1)}=1.47`$, and $`c_0^{(\lambda ,0)}=8.69`$. The fit yields a better threshold behavior but a worse high-energy behavior compared to the $`g_0+g_1`$ fit.
* For the range 1.1 GeV to 1.6 GeV, fitting with nonzero $`g_0`$ and $`g_1`$ yields $`\chi ^2/\mathrm{\#}=82`$, $`g_0=0.060`$, $`g_1=0.79`$, $`c_1^{(0)}=0.69`$, and $`c_1^{(1)}=0.69`$. Both renormalization constants take about half the resonance mass, changing this fraction does not change the $`\chi ^2`$ value. Changing the range of the data set changes the results substantially, as can be seen in Figure 3. At threshold till 1.25 GeV the $`g_1`$ coupling fits the data well.
* For the range 1.1 GeV to 1.6 GeV, fitting only with nonzero $`g_{\frac{1}{2}}`$ yields $`\chi ^2/\mathrm{\#}=56`$, $`g_{\frac{1}{2}}=0.44`$, and $`c_0^{(\frac{1}{2})}=1.227`$. It is by far the best fit archieved.
* For the extended range 1.1 GeV to 1.8 GeV, the fit even improves, due to poorer quality data: $`\chi ^2/\mathrm{\#}=45`$ and $`g_{\frac{1}{2}}=0.44`$, and $`c_0^{(\frac{1}{2})}=1.227`$. The actual values of the coupling constant and the renormalization constant do not change. For the whole range the two parameters $`g_{\frac{1}{2}}`$ and $`m_\mathrm{\Delta }`$ fit the data with a consistent coupling strength of $`g_{\frac{1}{2}}=47\pm 2`$, as can be seen in Figure 4. Fits with data below the inelastic threshold the $`\chi ^2/\mathrm{\#}`$ ranges between 16 and 23.
* A further dramatic improvement is archieved with an additional resonance. For the range 1.1 GeV to 2.0 GeV, the fit yields: $`\chi ^2/\mathrm{\#}=13`$, $`g_{a\frac{1}{2}}=0.48`$, $`c_0^{(a,\frac{1}{2})}=1.231`$, $`g_{b\frac{1}{2}}=1.75`$, and $`c_0^{(b,\frac{1}{2})}=3.74`$. The second resonance is around 2.5-4.0 GeV, however, its position cannot be determined accurately from the data, since it lies outside the range.

Although, consistency with non-relativistic treatment, where the power $`n`$ in $`W_n`$ is directly related to the partial wave, of the $`\mathrm{\Delta }`$ resonance would require that $`g_0`$ is zero, it is impossible to fit the $`P_{33}`$ data properly with that assumption. Attempts to make sensible fits with $`g_1`$ and two bare states failed. The accurate data at low energies forces a $`g_{\frac{1}{2}}`$ interaction in the case of two bare states. However, the $`P_{31}`$ and the $`P_{13}`$ channels do give results consistent with $`g_0=0`$. A number of fits were made, which are listed below, and show respectively in Fig. 5 and Fig. 6.

* For the $`P_{31}`$ a fit with $`g_1`$ only gives excellent results from the threshold up to the first structure in the data at 1.55 GeV. The $`\chi `$-square is comparable to unity.
* Similarly for the $`P_{31}`$ channel, a fit with $`\lambda `$ only gives good results too. However, the second solution of the quadratic equation starts to interfere with the physical solution at 1.6 GeV.
* Combining $`g_1`$ and $`\lambda `$ interactions does not give a significant improvement with respect to either fit.
* For the $`P_{13}`$ the data is rather coarse. Some structure appears at 1.4 GeV. Both $`g_1`$ on its own, as $`\lambda `$ on its own give a fit with a $`\chi `$-square comparable to unity. The shape of the fit depends on the range of the data points, however, as the $`\chi `$-square is comparable to unity, no conclusions can be drawn.
* Combining both fits for the $`P_{13}`$ data leads to an improvement as some of the structure at 1.7 GeV is reproduced, however, as the threshold behavior fails, it is considered an artificial improvement.

Since in the case of the $`P_{31}`$ and the $`P_{13}`$ channels there is no resonance that fixed a scale, there is considerable freedom to scale the coupling constants and renormalization constants and still produce the same results. However, especially for a non-zero $`\lambda `$, there seems to be no obvious rebundant parameters in the model, wich can be fixed a priori. The different coupling constants and renormalization constants are not quoted as they, on their own, have lost their significance.

## V Conclusions

A local field theory without cut-offs or form factors should fit the data at all energies. Deviations can only occur if other decay channels open, or higher-order terms in the interaction are important. Therefore we expect the data to be fitted well with our model below the two-pion production threshold, which is the case for the $`P_{31}`$ and the $`P_{13}`$ channels. In the same theory resonance states can be pushed out of the physical energy region, and resonances can be generated dynamically through a contact interaction. Eventually it is the quality of the fit that determines which is the likely process. For the $`P_{33}`$ channel the bare state gives a better results, for the $`P_{31}`$ and the $`P_{13}`$ channels both the contact interaction and the bare resonance state whose pole is pushed away from the real axis through a strong interaaction do a good job fitting the data below 1.5 GeV.

From the results it is clear that one cannot simply add a resonance structure to a contact interaction incoherently. The effects of the contact interaction are not like what one would expect from perturbation theory. When the interaction is summed to all orders, the resonance dominate the results. Any other, perturbative, effects will be irrelevant for the phase shift around the resonance peak. In all the different versions of the model, where some interactions are switched off, it is impossible to create the effect of a rising phase shift beyond $`180^{}`$ after the $`\mathrm{\Delta }`$ resonance, usually associated with a background, without distorting the resonance itself. It can only consistently be explained by a second resonance, which is found around 2.5-4.0 GeV.

Furthermore, from the data the $`W_{\frac{1}{2}}`$ coupling function prevails over the $`W_0`$ and the $`W_1`$ coupling functions, or even a combination of both with five fitting parameters, compared to the two of the $`W_{\frac{1}{2}}`$ coupling alone.

The same model, with different values for the coupling constants and the renormalization constants, fit the low-energy $`P_{13}`$ and $`P_{31}`$ data well, which have no resonance structure. The fits give a model in with the derivative coupling dominates, associated with the $`g_1`$ coupling constant.


text:
# 1 Introduction

## 1 Introduction

Evidence of antineutrino disappearance in a beam of antineutrinos in the Kamland experiment has been recently presented . The analysis of these results  in terms of neutrino oscillations have largely improved our knowledge of neutrino mixing on the LMA region. The results appear to confirm in a independent way that the observed deficit of solar neutrinos is indeed due to neutrino oscillations. The ability to measure the LMA solution, the one preferred by the solar neutrino data at present, “in the lab” puts Kamland in a pioneering situation: after these results there should remain little doubt of the physical reality of neutrino mass and oscillations.

The publication of the SNO results  has already made an important breakthrough towards the solution of the long standing solar neutrino  problem (SNP) possible. These results provide the strongest evidence so far (at least until Kamland improves its statistics) for flavor oscillation in the neutral lepton sector. From the combined analysis of SNO and other solar evidence, one obtains, in the framework of two active neutrino oscillations, the following set of parameters, $`\mathrm{\Delta }m_{solar}^2=4.5_{1.4}^{+2.7}\times 10^5eV^2,\mathrm{tan}^2\theta _{solar}=0.40_{0.08}^{+0.10}.`$ We will see in this work how the evidence from the Kamland measurements improves or modifies these values.

The previous generation of reactor experiments (CHOOZ , PaloVerde ), performed with a baseline of about 1 km. They have attained a sensitivity of $`\mathrm{\Delta }m^2<10^3eV^2`$  and, not finding any disappearance of the initial flux, they demonstrated that the atmospheric neutrino anomaly  is not due to muon-electron neutrino oscillations. The Kamland experiment is the successor of such experiments at a much larger scale in terms of baseline distance and total incident flux. This experiment relies upon a 1 kton liquid scintillator detector located at the old, enlarged, Kamiokande site. It searches for the oscillation of antineutrinos emitted by several nuclear power plants in Japan. The nearby 16 (of a total of 51) nuclear power stations deliver a $`\overline{\nu }_e`$ flux of $`1.3\times 10^6cm^2s^1`$ for neutrino energies $`E_\nu >1.8`$ MeV at the detector position. About $`85\%`$ of this flux comes from reactors forming a well defined baseline of 139-344 km. Thus, the flight range is limited in spite of using several reactors, because of this fact the sensitivity of Kamland increases by nearly two orders of magnitude compared to previous reactor experiments.

The aim of this work is to study the implications of the recent Kamland results on the determination of the neutrino oscillation parameters, to understand which regions of the parameter space still allowed by the solar neutrino experiments are favored by them. The structure of this work is the following. In section 2 we discuss the main features of Kamland experiment that are relevant for our analysis: The salient aspects of the procedure we are adopting and the results of our analysis are presented and discussed in section 3. Finally, in section 4 we draw our conclusions and discuss possible future scenarios.

## 2 The computation of the expected signals

### 2.1 The Kamland signal

Electron antineutrinos from nuclear reactors with energies above 1.8 MeV are measured in Kamland by detecting the inverse $`\beta `$-decay reaction $`\overline{\nu }_e+pn+e^+`$. The time coincidence, the space correlation and the energy balance between the positron signal and the 2.2 MeV $`\gamma `$-ray produced by the capture of a already-thermalized neutron on a free proton make it possible to identify this reaction unambiguously, even in the presence of a rather large background.

The two principal ingredients in the calculation of the expected signal in Kamland are the reactor flux and the antineutrino cross section on protons. A number of short baseline experiments (See Ref. and references therein) have previously measured the energy spectrum of reactors at distances where oscillatory effects have been shown to be nonexistent. They have shown that the theoretical neutrino flux predictions are reliable within 2% . The effective flux of antineutrinos released by the nuclear plants is a rather well understood function of the thermal power of the reactor and the amount of thermal power emitted during the fission of a given nucleus, which gives the total amount, and the isotopic composition of the reactor fuel which gives the spectral shape. Detailed tables for these magnitudes can be found in Ref. . For a given isotope the energy spectrum can be parametrized by an exponential expression  where the coefficients depend on the nature of the fissionable isotope (see Ref. for explicit values). Along the year, between periods of refueling, the total effective flux changes with time as the fuel is expended and the isotope relative composition varies. We take the average of the relative fission yields over the live time as given by the experiment: $`{}_{}{}^{235}U=57\%`$, $`{}_{}{}^{238}U=7.8\%`$, $`{}_{}{}^{239}Pu=30\%`$, $`{}_{}{}^{241}Pu=5.7\%`$.

In order to obtain the expected number of events at Kamland, we sum the expectations for all the relevant reactor sources weighting each source by its power and distance to the detector (table II in Ref.  ), assuming the same spectrum originated from each reactor. We sum over the nearby power reactors, we neglect farther Japanese and Korean reactors and even farther rest-of-the-world reactors which give only a minor additional contribution. The average number of positrons $`N_i`$ which are detected per visible energy bin $`\mathrm{\Delta }E_i`$ is given by the convolution of different quantities: $`\overline{P}`$, the oscillation probability averaged over the distance and power of the different reactors. Expressions for the antineutrino capture cross section are taken from the literature . The matrix element for this cross section can be written in terms of the neutron half-life, we have used the latest published value $`t_{1/2}=613.9\pm 0.55`$ . The antineutrino flux spectrum, the relative reactor-reactor power normalization which is included in the definition of $`\overline{P}`$ and the energy resolution of Kamland are used in addition. We use in our analysis the following expression for the energy resolution in the prompt positron detection

$`\sigma (E_e)`$ $`=`$ $`0.0062+0.065\sqrt{E_e}.`$ (1)
This expression is obtained from the raw calibration data presented in Ref.. Note that we prefer to use this expression instead of the much less accurate one given in Ref.. Moreover, we assume a 408 ton fiducial mass and standard nuclear plant power and fuel schedule, we take an averaged, time-independent, fuel composition equal for each detector as given above. Detection efficiency is taken close 100% and independent of the energy .

We will not add any background events as they (a total of $`0.95\pm 0.99`$ events which include random coincidences from radioactive decays, $`<0.01`$ evts, and correlated background from cosmic ray muons and neutrons, $`<0.5`$ events,) can be distinguished from the signal with sufficiently high efficiency or are negligible above the 2.6 MeV analysis threshold (table I in Ref. ). We also consider negligible the background from geological neutrinos above the 2.6 MeV analysis threshold.

### 2.2 The Solar Signal

We also need the expected signals in the different solar neutrino experiments. These are obtained by convoluting solar neutrino fluxes, sun and earth oscillation probabilities, neutrino cross sections and detector energy response functions. We closely follow the same methods already well explained in previous works , we will mention here only a few aspects of this computation. We determine the neutrino oscillation probabilities using the standard methods found in literature , as explained in detail in  and in . We use a thoroughly numerical method to calculate the neutrino evolution equations in the presence of matter for all the parameter space. For the solar neutrino case the calculation is split in three steps, corresponding to the neutrino propagation inside the Sun, in the vacuum (where the propagation is computed analytically) and in the Earth. We average over the neutrino production point inside the Sun and we take the electron number density $`n_e`$ in the Sun by the BPB2001 model . The averaging over the annual variation of the orbit is also exactly performed. To take the Earth matter effects into account, we adopt a spherical model of the Earth density and chemical composition . The joining of the neutrino propagation in the three different regions is performed exactly using an evolution operator formalism . The final survival probabilities are obtained from the corresponding (non-pure) density matrices built from the evolution operators in each of these three regions.

In this analysis in addition to night probabilities we will need the partial night probabilities corresponding to the 6 zenith angle bin data presented by SK . They are obtained using appropriate weights which depend on the neutrino impact parameter and the sagitta distance from neutrino trajectory to the Earth’s center, for each detector’s geographical location.

## 3 Analysis and Results

In order to study power of the Kamland results for resolving the neutrino oscillation parameter space, we have developed two kind of analysis. In the first case we deal with the Kamland measured global signal. In the second case we include the full Kamland spectrum information. We perform a complete $`\chi ^2`$ statistical analysis before and after including in addition the up-date solar evidence obtaining the regions in the parameter space which are favored.

### 3.1 Analysis of the Global rate

The total $`\chi ^2`$ value is given by the sum of two distinct contributions, that is the one coming from all the solar neutrino data and the contribution of the Kamland experiment:

$$\chi ^2=\chi _{}^2+\chi _{glob,KL}^2,$$

with

$`\chi _{glob,KL}^2`$ $`=`$ $`\left({\displaystyle \frac{R^{exp}R^{teo}(\mathrm{\Delta }m^2,\theta )}{\sigma _{stat+sys}}}\right)^2.`$ (2)
The “experimental” signal ratio is $`R^{exp}=0.611\pm 0.085\pm 0.041`$ , where the first error is statistical and the second one systematic.

The solar neutrino contribution can be written in the following way:

$`\chi _{}^2`$ $`=`$ $`\chi _{\mathrm{glob}}^2+\chi _{\mathrm{SK}}^2+\chi _{\mathrm{SNO}}^2.`$ (3)
The function $`\chi _{\mathrm{glob}}^2`$ correspond to the total event rates measured at the Homestake experiment  and at the gallium experiments SAGE , GNO  and GALLEX . We follow closely the definition used in previous works (see Ref. for definitions and Table (1) in Ref. ‘ for an explicit list of results and other references).

The contribution to the $`\chi ^2`$ from the SuperKamiokande data ($`\chi _{\mathrm{SK}}^2`$) has been obtained by using double-binned data in energy and zenith angle (see table 2 in Ref. and also Ref.): 8 energy bins of variable width and 7 zenith angle bins which include the day bin and 6 night ones. The definition is given by:

$`\chi _{\mathrm{SK}}^2`$ $`=`$ $`(𝑹^{\mathrm{th}}𝑹^{\mathrm{exp}})^t\left(\sigma _{\mathrm{unc}}^2+\sigma _{\mathrm{cor}}^2\right)^1(𝑹^{\mathrm{th}}𝑹^{\mathrm{exp}}).`$ (4)
The theoretical and experimental $`𝑹`$ quantities are this time matrices of dimension 8$`\times `$7. The covariance quantity $`\sigma `$ is a 4-rank tensor constructed in terms of statistic errors, energy and zenith angle bin-correlated and uncorrelated uncertainties. The data and errors for individual energy bins for SK spectrum has been obtained from Ref. .

The contribution of SNO to the $`\chi ^2`$ is given by

$`\chi _{\mathrm{spec}\mathrm{SNO}}^2`$ $`=`$ $`{\displaystyle \underset{d,n}{}}(𝑹^{\mathrm{th}}𝑹^{\mathrm{exp}})^t\left(\sigma _{\mathrm{stat}}^2+\sigma _{\mathrm{syst}}^2\right)^1(𝑹^{\mathrm{th}}𝑹^{\mathrm{exp}}),`$ (5)
where the day and night $`𝑹`$ vectors of dimension 17 are made up by the values of the total (NC+CC+ES) SNO signal for the different bins of the spectrum. The statistical contribution to the covariance matrix, $`\sigma _{\mathrm{stat}}`$, is obtained directly from the SNO data. The part of the matrix related to the systematical errors has been computed by us studying the influence on the response function of the different sources of correlated and uncorrelated errors reported by SNO collaboration (see table II of Ref. ), we assume full correlation or full anticorrelation according to each source.

To test a particular oscillation hypothesis against the parameters of the best fit and obtain allowed regions in parameter space we perform a minimization of the two dimensional function $`\chi ^2(\mathrm{\Delta }m^2,\mathrm{tan}^2\theta )`$. A given point in the oscillation parameter space is allowed if the globally subtracted quantity fulfills the condition $`\mathrm{\Delta }\chi ^2=\chi ^2(\mathrm{\Delta }m^2,\theta )\chi _{\mathrm{min}}^2<\chi _n^2(CL)`$. Where $`\chi _{n=2}^2(90\%,95\%,\mathrm{})`$ are the quantiles for two degrees of freedom.

In Figs.1 we graphically show the results of this analysis. The favored region by the Kamland global rate alone is presented in fig.1(Left). Fig.1(right) includes the full Kamland and Solar contributions. Kamland global rate alone restricts the value of $`\mathrm{tan}^2\theta `$ to be in the range $`\mathrm{tan}^2\theta >4\times 10^11`$ (symmetric with respect the value $`\theta =\pi /4`$) except for some small region corresponding to $`\mathrm{tan}^2\theta 1`$ and $`\mathrm{\Delta }m^223\times 10^5`$ eV2. The overall effect of the measured observed ratio $`R0.6\%`$ is that the LMA region remains the only one which is still favored. Only if Kamland should had measured a ratio $`0.10.3`$, more than three sigma away from the actual measurement, other regions as the LOW solar solution would have some statistical chance to survive. In Table (I.a) we present the best fit parameters or local minima obtained from the minimization of the $`\chi ^2`$ function. Also shown are the values of $`\chi _{\mathrm{min}}^2`$. According to these results, regardless of the individual Kamland results, a maximal mixing solution is excluded when solar data is included in the analysis.

### 3.2 Analysis including the full Kamland signal

Here we use the binned Kamland signal (See table 2 extracted from Fig.5 in Ref.) combined with the evidence of the up-date solar experiments (CL,GA,SK and SNO). The total $`\chi ^2`$ value is given now by the sum of two distinct contributions, that is the one coming from all the solar neutrino data and the contribution of the Kamland experiment which includes both, shape and total signal information:

$`\chi ^2`$ $`=`$ $`\chi _{}^2+\chi _{spec,KL}^2+\chi _{glob,KL}^2.`$ (6)
The contributions of the solar neutrino experiments $`\chi _{}^2`$ and global Kamland signal are described in detail in the previous section. The contribution of the Kamland spectrum is now as follows:

$`\chi _{\mathrm{spec},\mathrm{KL}}^2`$ $`=`$ $`(\alpha 𝑹^{\mathrm{th}}𝑹^{\mathrm{exp}})^t\left(\sigma _{unc}^2+\sigma _{corr}^2\right)^1(\alpha 𝑹^{\mathrm{th}}𝑹^{\mathrm{exp}})`$ (7)
The total error matrix $`\sigma `$ is computed as a sum of assumed systematic deviations, $`\sigma _{sys}/S6.5\%`$, mainly coming from flux uncertainty ($`3\%`$), energy calibration and threshold (see table II of Ref.  for a total systematic error $`6.4\%`$), see also Ref.) and statistical errors. The parameter $`\alpha `$ is a free normalization parameter. The effect of systematic sources on individual bin deviations has been computed by us studying the influence on the response function, furtherly we have assumed full correlation among bins.

Note that, as an alternative to the previous formula 7 based on ratios and using and standard $`\chi ^2`$ expression, we could have directly used the absolute number of events in each bin and applied poissonian statistics formulas. The Kamland collaboration in its published analysis uses a mixed strategy and divides the $`\chi ^2`$ expression in two parts: it adds a standard $`\chi ^2`$ contribution corresponding to the global rate to a likelihood contribution which models the statistical significance of the spectrum shape. In this work, where we are interested not only on the Kamland data itself but in combining it with the solar data, we prefer to continue using the familiar expression 7 for various reasons: first, for simplicity, clarity and comparability with previous and future results, second, in the major part of the bins, certainly in the most significant ones, the signal is high enough, and more importantly the errors are small enough, to use the Gaussian approximation. For the large energy bins the signal is strictly zero: its contribution is anyway effectively zero in both approaches. Finally, not less importantly, the possibility of introducing in our approach in a straightforward way the effect of correlated systematic deviations among bins.

The $`𝑹`$ are length 13 vectors containing the binned spectrum (0.425 MeV bins ranging from 2.6 to 8.125 MeV) normalized to the non-oscillation expectations. Theoretical vectors are a function of the oscillation parameters: $`𝑹^{\mathrm{th}}=𝑹^{\mathrm{th}}(\mathrm{\Delta }m^2,\theta )`$. The experimental vectors $`𝑹^{\mathrm{exp}}`$ contain the Kamland binned signal. We generate acceptance contours (at 90,95 and 99 % CL) in the $`(\mathrm{\Delta }m^2,\mathrm{tan}^2\theta )`$ plane in a similar manner as explained in the previous section including now a minimization respec the parameter $`\alpha `$ for any of the other oscillation parameters. For the sake of comparison we have also obtained exclusion regions derived from the consideration of the Kamland contribution alone ($`\chi _{KL}^2=\chi _{spec}^2+\chi _{gl}^2`$).

In Figs.(2) we graphically show the results of this analysis. The first case, study of the Kamland data alone ($`\chi _{KL}^2`$) is represented by the Fig.(2)( left). The allowed regions in parameter space corresponding to each particular point is formed by a number of regions symmetric with respect the line $`\mathrm{tan}\theta =1`$. The general effect of the inclusion of the solar evidence, showed in Fig.2 (right), in the $`\chi ^2`$ is the breaking of the symmetry in $`\mathrm{tan}^2\theta `$, as expected, and the strong contraction of the allowed area to well defined regions. The best fit parameters obtained from the minimization of the $`\chi ^2`$ function are presented in Table (I.a) where two well separated solutions appear LMAI,LMAII. In both cases the mixing angle is around $`\mathrm{tan}^2\theta 0.450.48`$.

Again, the introduction of the solar data strongly diminishes the favored value for the mixing angle. The final value is more near to those values favored by the solar data alone than to the Kamland ones. This effect could be simply due to the present low Kamland statistics or, more worrying, to some statistical artifact derived from the complexity of the analysis and of the heterogeneity of binned data involved.

In this respect, we perform additionally a second kind of analysis in order to obtain concrete values for the individual oscillation parameters and estimates for their uncertainties. We study the marginalized parameter constraints where the $`\chi ^2`$ quantity is converted into likelihood using the expression $`=e^{(\chi ^2\chi _{min}^2)/2}`$. This normalized marginal likelihood, obtained from the integration of $``$ for each of the variables, is plotted in Figs. (3) for each of the oscillation parameters $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$. For $`\mathrm{tan}^2\theta `$ we observe that the likelihood function is strongly peaked in a region $`\mathrm{tan}^2\theta 0.40.6`$. The situation for $`\mathrm{\Delta }m^2`$ is similar except for the existence of additional, somehow narrower, secondary peaks. Concrete values for the parameters are extracted by fitting one- or two-sided Gaussian distributions to any of the peaks (fits not showed in the plots). In both cases, for angle and the mass difference distributions the goodness of fit of the Gaussian fit to each individual peak is excellent (g.o.f $`>99.8\%`$) thus justifying the consistency of the procedure. The values for the parameters obtained in this way appear in Table (I.b). They are fully consistent and very similar to the values obtained from simple $`\chi ^2`$ minimization. In particular, the maximal mixing solution is again excluded at the $`3\sigma `$ level.

Although both are mutually compatible, the slight difference of the value obtained for the mixing angle is well explained by the shape of the allowed regions in Fig 2 (right): the right elongation of these makes the value of the integral which defines the marginal distribution for $`\mathrm{tan}^2\theta `$ to be shifted. Additional variability can be easily introduced if would have used different prior information or different parameter definition (i.e. $`\mathrm{sin}2\theta `$ instead of $`\mathrm{tan}^2\theta `$). This an example of how the details of statistical analysis can significantly modify the values of the physical parameters extracted from data and how important is that different methods are explored, specially in the present context of neutrino physics where heterogeneous data coming from very different sources is jointly used.

## 4  Summary and Conclusions

We have analyzed the present experimental situation of our knowledge of the neutrino mixing parameters in the region of the parameter space that is relevant for solar neutrinos and we have studied in detail how this knowledge has improved with the recent results presented by the reactor experiment Kamland. We show, how, in general, the regions selected by Kamland alone, all symmetric with respect to $`\mathrm{tan}^2\theta =1`$, have a large spreadth in the mixing angle. The experiment has, however, a much higher sensitivity to the mass difference parameter.

Kamland global rate alone restricts the value of $`\mathrm{tan}^2\theta `$ to be in the range $`\mathrm{tan}^2\theta >4\times 10^11`$ (symmetric with respect the value $`\theta =\pi /4`$) except for some small region corresponding to $`\mathrm{tan}^2\theta 1`$ and $`\mathrm{\Delta }m^223\times 10^5`$ eV2. The overall effect of the measured observed ratio $`R0.6\%`$ is that the LMA region remains the only one which is still favored. Only if Kamland should had measured a ratio $`0.10.3`$, more than three sigma away from the actual measurement, other regions as the LOW solar solution would have some statistical chance to survive.

From the analysis of the Kamland evidence alone, the allowed regions in parameter space corresponding to each particular point is formed by a number regions symmetric with respect the line $`\mathrm{tan}\theta =1`$. The general effect of the inclusion of the solar evidence in the analysis is the breaking of the symmetry in $`\mathrm{tan}^2\theta `$, the strong reduction of the area of the allowed regions and, more importantly, the shifting of the value of the best fit mixing angle from being compatible with maximal mixing to more than $`3\sigma `$ incompatibility.

In addition to parameter extraction from $`\chi ^2`$ minimization, we obtain concrete values for the individual oscillation parameters and estimates for their uncertainties from the marginalized likelihood distributions. For $`\mathrm{tan}^2\theta `$ we observe that the likelihood function is concentrated in a region slightly higher that that one obtained from minimization. According to this analysis, maximal mixing is again excluded at the $`3\sigma `$ level. The marginal distribution for $`\mathrm{\Delta }m^2`$ shows clearly the existence of an additional, somehow narrower and much less significant, secondary peaks. Kamland, after only 150 days of data taking, together with the rest of solar experiments is already able to resolve the neutrino mass difference with a high precision of $`\delta \mathrm{\Delta }m_{kl+solar}^2<\pm 0.8`$ to be compared with the solar only case where the precision obtained was $`\delta \mathrm{\Delta }m_{solar}^2<\pm 1.5`$. The situation for the other oscillation parameter is worse, from the present data we obtain $`\delta \mathrm{tan}^2\theta _{kl+solar}<\pm 0.15`$. This precision is comparable, in fact even worse, than the precision attainable from the solar data alone $`\delta \mathrm{tan}^2\theta _{solar}<\pm 0.10`$ . The perspectives of Kamland, after 1-3 years of data taking, together or not with solar evidence, for a much better determination of the mixing angle are not very optimistic (see the results in Ref.). The ability of experiments as BOREXINO to improve the determination of the mixing angle with a smaller error in the near future remains an open question (see Ref.).

In summary, in the framework of two active neutrino oscillations we obtain, combining Kamland and Solar data and assuming CPT invariance, i.e. the same mass matrix for neutrino and antineutrinos, the following values for neutrino mixing parameters $`\mathrm{\Delta }m_{kl}^2=7.1\times 10^5eV^2,\mathrm{tan}^2\theta _{kl}=0.47`$ from $`\chi ^2`$ minimization (LMAI solution).

### Acknowledgments

We acknowledge the financial support of the Italian MURST and the Spanish CYCIT funding agencies. The numerical calculations have been performed in the computer farm of the Milano University theoretical group.