text:
# Fine-Tuning Solution for Hybrid Inflation in Dissipative Chaotic Dynamics

## I Introduction

Among the various proposed models for implementing inflation, the hybrid inflation model  is one of the most attractive ones due to the possibility of its implementation in the context of supersymmetry and supergravity models . The principle behind these models for inflation are basically based on models where the inflaton is coupled to one or more scalar fields. The inflationary phase is characterized by an initial phase in the evolution of the fields where the inflaton field slowly moves towards zero vacuum expectation value, till it approaches a critical value which then induces a spontaneous symmetry breaking in one of the other fields which are coupled to the inflaton, after which the fields quickly evolve to their vacuum states and inflation ends.

It has been shown in the context of supergravity motivated models  that for inflaton field amplitudes larger than the Planck mass there will appear large quantum corrections in the inflaton’s effective potential, destroying the flatness of the inflaton’s potential as required by the slow-roll conditions for inflation. The question then turns to whether inflation can be achieved in these models for initial inflaton field amplitudes smaller than the Planck scale, where a reliable model could be constructed. But in this case it has been shown  that small fluctuations of the other fields coupled to the inflaton are efficiently able to prevent the onset of inflation, or quickly drive previous inflationary trajectories to non-inflationary trajectories and therefore making inflation to end before enough e-folds of inflation have been produced ($`N_e70`$ as required) to solve the usual cosmological horizon and flatness problems. This is related to the homogeneity requirement over Hubble size distances  in order to inflation to begin. In the context of hybrid inflation this problem is even more severe, requiring an extreme fine-tuning of the initial conditions in order to have sufficient homogeneity, with negligible spatial and time derivatives of the other fields that are coupled to the inflaton, over regions extending far beyond a few Hubble lengths.

A few solutions to this homogeneity, or fine-tuning problem, have been proposed. The authors in Ref.  for example propose two-stage inflation models where a first short inflationary phase would smooth a large enough region in phase space making then possible a later longer second inflationary phase. In Ref.  it is analyzed the possibility of the solution of the fine-tuning problem in the context of more complicate generalizations of the hybrid inflation model as motivated from higher-dimensional scalar field models or brane cosmology.

Here I will be interested in investigating the possible role of how effective field interactions may shed further light in this fine-tuning problem in the pre-inflationary phase. This is mostly motivated from several works on the effective dynamics of scalar field configurations, which have shown that dissipation is intrinsically present in the field dynamics . In special, the role of dissipation in all of these works has been emphasized. We can trace the origin of these effects by considering the dynamics of a given system, which interacts with a sufficiently large environment (in the sense of degrees of freedom). One well known example of this, in the context of statistical mechanics, is the model of an oscillator, taken as being the system, in interaction with a large number of harmonic oscillators which are taken as representing the environment. By functionally integrating over the environment degrees of freedom we can then show that the dynamics of the system oscillator will be dissipative, with an equation of motion of the form of a Langevin-like one .

The dynamic of interacting fields has being a subject of intense study, in both Minkowski space-time , as also in FRW space-times , with special interest in studying the dynamics of the inflaton field both during and immediately after the inflationary phase. In particular, in Refs.  it was shown the importance of considering dissipative effects (due to inflaton’s decay) during inflation. The authors in  have shown how dissipation can influence the usual scenario of inflation, as also the effect of dissipation on the field trajectories in phase space. Strong dissipative regimes for field evolution have in special motivated the implementation of new inflationary models called “warm inflation” . The viability of the construction of microscopic motivated models of strong dissipative inflaton models has been shown possible in the context of superstring inspired models in Ref.  and argued also possible for more general models in .

On very general grounds, we are then lead to inquire about the possible role of field dissipation also at the early stages of inflation and in special to its role in smoothing large regions in phase space and thus easing considerable the homogeneity requirements for the onset of inflation and, consequently, providing a natural solution for the problem of fine-tuning associated to hybrid inflation. It is easy to understand the effect of dissipation on the initial conditions problem. For the regime of initial inflaton amplitude below the Planck scale the energy density is small at the beginning of inflation, resulting in a small Hubble parameter $`H`$, which determines how fast the inflaton energy is converted to expansion. For small values of $`H`$ the coupling of the fields to the background metric of gravity is small and the corresponding friction like term in the equations of motion, of the form $`3H\dot{\varphi }`$, is small. This results in a very long evolution for the fields, where those fields coupled to the inflaton oscillate around zero several times, generating strong sensitivity to very small variations of the initial conditions and resulting in a complete indetermination of the final outcome of the fields’ trajectories, which can be evolution toward an inflationary regime or to the minima of the potential. The behavior of the fields’ trajectories in the early times before a required long inflationary regime is then extremally chaotic. In the language of dynamical systems we have two types of attractors in the system, represented by the inflationary regime and the minima of the potential. The presence of dissipation works in damping the fluctuations of the fields in this initial critical time, suppressing the chaotic motions, which we then expect, from the results of Ref. , will bias the inflaton field trajectory toward the inflationary one. In fact, recently a study done by the authors in Ref.  have indicated that additional damping terms in the inflaton equation of motion could alleviate many of the problems related to the homogeneity requirements before inflation. Another study done by the authors in Ref. , using a supersymmetric hybrid inflation model, have also indicated that particle production may be a way to relax the extreme homogeneity requirements for hybrid inflation.

Here we will then be mostly concerned in determining for what magnitude of dissipation will the outcome of the evolution of the fields tend mostly toward the inflationary regime, or in the dynamical system approach, determine how fast dissipation will change the chaotic regime, characterized by unstable inflationary trajectories, to a non-chaotic one, with stable inflationary trajectories. We study chaos in our dynamical system of equations by means of the measure of the fractal dimension (or dimension information) \[for a review and definitions, see e.g., Ref. \], which gives a topological measure of chaos for different space-time settings and it is a quantity invariant under coordinate transformations, providing then an unambiguous signal for chaos in cosmology and general relativity problems in general . The method we apply in this work for quantifying chaos is then particularly useful in this cosmological pre-inflationary scenario context we are studying, in which case other methods may be ambiguous, like, for example, the determination of Lyapunov exponents, which does not give a coordinate invariant measure for chaos, as discussed in . Also, other methods for studying chaotic systems, like for example by Poincaré sections, are not suitable in the case we are interested here, in which case chaos is mostly a transitory phenomenon (it ends by the time the fields reach the potential minima, or when the inflaton enters the inflationary region). By computing the fractal dimension characterizing this chaotic behavior, which is related to the uncertainty in the system parameter values to predict the final outcome from a given initial condition (we may say that the boundaries of initial conditions that lead to inflation or evolution towards the potential minima are mixed) we are then able to infer the naturalness of inflation for different settings of the initial conditions. To our knowledge, this is the first time that this kind of study is performed in the context of a dissipative dynamical system in cosmology and applied to the initial condition problem of inflation in particular.

The paper is organized as follows: in Sec. II I give the basic equations defining the dynamical system and I discuss some of their properties. In Sec. III I present the numerical analysis of the dynamical system and the computation of the fractal dimension as a function of the magnitude of dissipation. From this analysis one will be able to conclude about the general effect of dissipation in the evolution of the fields and how it works in favor of the inflationary regime. Finally in Sec. IV I give the concluding remarks.

## II The Model and Its Properties

The model I will study here consists of the simplest hybrid inflation model, with a scalar field $`\varphi `$ (the inflaton) coupled to another scalar field $`\sigma `$, which triggers the end of inflation. The potential is

$$V(\varphi ,\sigma )=\frac{\lambda }{4}(\sigma ^2M^2)^2+\frac{m^2}{2}\varphi ^2+\frac{g^2}{2}\varphi ^2\sigma ^2,$$
(1)
where the parameters $`M^2,m^2`$ and the couplings $`\lambda `$ and $`g^2`$ are positive. In the following I will treat only the case of homogeneous fields, $`\varphi \varphi (t)`$ and $`\sigma \sigma (t)`$, which is the usual case in the works dealing with the initial condition problem in hybrid inflation. The interpretation of inflation from the potential (1) is the standard one. For values of $`\varphi `$ larger than a critical value $`\varphi _{\mathrm{cr}}`$, where $`\varphi _{\mathrm{cr}}^2=\lambda M^2/g^2`$, there is no symmetry breaking in the $`\sigma `$-field direction and $`\sigma =0`$ is a local minimum of the potential. We are interested in the region where $`M_{\mathrm{pl}}^2\varphi ^2\varphi _{\mathrm{cr}}^2`$ where inflation takes place (the false vacuum dominated regime). After the inflaton field drops below $`\varphi _{\mathrm{cr}}`$, symmetry breaking occurs in the $`\sigma `$ field direction of the potential. At this point the fields quickly move towards the minima $`\sigma =\pm M`$, $`\varphi =0`$ and inflation ends.

The chaotic properties of the classical (homogeneous) equations of motion, in Minkowski spacetime, for a model with potential similar to the one given by Eq. (1) were studied in Ref. , while the full effective equations of motion in Minkowski spacetime were studied by Ramos and Navarro in Ref. . In Ref.  the general form expected for the equations of motion for the fields, for potentials of the form of Eq. (1), was derived and a detailed account for the dissipative terms appearing in those equations was given. The chaotic behavior of the equations of motion for the fields was quantified by means of the fractal dimension and the authors studied in details how dissipation of the fields, due to decaying modes, changed the chaotic behavior of the dynamics. I will here then extend the results of Ref.  to the case of an expanding background, a flat Friedman-Robertson-Walker background metric.

As shown in Refs.  dissipation comes from the coupling of the system fields $`\varphi `$ and $`\sigma `$ (here taken as background field configurations) to a bath made of a set of other fields, for example made of $`N_\psi `$ fermions $`\psi _i`$ and/or $`N_\chi `$ scalars $`\chi _i`$. The general form of the interactions can be of the form

$$\underset{i=1}{\overset{N_\chi }{}}(g_\varphi ^2\varphi ^2\chi _i^2+g_\sigma ^2\sigma ^2\chi _i^2),$$
(2)
for scalar fields $`\chi _i`$ and

$$\underset{i=1}{\overset{N_\psi }{}}(f_\varphi \overline{\psi }_i\varphi \psi _i+f_\sigma \overline{\psi }_i\sigma \psi _i),$$
(3)
for fermion fields $`\psi _i`$. The usual form expected for the coupled effective equations of motion for $`\varphi `$ and $`\sigma `$ are complicated non-local equations of motion with typical non-Markovian dissipative kernels. In  we have shown that at high temperature and in the large $`N_\chi `$, $`N_\psi `$ limit for the (thermal) bath fields we can find an approximate Markovian limit for the kernels, from which we can express the equations of motion for the fields in terms of dissipative local equations. In Ref.  we have studied in details the general equations for the non-Markovian dissipative kernels at zero temperature and shown that for a certain class of field decaying modes the Markovian approximation is also a valid assumption. The results obtained in Ref.  show that we can evade the assumption of an initial high temperature thermal bath as used in Ref.  to derive the local equations of motion for the fields and makes it, therefore, more appropriate to the application we intend for here, which is the study of the initial conditions problem for the onset of inflation in hybrid inflation models characterized by potentials like Eq. (1).

I will just use the general form of the coupled effective equations of motion obtained in Ref.  without a complete specification of the dissipative coefficients, which we will take as free parameters. This is a valid analysis here since I will mostly be interested in how the general dissipation expected for the fields $`\varphi `$ and $`\sigma `$ in their coupled effective equations of motion will change the chaotic properties of the dynamical system, which, as explained in the introduction, are direct related to the initial conditions problem. The detailed form of the dissipative coefficients depends on the microscopic physics of the specific model under study, like the coupling of the system fields with the bath fields, the decaying modes available and the expanding metric. In fact the magnitude of the dissipation terms may be controlled by these couplings, as shown in Refs. .

In analogy to the results obtained in Ref.  we write the equations of motion for $`\varphi `$ and $`\sigma `$ in the form***This should be compared with Eqs. (5.1) and (5.2) of Ref. . In the notation used here we have that in those equations that $`\varphi _c\varphi `$, $`\psi _c\sigma `$, $`\overline{m}_\varphi ^2m^2`$, $`\overline{m}_\psi ^2\lambda M^2`$, $`\overline{\lambda }_\varphi 0`$, $`\overline{\lambda }_\psi 6\lambda `$ and $`\overline{g}g`$.

$`\ddot{\varphi }+3H\dot{\varphi }+m^2\varphi +g^2\varphi \sigma ^2+\eta _1\varphi ^2\dot{\varphi }+\eta _3\varphi \sigma \dot{\sigma }=0`$ (4)
and

$`\ddot{\sigma }+3H\dot{\sigma }\lambda M^2\sigma +\lambda \sigma ^3+g^2\sigma \varphi ^2+\eta _2\sigma ^2\dot{\sigma }+\eta _3\varphi \sigma \dot{\varphi }=0,`$ (5)
where $`\eta _1`$, $`\eta _2`$ and $`\eta _3`$ denote the dissipation coefficients. Besides the coupled equations of motion (4) and (5) we also have the Friedman equation and the evolution equation for the Hubble parameter, $`H=\dot{a}/a`$, where $`a`$ is the scale factor, are given as usual by

$$H^2=\frac{8\pi G}{3}(\rho _m+\rho _r)\frac{k}{a^2},$$
(6)
$$2\dot{H}+3H^2+\frac{k}{a^2}=8\pi G(p_m+p_r),$$
(7)
where $`G=1/m_{\mathrm{Pl}}^2`$, with $`m_{\mathrm{Pl}}`$ the Planck mass. $`k=0,+1,1`$ for a flat, closed or open Universe, respectively. $`\rho _{m(r)}`$ and $`p_{m(r)}`$ are the energy density and pressure for matter (radiation), respectively. We also have the standard relations:

$$\rho _m=\frac{1}{2}\dot{\varphi }^2+\frac{1}{2}\dot{\sigma }^2+V(\varphi ,\sigma ),$$
(8)
$$p_m=\frac{1}{2}\dot{\varphi }^2+\frac{1}{2}\dot{\sigma }^2V(\varphi ,\sigma )$$
(9)
and $`p_r=\frac{1}{3}\rho _r`$.

The matter and radiation energy densities $`\rho _m`$ and $`\rho _r`$ evolve in time as:

$$\dot{\rho }_m+3H\left(\dot{\varphi }^2+\dot{\sigma }^2\right)+\eta _1\varphi ^2\dot{\varphi }^2+\eta _2\sigma ^2\dot{\sigma }^2+2\eta _3\varphi \sigma \dot{\varphi }\dot{\sigma }=0$$
(10)
and (from the energy conservation law)

$$\dot{\rho }_r+4H\rho _r\eta _1\varphi ^2\dot{\varphi }^2\eta _2\sigma ^2\dot{\sigma }^22\eta _3\varphi \sigma \dot{\varphi }\dot{\sigma }=0.$$
(11)
Assuming a flat universe ($`k=0`$), from Eq. (6), we can consider

$$\rho _r=\frac{3}{8\pi G}H^2\rho _m=\frac{3}{8\pi G}H^2\frac{\dot{\varphi }^2}{2}\frac{\dot{\sigma }^2}{2}V(\varphi ,\sigma )$$
(12)
as the first integral of Eq. (11). Using Eqs. (6)-(9), we can also express the equation for the acceleration in the following form

$$\frac{\ddot{a}}{a}=\frac{8\pi G}{3}(\rho _m\rho _r)4\pi G\left(\dot{\varphi }^2+\dot{\sigma }^2\right).$$
(13)
The field equations (4) and (5) with (13) form a dissipative dynamical system that we will study numerically. Using (12) in (13) and defining the following dimensionless variables: $`y=\sqrt{8\pi G}\varphi `$, $`x=\sigma /M`$, $`\tau =Mt`$ and the constants $`\alpha ^2=m^2/M^2`$, $`\beta ^2=8\pi GM^2`$ and the rescaled dissipative coefficients $`\eta _y=M\eta _1`$, $`\eta _x=M\eta _2`$ and $`\eta _{xy}=M\eta _3`$, the system of equations (4), (5) and (13) can then be rewritten in terms of dimensionless variables in the form of the following system of first order differential equations:

$`x^{}=z,`$ (14)
$`y^{}=w,`$ (15)
$`a^{}=u,`$ (16)
$`z^{}=3{\displaystyle \frac{u}{a}}z+\lambda x\lambda x^3{\displaystyle \frac{g^2}{\beta ^2}}xy^2\eta _xx^2z{\displaystyle \frac{1}{\beta ^2}}\eta _{xy}xyw,`$ (17)
$`w^{}=3{\displaystyle \frac{u}{a}}w\alpha ^2yg^2yx^2{\displaystyle \frac{1}{\beta ^2}}\eta _yy^2w\eta _{xy}xyz,`$ (18)
$`u^{}={\displaystyle \frac{u^2}{a}}{\displaystyle \frac{a}{6}}\left(\beta ^2z^2+w^2\beta ^2\lambda +2\lambda \beta ^2x^2\lambda \beta ^2x^42\alpha ^2y^22g^2x^2y^2\right),`$ (19)
where prime indicates derivative with respect to the dimensionless time $`\tau `$, e.g., $`x^{}=dx/d\tau `$.

For convenience we use the values for $`\lambda `$, $`g^2`$, $`M`$ and $`m`$ as given by Mendes and Liddle in Ref. : $`\lambda =g^2=1`$, $`M=2\times 10^2M_{\mathrm{Pl}}`$ and $`m=5\times 10^6M_{\mathrm{Pl}}`$, where $`M_{\mathrm{Pl}}=1/\sqrt{8\pi G}`$ is the reduced Planck mass. From these values of $`m`$ and $`M`$ we find for the constants $`\alpha `$ and $`\beta `$ in Eq. (19) the values: $`\alpha =2.5\times 10^4`$ and $`\beta =2\times 10^2`$. In terms of these values we also have that the critical value for the inflaton field, $`\varphi _{\mathrm{cr}}`$, in the dimensionless variables, is given by $`y_{\mathrm{cr}}=2\times 10^2`$.

## III Chaos and the Fractal Dimension

We next numerically solve the system of equations in Eq. (19) and we search for chaotic regimes and how they change with increasing dissipation. Lets initially consider the following initial conditions at $`t=0`$: $`\varphi =4\times \varphi _{\mathrm{cr}}`$, $`\dot{\varphi }=0`$, $`\sigma =2\times 10^2M`$, $`\dot{\sigma }=0`$ and $`H`$ is determined initially by Eq. (6), with zero initial radiation energy density. These values guarantee that the stable inflationary trajectories will correspond to at least $`N_e60`$ e-folds of inflation. For the dissipation coefficients $`\eta _x`$, $`\eta _y`$ and $`\eta _{xy}`$, for convenience, we will consider them all with the same value, which is consistent with the recent findings in Ref. , or the calculations in Ref. . Small variations of the magnitude among these coefficients are not critical here, since the dissipation in Eq. (19) is dominated by the term corresponding to $`\eta _1\varphi ^2\dot{\varphi }`$. We take the dissipation coefficients with initial values as given by $`\{\eta \}=(\eta _x,\eta _y,\eta _{xy})=(0,0,0)`$ and increase them till the chaotic behavior of the system changes to non-chaotic.

For $`\{\eta \}=0`$ a typical result for three trajectories around the initial values of the fields (in dimensionless units), separated by a $`\pm 10^5`$ variation, is shown in Fig. 1.

In Fig. 1 we can distinguish the inflationary trajectory represented by the straight line at $`x=0`$, which represents the evolution of the fields in phase space along the valley of local minima of the potential, $`\sigma =0`$ for $`\varphi >\varphi _{\mathrm{cr}}`$. This trajectory is characterized by a very long evolution by which the universe expands exponentially. The number of e-folds produced by the end of the inflationary phase, when the inflaton reaches $`\varphi \varphi _{\mathrm{cr}}`$, for this particular trajectory is $`N_e350`$. The other two oscillatory trajectories, that evolve towards $`x=\pm 1`$ represent the non-inflationary trajectories, where the fields quickly evolve towards the potential minima $`\sigma =\pm M`$ and $`\varphi =0`$. In Fig. 2 we show a blow-out of the initial time evolution of Fig. 1. It clearly shows highly oscillatory chaotic like behavior.

The chaotic behavior of the dynamical system is quantified by means of the determination of the fractal dimension of the boundary separating the inflationary trajectories from those that are not inflationary, i.e., the ones that evolve towards the minima of the potential, $`\varphi =0`$ and $`\sigma =\pm M`$. The fractal dimension is associated with the possible different exit modes under small changes of the initial conditions at $`t=0`$ and it will give a measure of the degree of chaos of our dynamical system. The exit modes we refer to above are ones of the symmetry breaking minima in the $`\sigma `$-field direction, $`\sigma =\pm M`$, $`\varphi =0`$ and the inflationary one, which are attractors of field trajectories in phase space. The method we employ to determine the fractal dimension is the box-counting method which is a standard method for determining the fractal dimension of boundaries . Its definition and the specific numerical implementation we use here have been described in details in Ref. .

The basic procedure is, given a set of initial conditions $`𝐱_0`$ at $`t=t_0`$, which leads to a certain outcome for the trajectory in phase space, by perturbing them by an amount $`\delta `$ we then study whether there will be a change of outcome for the trajectory or not (whether the perturbation will lead to a different attractor or not). Given a volume region in phase space around a boundary between different attractors and perturbing a large set of initial conditions inside that region, the fraction of uncertain trajectories, $`f(\delta )`$, which result in a different outcome under a small perturbation, can be shown to scale with the perturbation $`\delta `$ as  $`f(\delta )\delta ^ϵ`$, where $`ϵ`$ is called the uncertainty exponent. The box-counting dimension of the boundary in phase space separating different attractors, or fractal dimension $`f_d`$, is given by  $`f_d=dϵ`$, where $`d`$ denotes the dimension of the phase space. For a fractal boundary $`f_d>d1`$, implying that $`ϵ<1`$, whereas for a non-fractal boundary, $`f_d=d1`$, and $`ϵ=1`$.

Following the method of box-counting we then consider a box in phase space (for the dimensionless variables) of size $`10^5`$, around the initial conditions used in Figs. 1 and 2, inside which a large number of random points are taken (a total of $`100.000`$ random points were used in each run). All initial conditions are then numerically evolved by using an eighth-order Runge-Kutta integration method and the fractal dimension is obtained by statistically studying the outcome of each initial condition for each run of the large set of points. Special care is taken to keep the statistical error in the results always below $`1\%`$. From these numerical simulations we obtain for the zero dissipative regime of Fig. 1 the result for the fractal dimension as given by $`f_d5.80\pm 0.05`$. This corresponds to the dimension of the dynamical system, Eq. (19), corresponding here to $`d=6`$, minus the uncertainty coefficient $`ϵ`$ (which gives a measure of how chaotic is the system ), $`f_d=6ϵ`$. The results obtained by increasing the value of the dissipation coefficients $`\{\eta \}`$, while keeping the same initial condition used above, around which the perturbations are taken, are shown in Fig. 3.

The results in Fig. 3 show that the system quickly changes to a non-chaotic behavior for a dissipation coefficient around $`\{\eta \}3.1\times 10^4`$. For those values of dissipation and higher all trajectories are inflationary ones with $`N_e60`$. At that point we then have the breakdown of the fractal structure of the boundary between inflationary and non-inflationary trajectories. The boundary becomes smooth and we can reliable predict that the initial conditions will evolve towards the inflationary region.

In order to study how efficient is dissipation to drive non-inflationary trajectories towards the inflationary region, which is characterized by $`\ddot{a}>0`$, we show in Fig. 4 the results for the number of e-folds of inflation as a function of dissipation for an initial trajectory that in the absence of dissipation would be a non-inflationary one. The particular initial conditions we take purposefully correspond to fields well below the Planck scale, $`\varphi =3.5M`$, $`\sigma =M`$, and non-vanishing field derivatives $`\dot{\varphi }=2.5\times 10^2M^2`$ and $`\dot{\sigma }=5\times 10^4M^2`$, which in the absence of dissipation would be quickly driven to one of the minima of the potential (and therefore corresponding to a non-inflationary behavior). We see that above some value of dissipation coefficients $`\eta 1`$ the number of e-folds of inflation increases fast for rather very small increase of dissipation. We then clearly see using the set of initial conditions above how dissipation acts in turning non-inflationary trajetories into inflationary ones, shifting the system deep into the inflationary region as dissipation is increased. The initial radiation energy density in this case, as in the previous one, is taken initially zero and remains always smaller than the vacuum energy density during the whole evolution of the inflaton. In particular, for the values of dissipation coefficients we have studied here, after a few e-folds of inflation the radiation energy density is a negligible fraction of the vacuum energy density and right before the end of inflation the radiation energy density is a tiny fraction of the vacuum energy density, $`\rho _r/\rho _v10^{14}`$. A detailed treatment of the liberated radiation to the medium and its eventual translation to a temperature depends of course on how the radiation thermalization rate $`\mathrm{\Gamma }_{\mathrm{rad}}`$, which is determined by the microscopic physics, compares to the expansion rate $`H`$.

It is also useful to compare the regimes of field dissipation corresponding for instance to the ones studied in Fig. 3 and Fig. 4 and their relation to the warm inflation scenario. In warm inflation  we are usually interested in the regime of strong dissipation in the effective equation of motion for the inflaton, as compared with the expansion term $`3H\dot{\varphi }`$. From Eq. (4) dissipation is dominated by the term $`\eta _1\varphi ^2\dot{\varphi }=\eta (\varphi )\dot{\varphi }`$. From the values of parameters and initial conditions used in Fig. 3 we find that $`3H/\eta (\varphi )5`$ for the largest value of dissipation studied, while for the parameters and initial conditions used in Fig. 4 we have that $`3H/\eta (\varphi )0.008`$ and therefore deep in the regime of warm inflation. We can easily see then that warm inflation kind of models does not seem to suffer from the initial conditions problem seen in standard hybrid inflation models. In particular, for the parameters used in Fig. 3, for values of dissipation $`\eta 0.0015`$ we enter in the warm inflation regime and the system has well before became non-chaotic, with stable inflationary trajectories.

## IV Conclusions

We have examined a dynamical system describing hybrid inflation in the presence of dissipation for both the inflaton and the field that triggers the end of inflation. Dissipation is taken as being present throughout the system’s evolution. In the absence of dissipation the dynamical system is shown to be highly chaotic by means of the measure of the fractal dimension. The fractal dimension gives a coordinate invariant measure for chaos and therefore is an appropriate way to quantify chaos in cosmology.

For a given initial condition that evolves to an inflationary region, by varying it by a small perturbation we completely loose track of its final outcome, which may still evolve towards the inflationary region, characterized by evolution along the valley of local minima of the potential, characterizing the inflationary regions, or fast evolution towards one of the symmetry breaking minima. The time evolution of the fields are shown to be highly oscillatory for a very long time, which eventually change the inflationary trajectories to non-inflationary ones. This is the main characteristics in the evolution of the fields and that gives rise to the initial condition problem in hybrid inflation (and in a lesser degree in any inflationary model). The long oscillatory behavior of the fields is highly chaotic. By measuring the magnitude of chaos in the system we can then infer the severity of the initial conditions problem.

However, by changing the magnitude of dissipation we have then shown that the highly chaotic evolution of the system, with unstable inflationary trajectories, can change quickly to a non-chaotic, stable one, where all trajectories eventually evolve toward the inflationary region. Dissipation acts by damping the initial oscillatory behavior of the fields and then making it possible to provide a solution for the initial condition problem. We have shown that this solution can be obtained with relative very small dissipation of the initial inflaton amplitude, with rates of conversion of vacuum energy to radiation as small as $`10^{14}`$. Here dissipation works in two ways. By damping the fluctuations associated to the fields coupled to the inflaton the fine-tuning problem associated with hybrid inflation models is avoided and at the same time it makes possible for initial inflationary trajectories to enter faster in the inflationary region (characterized by $`\ddot{a}>0`$) and stay longer in there.

The method used here to study the initial condition problem can also be easily extented to include other variants of the hybrid inflation model, for example in the context of supersymmetric motivated potentials, where more than one field is coupled to the inflaton. The additional equations of motions will also generically be dissipative equations and the same changing of behavior due to dissipation observed here is also expected to appear in those more complicated models.

There have been previous studies on chaos in hybrid inflation models , but these studies have concentrate on the reheating period after the inflationary phase. In that context, the effect of dissipation, that would be inherently present in the field equations also during this period, has been neglected. It would then be interesting to investigate in those cases also the effect of dissipation and its role, for example, in the phenomenon of parametric resonance during preheating in hybrid inflation models . We are currently investigating this and other applications and we expect to report on them soon.

###### Acknowledgements.

I would like to thank R. Caldwell and M. Gleiser for many discussions in the subject of this paper and on the applications of dissipative regimes for scalar fields in cosmology. The author was supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq (Brazil) and SR-2 (UERJ).


format: text

text:
# 1 Introduction

## 1 Introduction

It has long been known that in certain restricted regions of phase space the coefficients of the QCD perturbation expansion can be parametrically large. For semi-inclusive observables this is often due to the incomplete cancellation of soft and/or collinear logarithms.

A set of variables particularly sensitive to this kind of problem is event-shape distributions $`\frac{1}{\sigma }\frac{d\sigma }{dv}`$ which we shall use as an illustration throughout. Typically one finds that the perturbative series is dominated by terms $`(\alpha _\text{s}^n\mathrm{ln}^{2n1}v)/v`$, where $`v1`$ denotes the value of the event-shape variable. In order to have a meaningful perturbative description these terms have to be resummed to all orders. This has been done for a range of observables .

Most of the observables considered so far have the property of exponentiation, implying the following structure for the resummation,

$$\mathrm{\Sigma }(v)=\left(1+\underset{n}{}C_n\overline{\alpha }_\text{s}^n\right)e^{Lg_1(\alpha _\text{s}L)+g_2(\alpha _\text{s}L)+\mathrm{}}+D(v),$$
(1)
where the distribution is given by $`d\mathrm{\Sigma }/dv`$. The remainder function $`D(v)`$ has the property that it goes to zero for $`v0`$; $`L=\mathrm{ln}1/v`$ and $`\overline{\alpha }_\text{s}=\alpha _\text{s}/2\pi `$. The terms resummed by $`Lg_1(\alpha _\text{s}L)`$ are referred to as leading (or double) logarithms (LL), while those in $`g_2(\alpha _\text{s}L)`$ are next-to-leading (or single) logarithms (NLL). For certain variables, rather than a single exponent one finds a sum or integral over exponents .

A property common to most of the observables considered until recently is that they are sensitive to emissions anywhere in phase space. We call these *global* event shapes. Examples are the thrust, the heavy-jet mass, the $`C`$-parameter. In this paper we are interested in *non-global* observables, those sensitive to radiation in only a part of phase space. We shall discuss mainly single-hemisphere event-shapes in $`e^+e^{}`$ (especially the jet-mass). Other examples include a number of DIS current-hemisphere observables and quantities such as the mass and broadening of single jets in multi-jet ensembles.

The literature on the subject of global event shapes often makes use of a combination of results in separate hemispheres to obtain the properties of the whole event . When the separate hemisphere results are combined in this fashion the answer obtained is correct. However this can be misleading, since it turns out that the existing single-hemisphere results, applied literally to just a single hemisphere event-shape, are incomplete at NLL level. This has not been taken into account in the recently presented results for the light-jet mass and the narrow-jet broadening  and for the narrow-jet thrust-minor . In the case of the first two variables the incompleteness of the results can be seen from a comparison of the $`𝒪\left(\alpha _\text{s}^2L^2\right)`$ term of the expansion of the resummation with fixed-order predictions from a Monte Carlo program such as Event2 .222Here we refer to the original version of  where it was stated that terms down to $`\alpha _\text{s}^nL^{2n2}`$ of the expanded result were under control. In a revised version (v2) which appeared on the hep-ph arXiv shortly before the submission of our paper, the claims regarding the accuracy were revised to state that only terms down to $`\alpha _\text{s}^nL^{2n1}`$ are controlled. In our notation that corresponds to controlling only the $`\alpha _\text{s}L`$ term of $`g_2`$, but not terms $`\alpha _\text{s}^nL^n`$ with $`n2`$, nor $`C_1`$. Here we follow the convention adopted in  rather than that in , and use the term next-to-leading to refer to logarithms in the exponent(s), not in the expansion of the exponent.

To illustrate the point it is convenient to study a simpler variable, the hemisphere jet-mass, defined as

$$\rho =\frac{\left(_{i\mathrm{hemisphere}}k_i\right)^2}{\left(_iE_i\right)^2},$$
(2)
where the hemisphere is taken with respect to the thrust axis. Once the distribution of $`\rho `$ has been determined, it can be combined with the known heavy-jet mass distribution  to obtain the light-jet mass distribution, using the following exact relation

$$\frac{d\sigma }{d\rho _{\mathrm{light}}}+\frac{d\sigma }{d\rho _{\mathrm{heavy}}}=\frac{d\sigma }{d\rho _{\mathrm{left}}}+\frac{d\sigma }{d\rho _{\mathrm{right}}}2\frac{d\sigma }{d\rho },$$
(3)
where $`\rho `$ is the hemisphere jet-mass. An analogous relation applies for the jet-broadenings. For the narrow-jet thrust minor the situation is more complex but considerations similar to those addressed here will apply.

One can understand the origin of the qualitative difference between the heavy-jet and single-hemisphere mass resummations by examining figure 1, which depicts two soft large-angle gluons, whose energies are ordered $`Qk_1k_2`$. For the heavy-jet mass there is a cancelation between the two contributions because the emission of $`k_2`$ does not affect the value of $`\rho _{\mathrm{heavy}}`$ (which is determined by the harder emission, $`k_1`$). In the case of the right-jet mass, $`k_1`$ has no effect on $`\rho _{\mathrm{right}}`$, so the event shape receives a contribution only from the softer emission, $`k_2`$. This spoils the cancelation with the virtual diagram, figure 1a, thus leading to a term $`\alpha _\text{s}^2L^2`$, i.e. a NLL effect. Such effects are a general feature of non-global observables.

At all orders, to NLL accuracy, the form of the resummed jet-mass distribution will be

$$\mathrm{\Sigma }(\rho )=\left(1+\overline{\alpha }_\text{s}C_1^{(q)}\right)𝒮(\alpha _\text{s}L)\mathrm{\Sigma }_q(\alpha _\text{s},L)+\overline{\alpha }_\text{s}C_1^{(g)}\mathrm{\Sigma }_g(\alpha _\text{s},L)+D(\rho ).$$
(4)
$`\mathrm{\Sigma }_q`$ and $`\mathrm{\Sigma }_g`$ are the resummed jet-mass cross sections which can be obtained by integrating the quantities $`J_q`$ and $`J_g`$ first derived in  and widely used in the literature. They have structures like that of the exponential factor in eq. (1). The sum of two resummed contributions arises because in a fraction $`𝒪\left(\alpha _\text{s}\right)`$ of events the hard particle in the right hemisphere ($`_\mathrm{R}`$) is a gluon.

The new effect considered in this article is embodied by the function $`𝒮`$.333In principle a similar function multiplies $`\mathrm{\Sigma }_g`$, but at NLL accuracy it can be neglected. We write its expansion as

$$𝒮(\alpha _\text{s}L)=1+\underset{n=2}{}S_n(\overline{\alpha }_\text{s}L)^n.$$
(5)
It is straightforward to exactly compute the first non-trivial term $`S_2`$ and this is done in the following section. The full computation of $`𝒮`$ involves considering an ensemble of an arbitrary number of large-angle energy-ordered soft gluons in $`_\mathrm{L}`$, which coherently emit a single, softer gluon into $`_\mathrm{R}`$. For reasons elucidated later it is difficult to carry out an all-orders treatment of such an effect analytically. We therefore opt to treat these effects using a Monte Carlo algorithm valid in the large-$`N_C`$ limit. This is outlined in section 3 and further details are given in the appendix.

Finally in section 4 we compare our results to the $`𝒪\left(\alpha _\text{s}^2\right)`$ predictions of Event2. Phenomenological predictions including this effect will be shown elsewhere .

## 2 Fixed order calculation

First we calculate the contribution to the jet-mass distribution from the configuration in figure 1b, considering the right-hemisphere jet for concreteness. We introduce the following particle four-momenta

$`k_a`$ $`={\displaystyle \frac{Q}{2}}(1,0,0,1),`$ (6a)
$`k_b`$ $`={\displaystyle \frac{Q}{2}}(1,0,0,1),`$ (6b)
$`k_1`$ $`=x_1{\displaystyle \frac{Q}{2}}(1,0,\mathrm{sin}\theta _1,\mathrm{cos}\theta _1),`$ (6c)
$`k_2`$ $`=x_2{\displaystyle \frac{Q}{2}}(1,\mathrm{sin}\theta _2\mathrm{sin}\varphi ,\mathrm{sin}\theta _2\mathrm{cos}\varphi ,\mathrm{cos}\theta _2),`$ (6d)
where we have labelled the quark and antiquark as $`a`$ and $`b`$ and defined energy fractions $`x_{1,2}1`$ for the two gluons. We have ignored recoil in the kinematics, because the jet-mass is insensitive to it.

When gluon $`2`$ is in $`_\mathrm{R}`$ the jet mass has the value $`\rho =x_2(1\mathrm{cos}\theta _2)/2`$. When only the quark is in $`_\mathrm{R}`$, $`\rho =0`$.

We write the matrix element for ordered two-gluon emission as (see for example )

$`W`$ $`=4C_F{\displaystyle \frac{(ab)}{(a1)(1b)}}\left({\displaystyle \frac{C_A}{2}}{\displaystyle \frac{(a1)}{(a2)(21)}}+{\displaystyle \frac{C_A}{2}}{\displaystyle \frac{(b1)}{(b2)(21)}}+\left(C_F{\displaystyle \frac{C_A}{2}}\right){\displaystyle \frac{(ab)}{(a2)(2b)}}\right)`$ (7a)
$`=C_F^2W_1+C_FC_AW_2,`$ (7b)
where $`(ij)=k_ik_j`$. The result is valid for $`1x_1x_2`$ as well as for the opposite ordering of the gluons, and furthermore is completely symmetric under interchange of $`k_1`$ and $`k_2`$. We have however chosen to write it in an asymmetric form to emphasise the dipole structure of the emissions, namely radiation of gluon $`k_1`$ from the $`ab`$ dipole, followed by the radiation of gluon $`k_2`$ from the $`a1`$, $`1b`$ and $`ab`$ dipoles.

The $`C_F^2`$ piece of the matrix element, $`W_1`$ corresponds to independent gluon emission and is included in the usual resummation of the quark jet-mass, *i.e.* in $`\mathrm{\Sigma }_q`$ of eq. (4). To study specifically the modification relative to the standard quark-jet mass result due to configurations like that shown in figure 1b, one must consider the $`C_FC_A`$ part, $`W_2`$, of the emission probability above.

The integral to be considered, related to the probability for the jet mass to be less than $`\rho `$, is then

$$\begin{array}{c}C_FC_A\left(\frac{\alpha _s}{2\pi }\right)^2_1^0d\mathrm{cos}\theta _1_0^1d\mathrm{cos}\theta _2_0^{2\pi }\frac{d\varphi _1}{2\pi }_0^{2\pi }\frac{d\varphi _2}{2\pi }\hfill \\ \hfill \frac{Q^4}{16}_0^1x_2𝑑x_2_{x_2}^1x_1𝑑x_1\mathrm{\Theta }\left(\frac{x_2}{2}(1\mathrm{cos}\theta _2)\rho \right)W_2,\end{array}$$
(8)
whose dominant term at small $`\rho `$ is $`S_2\overline{\alpha }_\text{s}^2\mathrm{ln}^21/\rho `$.

The integrals over the energy fractions are straightforward. Keeping only the $`\mathrm{ln}^21/\rho `$ piece and performing the azimuthal average we have

$$S_2=2C_FC_A_0^1d\mathrm{cos}\theta _2_1^0d\mathrm{cos}\theta _1\mathrm{\Omega }_2(\mathrm{cos}\theta _1,\mathrm{cos}\theta _2),$$
(9)
where the angular function $`\mathrm{\Omega }_2`$ is

$$\mathrm{\Omega }_2=\frac{2}{(\mathrm{cos}\theta _2\mathrm{cos}\theta _1)(1\mathrm{cos}\theta _1)(1+\mathrm{cos}\theta _2)}.$$
(10)
Integrating over the polar angles we obtain

$$S_2=C_FC_A\frac{\pi ^2}{3}.$$
(11)
We note that possibly related “$`\pi ^2`$” terms have been discussed in .

## 3 Resummation

To perform the all orders resummation of $`𝒮(\alpha _\text{s}L)`$ we need to consider general configurations of soft gluons. We define the probability $`P_𝒞(L)`$ for a given configuration $`𝒞`$ at a resolution scale $`L`$ (*i.e.* not resolving gluons with energies less than $`Qe^L`$). Given $`P_𝒞`$, the probability $`P_𝒞^{}`$ of a configuration $`𝒞^{}`$ with one extra gluon with scale $`L^{}>L`$, polar angle $`\theta ^{}`$ and azimuth $`\varphi ^{}`$ is (see for example )

$$P_𝒞^{}(L^{})=\overline{\alpha }_\text{s}\left(L^{}\right)\mathrm{\Delta }_𝒞(L,L^{})F_𝒞(\theta ^{},\varphi ^{})P_𝒞(L),$$
(12)
where the form factor is given by

$$\mathrm{ln}\mathrm{\Delta }_𝒞(L,L^{})=_L^L^{}𝑑L^{\prime \prime }d\mathrm{cos}\theta d\varphi \overline{\alpha }_\text{s}\left(L^{\prime \prime }\right)F_𝒞(\theta ,\varphi ),$$
(13)
and $`F_𝒞(\theta ,\varphi )`$ describes the angular and colour-structure of the radiation pattern from configuration $`𝒞`$. The angular integrations in (13) diverge whenever the angle is collinear to one of the emitting particles. In practice it is therefore convenient to introduce an angular cutoff $`ϵ`$ to regulate these divergences both in the real emissions and the virtual corrections.

We obtain $`𝒮`$ by calculating the probability of there being no emissions in $`_\mathrm{R}`$ down to a scale $`L`$, divided by the corresponding probability had the only source of emissions been the original $`q\overline{q}`$ pair. This gives

$$𝒮(\alpha _\text{s}L)=\frac{1}{\sqrt{\mathrm{\Delta }_{ab}(L)}}\underset{𝒞|_\mathrm{R}\mathrm{empty}}{}P_𝒞(L),$$
(14)
where the sum runs over all configurations $`𝒞`$ which contain no emissions in $`_\mathrm{R}`$.

### 3.1 Large-$`𝑵_𝑪`$ limit and Monte Carlo implementation

In practice two considerations make it difficult to implement the above approach analytically. One is that the colour algebra involved in the determination of $`F_𝒞`$ becomes progressively more complicated as the number of gluons increases. The other problem is simply that the treatment of the geometry of the many large-angle gluon ensemble quickly becomes prohibitive. The first problem can be partially solved by taking the large-$`N_C`$ limit. To address the issue of the geometry we shall use a Monte Carlo approach.

In the large-$`N_C`$ limit, one can represent gluons by pairs of colour-anticolour lines, as illustrated in figure 2. When squaring the amplitude one ignores contributions that in terms of their colour flow are topologically non-planar, because they are suppressed by powers of $`1/N_C^2`$ . As a result, for an ensemble of $`n`$ gluons, $`F_𝒞`$ just reduces to a sum of independent emission intensities from $`n+1`$ independent dipoles:

$$F_𝒞(\theta _k,\varphi _k)=\underset{\mathrm{dipoles}ij}{}\frac{2C_A}{(1\mathrm{cos}\theta _{ik})(1\mathrm{cos}\theta _{kj})}.$$
(15)
When the dipole $`ij`$ radiates a gluon $`k`$ it splits into two dipoles, $`ik`$ and $`kj`$. Thus the dipole structure is determined by the history of the gluon branching.

Such a branching pattern can be very naturally implemented using a Monte Carlo algorithm. At first sight one might envisage calculating the two factors in (14) separately, using the $`F_𝒞`$ to generate the distribution of radiation for each new configuration. However because of the collinear divergence along the direction of the quark in $`_\mathrm{R}`$, only a tiny fraction of events would be free of emissions in $`_\mathrm{R}`$ and so contribute to the sum in (14). The sum would therefore have a large relative error, which would translate to a large absolute error on $`𝒮`$ because of the division by the small quantity $`\sqrt{\mathrm{\Delta }_{ab}(L)}`$.

Instead a more efficient procedure involves moving the division by $`\sqrt{\mathrm{\Delta }_{ab}(L)}`$ directly into the calculation of the $`P_𝒞`$. This can be achieved using a modified radiation intensity, $`\stackrel{~}{F}_𝒞`$ (for both the emissions and the virtual corrections),

$$\stackrel{~}{F}_𝒞(\theta ,\varphi )=F_𝒞(\theta ,\varphi )F_{ab}(\theta ,\varphi )\mathrm{\Theta }(\theta ),$$
(16)
where one subtracts out the radiation intensity $`F_{ab}`$ which would have been produced by the original $`q\overline{q}`$ pair (in the large-$`N_C`$ limit). One calculates quantities $`\stackrel{~}{P}_𝒞`$ using analogs of eqs. (12) and (13) with $`F_𝒞`$ replaced by $`\stackrel{~}{F}_𝒞`$ and then $`𝒮`$ is simply given by

$$𝒮(\alpha _\text{s}L)=\underset{𝒞|_\mathrm{R}\mathrm{empty}}{}\stackrel{~}{P}_𝒞(L).$$
(17)
It should be kept in mind that since $`\stackrel{~}{F}_𝒞`$ is negative in certain regions of phase space one loses a strict probabilistic interpretation for the $`\stackrel{~}{P}_𝒞`$. Nevertheless the sum over configurations is well-defined and meaningful.

The exact details of the Monte Carlo algorithm are given in the appendix. Here we restrict ourselves to giving a parameterisation for $`𝒮`$ obtained by fitting to the Monte Carlo results:

$$𝒮(\alpha _\text{s}L)\mathrm{exp}\left(C_FC_A\frac{\pi ^2}{3}\left(\frac{1+(at)^2}{1+(bt)^c}\right)t^2\right),$$
(18)
with

$$t(\alpha _\text{s}L)=\frac{1}{2\pi }_{e^L}^1\frac{dx}{x}\alpha _\text{s}(xQ)=\frac{1}{4\pi \beta _0}\mathrm{ln}\frac{1}{12\beta _0\alpha _\text{s}L},$$
(19)
where $`\beta _0=(11C_A2n_f)/(12\pi )`$ and

$$a=0.85C_A,b=0.86C_A,c=1.33.$$
(20)
The parameterisation should be accurate to the order of a few percent (better in most of the region) for $`t<0.7`$, corresponding to $`12\alpha _\text{s}\beta _0L0.005`$.444The accessible range of $`t`$ is limited by two issues: firstly only a small fraction of events are generated at large $`t`$, requiring considerable statistics in order to investigate that region; and secondly because an accurate determination of $`𝒮`$ at large $`t`$ requires a very small angular cutoff, which leads to there being many dipoles in an event, and a consequent slowing down of the evolution.

Actually, for the purposes of the fit one replaces $`C_FC_A`$ in (18) with $`C_A^2/2`$ since the Monte Carlo works in the large-$`N_C`$ limit. But for use in phenomenology one wishes to have the exact colour structure at least at $`𝒪\left(\alpha _\text{s}^2\right)`$, hence the use of $`C_FC_A`$ in (18).

## 4 Checks and conclusions

It is useful to check our results against fixed order results from the next-to-leading order Monte Carlo program Event2 . First it is necessary to determine the constant terms $`C_1^{(q)}`$ and $`C_1^{(g)}`$, which are obtained by requiring consistency between (4) and a full $`𝒪\left(\alpha _\text{s}\right)`$ calculation. It is straightforward to show that they are given by

$$C_1^{(q)}=\frac{1}{2}\left(C_1^\tau r_3\right),C_1^{(g)}=\frac{r_3}{2},$$
(21)
where $`C_1^\tau `$ is the constant determined for the thrust in  and $`\overline{\alpha }_\text{s}r_3`$ is the probability of there being a hard gluon along the thrust axis, as given in :555Only in the hep-ph arXiv version 1.

$$r_3=C_F\left(2\mathrm{ln}^22\frac{5}{4}\mathrm{ln}3+4\mathrm{L}\mathrm{i}_2\left(\frac{1}{2}\right)+\frac{\pi ^2}{3}\frac{1}{6}\right)0.917C_F.$$
(22)
The plots in figure 3 show the difference between the $`𝒪\left(\alpha _\text{s}^2\right)`$ exact and resummed distributions. To NLL accuracy we aim to account for terms down to $`\alpha _\text{s}^2L^2`$ in $`\mathrm{\Sigma }`$ and therefore $`\alpha _\text{s}^2L`$ in the distribution. Accordingly the difference should be independent of $`\rho `$ for small $`\rho `$. One can see that this is the case for each individual colour factor. We also show that without the contribution $`𝒮`$ calculated in this paper the $`C_FC_A`$ part of the result is not consistent with Event2.

Using eq. (3) it is possible to combine our results for the single-hemisphere mass together with those of  for the heavy-jet mass in order to obtain predictions correct to single-log level for the light-jet mass (including the appropriate $`C_1`$ terms).

We point out that the form computed for $`𝒮`$ applies to all single-hemisphere 2-jet event-shapes that for large-angle gluons have a value of the order of the transverse momentum: in $`e^+e^{}`$ this means the hemisphere broadening (for the constant terms a relation analogous to (21) applies, with $`C_1^\tau `$ replaced with $`C_1^{B_T}`$). In DIS our form for $`𝒮`$ applies to the current hemisphere jet-mass, the $`C`$-parameter, and the thrust and broadening with respect to the thrust axis. The detailed phenomenology of the DIS variables will be considered elsewhere . It indicates that the inclusion of $`𝒮`$ reduces the peak height (in $`d\sigma /d\rho `$) by about $`30\%`$ before matching to the $`𝒪\left(\alpha _\text{s}^2\right)`$ calculation, and by about $`20\%`$ after matching (using an $`R`$-type matching ).

For the narrow-jet $`K_{out}`$ and jet masses in multi-jet events (including hadron-hadron event shapes) different forms will apply because of the more complicated geometry of the underlying hard event, but the Monte Carlo algorithm described here can be adapted to those cases as well.

Acknowledgments

We wish to thank Stefano Catani for helpful discussions and Andrea Banfi, Yuri Dokshitzer, Pino Marchesini and Giulia Zanderighi for helpful comments and suggestions. One of us (GPS) would like to acknowledge the hospitality of the DESY theory group while part of this work was carried out. We are also grateful to the INFN, Sezione di Milano and the Università di Milano Bicocca for the use of computing facilities.

## Appendix

Below we give details of the Monte Carlo algorithm used to calculate $`𝒮`$.666We present a version of the algorithm which has weighted events; there exists also a slightly more efficient version of the algorithm with unweighted events, but it requires some integrations to be done analytically and is therefore less easily generalisable to observables with more complex boundaries in phase space such as jet masses in multi-jet events.

For each dipole $`𝒟_i`$, the collinear angular cutoff $`ϵ`$ introduced in section 3 defines a rapidity interval $`\mathrm{\Delta }\eta _i`$ for emission in the dipole c.o.m. frame. We shall denote by $`𝒟_1`$ the dipole which has as one of its sources the quark in $`_\mathrm{R}`$. The evolution variable used will be $`t`$, defined in (19). This will account for the running of the coupling.

One starts an event with a quark-antiquark dipole (labelled $`ab`$), $`t=0`$ and a weight $`w=1`$. One then goes through the following steps to determine an event’s contribution to the quantity $`d𝒮/dt`$:

* Establish $`\mathrm{\Delta }\eta _{\mathrm{tot}}=_i\mathrm{\Delta }\eta _i`$. Increase $`t`$ by an amount $`\mathrm{\Delta }t`$ chosen according to a distribution proportional to

$$\mathrm{exp}\left(2C_A\mathrm{\Delta }\eta _{\mathrm{tot}}\mathrm{\Delta }t\right).$$
(23)
* Choose a dipole at random, such that the probability of obtaining $`𝒟_i`$ is $`\mathrm{\Delta }\eta _i/\mathrm{\Delta }\eta _{\mathrm{tot}}`$. Create a gluon uniformally in rapidity and azimuth in the dipole centre of mass frame, such that in the lab frame its angle relative to either of the dipole sources is larger than $`ϵ`$.
* There are now various possibilities:

+ The gluon is in $`_\mathrm{L}`$: insert it into dipole $`i`$ so as to create two new dipoles (removing the old one from the list). Go to step 1.
    + The gluon is in $`_\mathrm{R}`$ and was not emitted from $`𝒟_1`$: add $`w`$ to the bin corresponding to the value of $`t`$ and start a new event.
    + The gluon is in $`_\mathrm{R}`$ and was emitted from $`𝒟_1`$: denote the momenta of its sources by $`p_{\mathrm{}}`$ and $`p_a`$ (the latter is just the $`_\mathrm{R}`$ quark momentum) and the new gluon momentum by $`p_g`$. One determines the value of

    $$X=1\frac{F_{ab}}{F_a\mathrm{}}=1\frac{(ab)}{(gb)}\frac{(g\mathrm{})}{(a\mathrm{})}.$$

    If $`X>0`$ then with probability $`X`$ we add $`w`$ to current $`t`$-bin and start a new event; with probability $`1X`$ we throw away the gluon that has just been generated and go to step 1. If $`X<0`$ then we add $`Xw`$ to the current $`t`$ bin, multiply $`w`$ by $`1X`$, and go to step 1. This procedure is a simple extension to negative probabilities of the veto algorithm and corresponds to replacing $`F_𝒞`$, used implicitly in the previous steps, with $`\stackrel{~}{F}_𝒞`$.


text:
# Status of Neutrino Oscillations

## 1. Introduction

Undoubtedly the solar (Suzuki 2000) and atmospheric (Sobel 2000, Becker-Szendy 1992) neutrino problems provide the two most important milestones indicating physics beyond the Standard Model (SM). Of particular importance has been the confirmation in 1998 by the Super-Kamiokande (SK, for short) collaboration of the zenith-angle-dependent deficit of atmospheric neutrinos. Altogether the data provide a strong evidence for $`\nu _e`$ and $`\nu _\mu `$ conversions, respectively. Neutrino conversions are a natural consequence of theories beyond the Standard Model (Valle 1991). The first example is oscillations of low-mass neutrinos. While the theoretical understanding of the origin of neutrino masses is still lacking, there is a variety of attractive options available. Most likely, the exceptional nature of neutrinos as the only electrically neutral fermions in the SM underlies the smallness of their mass, as it would be associated with the violation of lepton number. Indeed in gauge theories one expects, on fundamental grounds, neutrinos to be Majorana fermions (Schechter 1980a). This is the generic situation in actual models. It will be surprising indeed if massive neutrinos turn out to be Dirac particles, like the quarks. Lepton number violation would imply processes such as neutrino-less double beta decay (Schechter 1982), novel CP violation effects (Schechter 1980a and 1981a), and/or neutrino electromagnetic properties (Schechter 1981b), so far unobserved. Present data and theoretical considerations suggest either hierarchical or quasi-degenerate neutrino masses. While solar neutrino rates favour the small mixing angle MSW oscillations (Wolfenstein 1978, Smirnov 1986), present data on the recoil-electron spectrum prefer the large mixing solutions. When interpreted in terms of neutrino oscillations, the observed atmospheric neutrino zenith-angle-dependent deficit clearly indicates that the mixing involved is maximal (Gonzalez-Garcia 2001). Adding information from reactor experiments one concludes that the third angle amongst the three neutrinos is small (Apollonio 1999). Thus, altogether, we have the intriguing possibility that, unlike the case of quarks, neutrino mixing is bi-maximal (Barger 1998, Davidson 1998, de Gouvea 2000, Chankowski 2000fp, Hirsch 2000) which could be tested at the upcoming long-baseline experiments or at a neutrino factory experiment (Quigg 1999) or at the proposed KamLAND experiment (De Braeckeleer 2000).

In addition to the above, there is also a long history of searches for neutrino oscillations at accelerators. Except for the unconfirmed hint provided by the LSND experiment (Athanassopoulos 1998, Smith 2000), accelerator searches have so–far been negative. The resulting limits, however, are not very restrictive on the scale of the indications from underground experiments and I will not discuss them any further. Barring exotic neutrino conversion mechanisms the hint of the LSND experiment together with the solar and atmospheric data require three mass scales, hence the need for a fourth light neutrino, which must be sterile (Peltoniemi 1993, Caldwell 1993, Liu 1998, Hirsch 2000, Giunti 2000). The most attractive possibility is to have, out of the four neutrinos, two of them lie at the solar neutrino scale, with the other two maximally-mixed neutrinos at the LSND scale (Peltoniemi 1993, Caldwell 1993, Hirsch 2000). These schemes have distinct implications at future solar & atmospheric neutrino experiments with sensitivity to neutral current neutrino interactions such as SNO. Cosmology can also place restrictions on these four-neutrino schemes (Raffelt 1999).

## 2. Indications for New Physics

The most solid hints in favour of new physics in the neutrino sector come from underground experiments on solar (Suzuki 2000) and atmospheric (Sobel 2000, Becker-Szendy 1992) neutrinos. The most recent SK data correspond to 1117–day solar and 1144–day (71 kton-yr) atmospheric data samples, respectively. There are also new data from Soudan-2 (5.1 kton-yr) and MACRO.

### 2.1. Solar Neutrinos

Our sun produces $`\nu _e`$ ’s through various nuclear reactions which take place in its interior. The predicted spectrum of solar neutrinos is illustrated in Fig. (1), taken from (Bahcall 1998). I will refer to this model as “the” SSM.

Solar neutrinos are detected either with geochemical methods (the $`\nu _e+^{37}Cl^{37}Ar+e^{}`$ reaction at the Homestake experiment and the $`\nu _e+^{71}Ga^{71}Ge+e^{}`$ reaction at the Gallex, Sage and GNO experiments) or through $`\nu _ee^{}`$ scattering on water, using Cerenkov techniques at Kamiokande and Super-Kamiokande. As summarized in Fig. (2) all experiments observe a deficit of 30 to 60 % whose energy dependence follows mainly from the lower Chlorine rate. Note that Fig. (2) includes the latest results from SK, SAGE & GNO presented at $`\nu `$2000, but not the first results from SNO.

It is convenient to present the predictions of various standard solar models in terms of the 7Be and 8B neutrino fluxes, normalized to the SSM predictions (Bahcall 1998), as seen in Fig. (3), which includes most of the existing solar models.

On the other hand the values of the fluxes indicated by measured neutrino event rates are shown by the contours in the lower-left part of the figure, with a negative best-fit 7Be neutrino flux! This discrepancy strongly suggests the need for new particle physics (Bahcall 1994). Since possible non-standard astrophysical solutions are rather constrained by helioseismology studies (Bahcall 1998) one is led to assume the existence of neutrino conversions, such as those induced by very small neutrino masses.

The high statistics of SK after 1117 days of data–taking also provides very useful information on the recoil electron energy spectrum with event rates given for 18 bins starting at 5.5 MeV ***They have also reported results of a lower energy bin 5 MeV $`<E_e<`$5.5 MeV, but due to systematic errors this is not yet included in the analysis..

The spectrum in Fig. (4) is well described by the flat hypothesis $`\chi _{flat}^2=13/(17\mathrm{d}\mathrm{o}\mathrm{f})`$, in contrast with hints from previous 825–day sample. Moreover, SK measures the zenith angle distribution (day/night effect) which is sensitive to the effect of the Earth matter in the neutrino propagation.

One sees a slight excess of events at night, but the corresponding day versus night asymmetry $`A_{D/N}=\frac{DN}{\frac{D+N}{2}}=0.034\pm 0.022\pm 0.013`$ is only 1.3$`\sigma `$ away from zero. In order to combine this day–night information with the spectral data, SK has also presented separately the measured recoil energy spectrum during the day and during the night. This will be referred in the following as the day–night spectra data and contains $`2\times 18`$ data bins.

Note that the absence of clear hints of spectral distortion, day-night or seasonal variation implies that, per se, they do not give any clear indication for physics beyond the standard model. From this point of view, despite the increasing weight of such rate-independent observables, the solar neutrino problem rests heavily on the rate discrepancy. Nevertheless, as we will see, rate-independent observables are already playing an important rôle in selecting amongst different solutions of the solar neutrino problem.

### 2.2. Atmospheric Neutrinos

Neutrinos produced as decay products in hadronic showers from cosmic ray collisions with nuclei in the upper atmosphere have been observed in several experiments (Sobel 2000, Becker-Szendy 1992). Although individual $`\nu _\mu `$ or $`\nu _e`$ fluxes are only known to within $`30\%`$ accuracy, their ratio is predicted to $`5\%`$ over energies that vary from 0.1 GeV to 100 GeV (Gaisser 1998). The long-standing discrepancy between the predicted and measured $`\mu /e`$ ratio of the muon ($`\nu _\mu +\overline{\nu }_\mu `$) over the electron atmospheric neutrino flux ($`\nu _e+\overline{\nu }_e`$) (Gaisser 1995) found both in water Cerenkov experiments (Kamiokande, SK and IMB) as well as in the iron calorimeter Soudan2 experiment is illustrated in Fig. (6)

This evidence has now been strengthened by the fact that it exhibits a strong zenith-angle dependence (Sobel 2000) as can be seen from Fig. (7). The zenith-angle distributions for the Super-Kamiokande $`e`$–like are shown in the left panels, both in the sub-GeV (upper panels) and multi-GeV (lower panels) energy range. The thick solid line is the expected distribution in the SM. It is consistent with the SM expectations. In contrast $`\mu `$–like events displayed in the right-panels show a clear deficit of neutrinos coming from below, which is very suggestive indeed of $`\nu _\mu `$ oscillations. In Fig. (7) we also give the predicted best-fit distributions obtained in a global 3-neutrino oscillation description of the data (Gonzalez-Garcia 2001). The thin full line is the prediction for the overall best-fit point of the contained plus up-going atmospheric data sample, with $`\mathrm{tan}^2\theta _{13}=0.025`$, $`\mathrm{\Delta }m_{32}^2=3.3\times 10^3`$ eV2 and $`\mathrm{tan}^2\theta _{23}=1.6`$.

Zenith-angle distributions have also been recorded for upward-going muon events in Super-Kamiokande and MACRO, as illustrated in Fig. (8).

The thick solid line is the expected distribution in the SM, while the thin full line corresponds, as in Fig. (7), to the prediction for the overall global best-fit point.

## 3. Three–neutrino fits

The most economical joint description of solar and atmospheric anomalies involves oscillations amongst all three known types of neutrinos. Here I summarize the results of the recent global analysis (Gonzalez-Garcia 2001) of the solar (Suzuki 2000), atmospheric (Sobel 2000, Becker-Szendy 1992) and reactor (Apollonio 1999) neutrino data in terms of three–neutrino oscillations. The present discussion goes beyond previous three–neutrino oscillation analyses including only solar (Fogli 2000) or only atmospheric neutrino data (Fogli 1999). It also updates some joint studies (Barger 1980) which do not take into account the latest SK data. As we saw the most recent solar neutrino rates include the 1117 day SK data sample on the recoil electron energy spectra for day and night periods. On the other hand the atmospheric data sample includes not only contained events but also the upward-going $`\nu `$-induced muon fluxes. In addition to previous Frejus, IMB, Nusex, and Kamioka data we use the most recent 71 kton-yr (1144 days) SK data set, the 5.1 kton-yr contained events of Soudan2, and the results on up–going muons from the MACRO detector.

The pattern of neutrino oscillations expected in any fundamental (gauge) theory of neutrino mass is determined by the structure of the lepton mixing matrix (Valle 1991, Bilenkii 1999). For the simplest three–neutrino theories this is in general characterized by three mixing angles and three CP violating phases (Schechter 1980a). The latter include, in addition to the Dirac-type phase analogous to that of the quark sector, two extra physical (Schechter 1981a) phases associated to the Majorana character of neutrinos. Conservation of CP implies that Dirac phases are zero modulo $`\pi `$, while Majorana phases are zero modulo$`\pi /2`$ (Wolfenstein 1981, Schechter 1981b). For our following discussion all three phases are set to zero. In this case the mixing matrix can be conveniently chosen in the form (Schechter 1980a)

$$\left(\begin{array}{ccc}c_{13}c_{12}& s_{12}c_{13}& s_{13}\\ s_{12}c_{23}s_{23}s_{13}c_{12}& c_{23}c_{12}s_{23}s_{13}s_{12}& s_{23}c_{13}\\ s_{23}s_{12}s_{13}c_{23}c_{12}& s_{23}c_{12}s_{13}s_{12}c_{23}& c_{23}c_{13}\end{array}\right)$$
(1)
The joint study of solar and atmospheric neutrino oscillations is characterized by a five-dimensional parameter space

$$\begin{array}{cc}\mathrm{\Delta }m_{}^2\hfill & \mathrm{\Delta }m_{21}^2=m_2^2m_1^2\hfill \\ \mathrm{\Delta }m_{atm}^2\hfill & \mathrm{\Delta }m_{32}^2=m_3^2m_2^2\hfill \\ \theta _{}\hfill & \theta _{12}\hfill \\ \theta _{atm}\hfill & \theta _{23}\hfill \\ \theta _{reactor}\hfill & \theta _{13}\hfill \end{array}$$
(2)
where all mixing angles are assumed to lie in the full range from $`[0,\pi /2]`$.

From the required hierarchy in the splittings $`\mathrm{\Delta }m_{atm}^2\mathrm{\Delta }m_{}^2`$ indicated by the solutions to the solar and atmospheric neutrino anomalies (see below) it follows that the analyses of solar data constrain three of the five independent oscillation parameters, namely, $`\mathrm{\Delta }m_{21}^2,\theta _{12}`$ and $`\theta _{13}`$ since for most cases oscillations over the atmospheric scale average out. Conversely, from the point of view of the atmospheric data analysis one can effectively assume that the lighter neutrinos become degenerate so that one can rotate away the corresponding angle $`\theta _{12}`$. The leptonic mixing matrix takes on the simplified form (Schechter 1980b)

$$𝐑=\left(\begin{array}{ccc}c_{13}& 0& s_{13}\\ s_{23}s_{13}& c_{23}& s_{23}c_{13}\\ s_{13}c_{23}& s_{23}& c_{23}c_{13}\end{array}\right);$$
(3)
As a result only three oscillation parameters: $`\mathrm{\Delta }m_{32}^2`$, $`\theta _{23}`$ and $`\theta _{13}`$ are necessary to describe the 3-neutrino propagation of atmospheric neutrinos.

It follows from the above discussion that $`\theta _{13}`$ is the only parameter common to both analyses and may potentially allow for some “cross-talk” between the two sectors. It is known that for $`\mathrm{\Delta }m_{32}^2\mathrm{\Delta }m_{12}^2`$ and $`\theta _{13}=0`$ the atmospheric and solar neutrino oscillations decouple in two 2–$`\nu `$ oscillation scenarios. In this respect our results also contain as limiting cases the pure two–neutrino oscillation scenarios and update previous analyses on atmospheric neutrinos (Fornengo 2000, Foot 1998) and solar neutrinos (Gonzalez-Garcia 2000a, Bahcall 1998).

In order to compute the solar neutrinos survival probabilities for any value of the neutrino mass and mixing the full expression for the survival probability has been used, without appealing to the usual approximations whose validity defines the MSW (Smirnov 1986) or the “ just-so”(Glashow 1987) regime. The treatment of neutrino oscillations is therefore unified, with MSW and vacuum oscillations considered on the same footing. Likewise, we include in our description conversions with $`\theta _{12}>\pi /4`$ (Gonzalez-Garcia 2000b). Results are found by numerically solving the Schrodinger neutrino evolution equation in the Sun and the Earth matter, using the electron number density of BP2000 model (Bahcall 2000) and the Earth density profile given in the Preliminary Reference Earth Model (PREM) (Dziewonski, 1981).

Reactor limits

To start the summary of the global 3-neutrino fits we first note that, of all laboratory searches for neutrino oscillation, reactors provide the most sensitive one when it comes to comparing with the strong hints from underground experiments. Thus we will first consider the region of oscillation parameters which can be excluded from reactor experiments. The restrictions on $`\mathrm{\Delta }m_{32}^2`$ and $`\mathrm{sin}^2(2\theta _{13})`$ that follow from the non observation of oscillations at the CHOOZ reactor experiment are shown in Fig. (9), taken from (Gonzalez-Garcia 2001).

The curves represent the 90, 95 and 99% CL excluded region defined with 2 d. o. f. for comparison with the CHOOZ published results. In what follows we will compare this direct bound with those obtained from a global analysis of solar and atmospheric data, as well as consider its effect in combination with the latter.

Solar data fit

We first present the allowed regions of solar oscillation parameters $`\theta _{12}`$, $`\mathrm{\Delta }m_{21}^2`$ as a function of $`\theta _{13}`$. All plots are taken from (Gonzalez-Garcia 2001) where all details can be found. In Fig. (10) we give the allowed three–neutrino oscillation regions in $`\mathrm{\Delta }m_{21}^2`$ and $`\mathrm{tan}^2\theta _{12}`$ from the measurements of the total event rates at the Chlorine, Gallium, Kamiokande and Super-Kamiokande experiments. The different panels represent the allowed regions at 99% (darker) and 90% CL (lighter) obtained as sections for fixed values of the mixing angle $`\mathrm{tan}^2\theta _{13}`$ of the three–dimensional volume defined by $`\chi ^2\chi _{min}^2`$=6.25 (90%), 11.36 (99%). The best–fit point is denoted as a star. It occurs, as expected, for the small mixing MSW solution (SMA) simply because this is the situation which most strongly suppresses the unwanted 7Be neutrino flux. It is characterized by a non-zero $`\theta _{13}`$ value. For higher $`\theta _{13}`$, the description worsens. Although relatively weak, the limit on $`\theta _{13}`$ from solar data is totally independent on the allowed range of the atmospheric mass difference $`\mathrm{\Delta }m_{32}^2`$.

The three–neutrino solar oscillation regions excluded by the measurement of the day–night spectra data in the Super-Kamiokande 1117-day data sample is illustrated in Fig. (11).

Finally, the allowed three–neutrino solar oscillation regions in $`\mathrm{\Delta }m_{21}^2`$ and $`\mathrm{tan}^2\theta _{12}`$ which follows from the global analysis of solar neutrino data is presented in Fig. (12). The best–fit point is denoted as a star.

The relative quality of the various oscillation solutions of the solar neutrino problem is illustrated in Fig. (13).

This figure gives $`\mathrm{\Delta }\chi ^2`$ as a function of $`\mathrm{tan}^2\theta _{13}`$ from the 3–neutrino analysis of the solar data. The dotted horizontal lines correspond to the 90%, 99% CL limits. The left panel corresponds the analysis of total rates only, while the right panel corresponds to the global analysis. The dotted horizontal lines correspond to the 90%, 99% CL limits. Though all the various SMA, LMA, LOW and vacuum solutions are still acceptable global descriptions of the solar data, they are not equally good. From the right panel in Fig. (13) we conclude that, for small enough $`\mathrm{tan}^2\theta _{13}`$ values (so that we gets effectively a two-neutrino scheme) the best solution is the LMA solution is the best, while SMA is the worst. One notices also that rate-independent observables, such as the electron recoil energy day-night spectra, are playing an increasing role in discriminating between different solutions to the solar neutrino problem, pushing the best–fit point towards the LMA solution (Gonzalez-Garcia 2001), a trend already noted in earlier 2-neutrino analyses due to the same reason (Gonzalez-Garcia 2000a).

Another issue which could play a more significant role in future investigations is that of seasonal variations, expected in the just-so regime and, more subtly, also in the MSW large mixing solutions, LMA and LOW. The first would result from the eccentricity of the Earth’s orbit around the Sun and could be tested by searching for seasonal variations in the 7Be neutrino flux at the Borexino experiment, and possibly at KamLAND (de Gouvea 1999). The latter would result from the regeneration effect at the Earth (day-night effect) and might be tested through time variations of event rates at GNO and Borexino (de Holanda 1999).

An interesting theoretical issue is the possible effect of random fluctuations in the solar matter density (Balantekin 1996, Nunokawa 1996, Bamert 1998) on the solar neutrino event rates. The existence of such noise fluctuations at a few percent level is not excluded by present helioseismology studies. The correlation length $`L_0`$ associated with the scale of the fluctuation can be assumed to lie between the mean free path of the electrons in the solar medium, $`l_{free}10`$ cm, and the neutrino oscillation length in matter, $`\lambda _m`$, e. g. $`l_{free}L_0\lambda _m`$. Even small fluctuations can have an important effect on averaged solar neutrino survival probabilities, especially for small solar mixing angles (Balantekin 1996, Nunokawa 1996). The fluctuations can affect the 7Be neutrino component of the solar neutrino spectrum, implying that Borexino can probe, at some level of precision, the magnitude of solar matter density fluctuations, thus an additional motivation for the experiment (Arpesella 1991).

Atmospheric and reactor data fit

Here I present the allowed regions of atmospheric oscillation parameters $`\theta _{23}`$, $`\mathrm{\Delta }m_{32}^2`$ for different values of $`\theta _{13}`$, common to solar and atmospheric analyses. All plots are taken from (Gonzalez-Garcia 2001) where details can be found. In Fig. (14) we display the allowed $`(\mathrm{tan}^2\theta _{23},\mathrm{\Delta }m_{32}^2)`$ regions for different $`\mathrm{tan}^2\theta _{13}`$ values, that follow from the combination of Super-Kamiokande atmospheric neutrino events. The regions refer to 90, 95 and 99% CL. The best–fit point is denoted as a star and corresponds to $`\mathrm{tan}^2\theta _{13}=0.026`$.

In Fig. (15) we present the three–neutrino regions in $`(\mathrm{tan}^2\theta _{23},\mathrm{\Delta }m_{32}^2)`$ allowed by the combination of all atmospheric neutrino data plus Chooz. The best–fit point is denoted as a star.

Comparing Fig. (14) and Fig. (15) one can see the weight of the reactor neutrino bound on the global analysis. The best global fit has $`\mathrm{tan}^2\theta _{13}=0.005`$ and for $`\mathrm{tan}^2\theta _{13}=0.075`$ even the 99% CL allowed disappears. Notice the important complementarity between atmospheric data and the reactor limits on $`\theta _{13}`$ since the latter apply only for $`\mathrm{\Delta }m_{32}^2\stackrel{>}{}\mathrm{\hspace{0.25em}10}^3`$ eV2.

Combining solar, atmospheric and reactor data

We also obtain the allowed ranges of parameters from the full five–dimensional combined analysis of all of the above neutrino data. In Fig. (16) we give the regions in $`\mathrm{\Delta }m_{21}^2`$ and $`\mathrm{tan}^2\theta _{12}`$ allowed by the global analysis of solar, atmospheric and reactor neutrino data. The left panel gives the regions for the unconstrained analysis defined in terms of the increases of $`\mathrm{\Delta }\chi ^2`$ for 5 d.o.f. from the global best fit point denoted as a star. The right panel shows the values of $`\mathrm{tan}^2\theta _{13}`$ beyond which the 99% CL region starts to disappear.

## 4. Putting the pieces together

Together with the solar neutrino data, the angle-dependent atmospheric neutrino deficits provide a strong evidence for physics beyond the Standard Model. Small neutrino masses provide the simplest and most generic explanation of the data. Theoretical neutrino mass models fall into two classes, which I call bottom-up and top-bottom. They can lead to either hierarchical (de Gouvea 2000, Hirsch 2000) or quasi-degenerate (Chankowski 2000) neutrinos. Top-bottom approaches inspired on the idea of Unification typically employ either a see-saw mechanism, or high dimension operators. As examples of minimalistic models I mention two we have recently investigated, both of which require a large mixing solution to the solar neutrino problem (de Gouvea 2000, Chankowski 2000).

Supersymmetry with bilinear breaking of R–parity provides a simple bottom-up-type model which allows for the possibility of probing the neutrino mixing, as indicated by the underground experiments, within the context of high–energy collider experiments such as the LHC (Hirsch 2000). For additional models and for models involving specific Yukawa textures see (Davidson 1998, Lola 1998, Altarelli 1998).

Last, but not least, I mention that it is impossible to reconcile solar, atmospheric and reactor data with the LSND hint in terms of neutrino oscillations without a light sterile neutrino (Peltoniemi 1993, Caldwell 1993) and in ref. (Hirsch 2000) we give an updated discussion in the context of an interesting model. For more references see (Liu 1998) and for a more complete discussion of limits on 4-neutrino models see (Giunti 2000). If light sterile neutrinos exist they should be probed at neutral-current-sensitive solar & atmospheric neutrino experiments such as SNO (Gonzalez-Garcia 2000c).

## 5. Conclusion

Within the neutrino oscillation framework present solar and atmospheric data suggest the intriguing possibility of bi-maximal neutrino mixing, which explicitly illustrates the sharp contrast between the lepton and quark sectors of the theory. With good luck this could be checked on the one hand at the upcoming long-baseline experiments or at a neutrino factory (Quigg 1999) and, on the other, via the search of seasonal effects in solar neutrinos, e. g. at the proposed KamLAND experiment (De Braeckeleer 2000). In certain models (Hirsch 2000) one may test the neutrino mixing angles involved in the explanation of the neutrino anomalies at high–energy collider experiments (Porod 2000), illustrating an amusing synergy between accelerator and underground experiments. On the other hand if the LSND result stands the test of time, this would be a strong indication for the existence of a light sterile neutrino which would be another radical difference between leptons and quarks. Let me mention however that the neutrino oscillation interpretation of solar neutrino anomalies is still far from unique. For example, resonant spin flavor precession (Akhmedov 1988) induced by transition magnetic moments (Schechter 1981b). For a recent fit of solar data see ref. (Miranda 2001). From this point of view, it is still too early for a precision determination of neutrino properties from underground experiments. Last, but not least, let me mention that in fifty years of weak interaction physics an answer to the most fundamental question about the nature - Dirac versus Majorana - of neutrinos has so far defied all experimental attempts.

Acknowledgements:

This work was supported by Spanish DGICYT under grant PB98-0693, by the European Commission TMR network HPRN-CT-2000-00148 and by the European Science Foundation network grant N. 86.

## References

Akhmedov, E. (1988) Phys. Lett. B213 64-68; C. S. Lim and W. Marciano, Phys. Rev. D37 1368

Altarelli, G. and Feruglio, F. (1999) Phys. Lett. B451 (1999) 388; S. Lola and G.G. Ross, hep-ph/9902283; R. Barbieri, L.J. Hall and A. Strumia, Phys. Lett. B445, 407 (1999); M.E. Gomez, G.K. Leontaris, S. Lola and J.D. Vergados, Phys. Rev. D59, 116009 (1999); G.K. Leontaris, S. Lola, C. Scheich and J.D. Vergados, Phys. Rev. D53, 6381 (1996)

Apollonio, M. et al., Phys.Lett. B466 (1999) 415; F. Boehm et al., hep-ex/9912050

Arpesella, C. at al., Proposal of the Borexino experiment (1991).

Athanassopoulos, C. et al \[LSND Collaboration\], Phys. Rev. Lett. 75 (95) 2650; Phys. Rev. Lett.  77 (96) 3082 ; C. Athanassopoulos et al, Phys. Rev. Lett. 81, 1774 (1998)

Bahcall, J., Basu, S. and Pinsonneault, M. (2000) Astrophys. J. 529, 1084 (2000).

Bahcall, J.N., Krastev, P.I., Smirnov, A.Yu. (1998) Phys. Rev. D58, 09601 (1998); S. Goswami, D. Majumdar, A. Raychaudhuri, hep-ph/0003163.

Bahcall, J. N. and Pinsonneault, M. (1998) http://www.sns.ias.edu/ jnb/ and Phys. Lett. B 433 (1998) 1.

Bahcall, J. N., (1994) Phys. Lett. B338 (1994) 276; V. Castellani, et al Phys. Lett. B324 (1994) 245; N. Hata, S. Bludman, and P. Langacker, Phys. Rev D 49 (1994)3622; V. Berezinsky, Comm. on Nucl. and Part. Phys. 21 (1994) 249

Balantekin, A.B et al (1996) Phys. Rev. D54 (1996) 3941; F. N. Loreti and A. B. Balantekin, Phys. Rev. D50 (1994) 4762; F. N. Loreti et al., Phys. Rev. D52 (1996) 6664.

Bamert, P., Burgess, C. and Michaud, D. (1998) Nucl. Phys. B513, 319 (1998); C.P. Burgess, hep-ph/9711425; C.P. Burgess and D. Michaud, Annals Phys. 256, 1 (1997) and hep-ph/9611368.

Barger, V., Whisnant, K., Phillips, R. J. N. (1980) Phys. Rev. D22, 1636 (1980); G.L. Fogli, E. Lisi, D. Montanino, Phys. Rev. D49, 3626 (1994); Astropart. Phys. 4, 177 (1995); R. Barbieri, L. J. Hall, D. Smith, A. Strumia, N. Weiner, JHEP 9812, 017 (1998); T. Teshima, T. Sakai, O. Inagaki, Int. J. Mod. Phys. A14, 1953 (1999); V. Barger, K. Whisnant, Phys. Rev. D59, 093007 (1999).

Barger, V., Whisnant, K. and Phillips, R. J. (1981) Phys. Rev. D24, 538 (1981); S.L. Glashow and L.M. Krauss, Phys. Lett. 190B, 199 (1987); S.L. Glashow, P.J. Kernan and L.M. Krauss, Phys. Lett. B445, 412 (1999)

Barger, V. et al Phys. Lett. B437, 107 (1998); E. Torrente-Lujan, Phys. Lett. B389 (1996) 557

Becker-Szendy, R. et al. (1992), Phys. Rev. D46, 3720 (1992). H.S. Hirata et al., Phys. Lett. B280, 146 (1992); Y. Fukuda et al., ibid B335, 237 (1994); W.W.M. Allison et al., Phys. Lett. B449, 137 (1999); A. Mann, http://nu2000.sno.laurentian.ca/; F. Ronga, hep-ex/0001058 and http://www.ichep2000.rl.ac.uk/; K. Daum et al. Z. Phys. C66, 417 (1995); M. Aglietta et al., Europhys. Lett. 8, 611 (1989).

Bilenkii, S. Giunti,C. and Grimus, W. (1999) Prog. Part. Nucl. Phys. 43 (1999) 100

Caldwell D. O. and Mohapatra, R. N. Phys. Rev. D48 (1993) 329

Chankowski, P. et al, Phys. Rev. Lett. 86, (2001) 3488 \[hep-ph/0011150\].

Davidson, S. and King, S. Phys. Lett. B445, 191 (1998); S. F. King, Nucl. Phys. B576 (2000) 85

De Braeckeleer, L. Nucl. Phys. Proc. Suppl. 87, 312 (2000).

de Gouvea, A. and Valle, J. W. F. (2000) Phys. Lett. B 501 (2001) 115 \[hep-ph/0010299\].

de Gouvea, A., Friedland, A. and Murayama, H. (1999) Phys. Rev. D60, 093011 (1999).

de Holanda, P. C. et al (1999) Phys. Rev. D60 (1999) 093010, hep-ph/9903473; G. L. Fogli, E. Lisi, D. Montanino and A. Palazzo, Phys. Rev. D61 (2000) 073009; see also J.N. Bahcall, P.I. Krastev and A.Y. Smirnov, hep-ph/9905220.

Dziewonski, A.M. and Anderson, D. Phys. Earth Planet. Inter. 25, 297 (1981).

Fogli, G. L. et al (1999) Phys. Rev. D59, 033001 (1999); G.L. Fogli, E. Lisi, D. Montanino and G. Scioscia Phys. Rev. D55, 4385 (1997); A. De Rujula, M.B. Gavela, P. Hernandez, hep-ph/0001124; T. Teshima, T. Sakai, Prog. Theor. Phys. 101, 147 (1999); Prog. Theor. Phys. 102, 629 (1999); hep-ph/0003038.

Fogli, G. L. et al (2000) Phys. Rev. D62, 013002 (2000); Phys. Rev. D54, 2048 (1996); A.M. Gago, H. Nunokawa, R. Zukanovich Funchal, hep-ph/0007270.

Foot, R. Volkas, R. and Yasuda, O. (1998) Phys. Rev. D58, 013006 (1998); O. Yasuda, Phys. Rev. D58, 091301 (1998); G.L. Fogli, E.Kh. Akhmedov, A. Dighe, P. Lipari and A.Yu. Smirnov, hep-ph/9808270.

Fornengo, N, Gonzalez-Garcia;M. and Valle, J. W. F. (2000) Nucl. Phys. B580 (2000) 58; M.C. Gonzalez-Garcia, H. Nunokawa, O.L. Peres and J. W. F. Valle, Nucl. Phys. B543, 3 (1999) and M. C. Gonzalez-Garcia, H. Nunokawa, O. L. Peres, T. Stanev and J. W. F. Valle, Phys. Rev. D58 (1998) 033004.

Gaisser, T. and Stanev, T. (1998) Phys. Rev. D57 1977 (1998); G. Barr, T. K. Gaisser and T. Stanev, Phys. Rev. D 39 (1989) 3532 and Phys. Rev. D38, 85; V. Agrawal et al., Phys. Rev. D53, 1314 (1996); L. V. Volkova, Sov. J. Nucl. Phys. 31, 784 (1980); M. Honda, T. Kajita, S. Midorikawa and K. Kasahara, Phys. Rev. D52, 4985 (1995).

Gaisser, T. Halzen, F. and Stanev, T., (1995) Phys. Rep. 258, 174 (1995)

Giunti, C., Gonzalez-Garcia, M. and Peña-Garay, C. (2000) Phys. Rev. D62, 013005 (2000); Yasuda, O. hep-ph/0006319

Gonzalez-Garcia, M. C. et al, (2000a) Nucl. Phys. B573, 3 (2000).

Gonzalez-Garcia, M. C. and Pena-Garay, C. (2000b) Phys. Rev. D62 (2000) 031301; A. de Gouvea, A. Friedland and H. Murayama, Phys. Lett. B490 (2000) 125 \[hep-ph/0002064\]; G. L. Fogli, E. Lisi and D. Montanino, Phys. Rev. D54, 2048 (1996)

Gonzalez-Garcia, M. and Peña-Garay, (2000c) hep-ph/0011245

Gonzalez-Garcia, M. C., et al (2001) Phys. Rev. D 63 (2001) 033005 \[hep-ph/0009350\]

Hannestad, S. and Raffelt, G. (1999) Phys. Rev. D59, 043001 (1999)

Hirsch, M. and Valle, J. W. F. Phys. Lett. B 495 (2000) 121 \[hep-ph/0009066\].

Hirsch, M. et al, Phys. Rev. D 62 (2000) 113008 \[hep-ph/0004115\]; J. C. Romão, M. A. Diaz, M. Hirsch, W. Porod and J. W. F. Valle, Phys. Rev. D61, 071703 (2000) \[hep-ph/9907499\]

Liu, Q. Y. and Smirnov, A. Nucl. Phys. B524 (1998) 505-523; V. Barger, K. Whisnant and T. Weiler, Phys.Lett. B427 (1998) 97-104; S. Gibbons, R. N. Mohapatra, S. Nandi and A. Raichoudhuri, Phys. Lett. B430 (1998) 296-302; Nucl.Phys. B524 (1998) 505-523; S. Bilenky, C. Giunti and W. Grimus, Eur. Phys. J. C 1, 247 (1998); S. Goswami, Phys. Rev. D 55, 2931 (1997); N. Okada and O. Yasuda, Int. J. Mod. Phys. A12 (1997) 3669-3694 E. J. Chun, A. Joshipura and A. Smirnov, in Elementary Particle Physics: Present and Future (World Scientific, 1996), ISBN 981-02-2554-7; P. Langacker, Phys. Rev. D58, 093017 (1998); M. Bando and K. Yoshioka, Prog. Theor. Phys. 100, 1239 (1998)

Lola, S. and Vergados, J. (1998) Prog. Part. Nucl. Phys. 40, 71 (1998); G. Altarelli and F. Feruglio, Phys. Lett. B439, 112 (1998).

Miranda, O. et al (2001) Nucl. Phys. B 595 (2001) 360 \[hep-ph/0005259\]. See also J. Pulido and E. Akhmedov, Phys. Lett. B485 (2000) 178; M. M. Guzzo and H. Nunokawa, Astropart. Phys. 12 (1999) 8; J. Derkaoui and Y. Tayalati, Astropart. Phys. 14 (2001) 351.

Nunokawa, H. et al (1996) Nucl. Phys. B472 495 \[see also talk at Neutrino 96, hep-ph/9610526\]

Peltoniemi, J. T. , Tommasini, D. and Valle, J. W. F. Phys. Lett. B298 (1993) 383; Peltoniemi, J. T. and Valle, J. W. F. Nucl. Phys. B406 (1993) 409 \[hep-ph/9302316\]; Ioannisian, A. and Valle, J. W. F. Phys. Rev. D 63 (2001) 073002 \[hep-ph/9911349\]

Porod, W. et al “Testing neutrino mixing at future collider experiments,” hep-ph/0011248, to appear in Phys. Rev. D.

Quigg, C. Introduction to NuFact ’99, the ICFA/ECFA Workshop on Neutrino Factories Based on Muon Storage Rings, hep-ph/9908357.

Schechter, J. and Valle, J. W. F. (1980a) Phys. Rev. D22, 2227 (1980).

Schechter, J. and Valle, J. W. F. (1980b) J. Schechter and J. W. F. Valle, Phys. Rev. D21 309 (1980)

Schechter, J. and Valle, J. W. F. (1981a) Phys. Rev. D23 (1981) 1666 showed that Majorana phases are physical, though their effects are small due to the smallness of neutrino masses.

Schechter, J. and Valle, J. W. F. (1981b) Phys. Rev. D24 1883, (1981), Err. ibid. D25 283, (1982).

Schechter, J. and Valle, J. W. F. (1982) Phys. Rev. D25, 2951 (1982); for a review on double beta decay experiments and comparison with theory see A. Morales, hep-ph/9809540

Smirnov A. and Mikheev, S. P. (1986) Neutrino Oscillations In Matter With Varying Density, In Proceedings, \*Tignes 1986, massive neutrinos\* 355

Smith, D. Proc. of the Third Dark Matter Conference DARK2000, Heidelberg, July 10-16, 2000; W. Louis, http://www.neutrino.lanl.gov/LSND/ W.C. Louis \[LSND Collaboration\], Prog. Part. Nucl. Phys. 40, 151 (1998).

Sobel, H. (2000) http://nu2000.sno.laurentian.ca; Y. Fukuda et al., Phys. Rev. Lett. 82 (1999) 2644 \[hep-ex/9812014\]; Phys. Lett. B433, 9 (1998); Phys. Lett. B436, 33 (1998); T. Toshito, http://www.ichep2000.rl.ac.uk

Suzuki, Y. (2000) Talk at XIX Int. Conference on Neutrino Physics and Astrophysics, Sudbury, Canada, http://nu2000.sno.laurentian.ca; T. Takeuchi, talk at the XXXth International Conference on High Energy Physics, ICHEP 2000 http://www.ichep2000.rl.ac.uk; Homestake Collaboration, B.T. Cleveland et al., Astrophys. J. 496, 505 (1998); R. Davis, Prog. Part. Nucl. Phys. 32, 13 (1994); K. Lande, http://nu2000.sno.laurentian.ca; SAGE Collaboration, J.N. Abdurashitov et al., Phys. Rev. C60, 055801 (1999); V. Gavrin, http://nu2000.sno.laurentian.ca. GALLEX Collaboration, W. Hampel et al., Phys. Lett. B447, 127 (1999). E. Belloti, http://nu2000.sno.laurentian.ca.

Valle, J. W. F. (1991) Gauge Theories and the Physics of Neutrino Mass, Prog. Part. Nucl. Phys. 26 (1991) 91

Wolfenstein, L. (1978) L. Wolfenstein, Phys. Rev. D 17 (1978) 236

Wolfenstein, L. (1981) Phys. Lett. B107 (1981) 77