text:
# Oscillation induced neutrino asymmetry growth in the early Universe

## 1 Introduction

Recent observations from atmospheric and solar neutrinos from Super Kamiokande (SK) and Sudbury Neutrino Observatory (SNO)  prefer conventional, mainly active-active mixings between ordinary electron, muon and tau neutrinos as the solutions to the observed flux deficits. However, while pure active-sterile solutions are disfavoured, sizable admixtures of sterile states in the observed fluxes are still allowed . New sterile states are also needed in order to explain the LSND anomaly  by neutrino physics. Moreover, in a recent analysis of the so called “2+2”-models111Alternative “3+1” models fit less well with short baseline reactor experiments . for four-neutrino mixing, the solutions with nonzero sterile neutrino components were found to provide best global fits to the solar, atmospheric and reactor data . Corresponding limits on active-sterile mixings are quite generous, for example

$`\mathrm{sin}^2\theta _{\mu s}`$ $`0.48`$ $`(\mathrm{Atmospheric})`$
$`\mathrm{sin}^2\theta _{es}`$ $`0.72`$ $`(\mathrm{Solar}\mathrm{LMA}),`$ (1)
where LMA refers to the Large Mixing Angle solution for the solar neutrino deficit.

Active-sterile mixing, if realized, would have several interesting effects in astrophysical settings  and in particular for the evolution of the early universe . In particular, sterile states could be brought into thermal equilibrium by mixing before nucleosynthesis, so that the resulting anomalous increase in the expansion rate of the universe would lead to overproduction of helium in disagreement with the observations . Excluding the parameter sets leading to equilibration provides useful bounds on active-sterile mixing. Indeed, from results of ref.  one can infer that the sterile components in mass eigenstates responsible for atmospheric anomaly and LMA are constrained by

$`\mathrm{sin}^2\theta _{\mu s}`$ $`0.013`$ $`(\mathrm{Atmospheric})`$
$`\mathrm{sin}^2\theta _{es}`$ $`0.026`$ $`(\mathrm{Solar}\mathrm{LMA})`$ (2)
where we used a SBBN-limit of $`N_\nu 3.4`$ for the number of effective neutrino degrees of freedom . These numbers correspond to the “light” case where the sterile component is the heavier of the mixing states; in the opposite, “dark” case, the corresponding limits are even stronger by one to two orders of magnitude. In either case the cosmological constraints are more stringent than those obtained in terrestrial laboratories.

While quite generic, the cosmological bounds (2) do depend on some simple prior assumptions, the most important of which is that the primordial lepton asymmetry is not anomalously large. It may be possible to circumvent them in more complicated mixing schemes involving more sterile states (for a recent discussion, see ). In a particular attempt it was shown in ref.  that for a large negative $`\delta m^2`$ and a small enough mixing angle (so that the equilibration bounds of  can be avoided), say on the $`\nu _\tau \nu _s^{}`$-sector, resonant oscillations may trigger a rapid growth of the tau lepton asymmetry. This asymmetry could then become very large, $`L_\tau 𝒪(1)`$ and change the natural prior condition of no significant initial lepton asymmetry for a mixing in any other sectors. Indeed, if created early enough, such an asymmetry could suppress $`\nu _\mu \nu _s`$-oscillations , obviating the bounds (2) for $`\nu _\mu \nu _s`$-mixing. However, as this scenario requires that the sterile state is the lighter one (resonance condition), and that the mass splitting $`|\delta m^2|`$ is large (to create $`L_\tau `$ early enough), it is essentially excluded by the recent solar and atmospheric neutrinos combined with the direct constraint on the electron neutrino mass .

Even with the mechanism of ref.  by and large excluded, there can be important effects due to $`L`$-growth with smaller $`|\delta m^2|`$. For example, a large homogeneous electron neutrino asymmetry $`L_{\nu _e}`$ would directly influence the nucleosynthesis prediction for the helium-4 abundance through the reactions

$`n+e^+`$ $``$ $`p+\overline{\nu }_e`$
$`p+e^{}`$ $``$ $`n+\nu _e,`$ (3)
which keep the neutron-to-proton ratio in equilibrium at early times. Changing the electron neutrino abundance could tilt the balance of these reactions with the possibility of either increasing $`Y_{He}`$ (negative $`L_{\nu _e}`$), and hence strengthening the bounds (2), or decreasing $`Y_{He}`$ (positive $`L_{\nu _e}`$), leading to weakening of the bounds (2).

The physics involved with the resonant growth of the asymmetry is very complicated because of several vastly different physical scales both in the temporal direction and in momentum variable, which effect the evolution of the ensemble in an essential way. The relevant QKE’s are nonlinear and strongly coupled internally via the asymmetry term, and therefore no useful analytical approximation exists for the problem222 Some confusion in the field was caused by a recent analytical treatment , whose results contradicted earlier numerical results . Disagreement was eventually clarified in favour of numerical work in .. Numerical solution is also greatly complicated by the various different physical scales, and while the phenomenon of asymmetry growth is fairly well established by now, the determinacy of the final sign of the asymmetry is still not well understood. The ambiguity was first observed in , and in  it was shown that this “chaoticity” occurs in certain well defined region of mixing parameters. References  employed a numerical solution for the momentum averaged approximation of the QKE’s however, and their results were challenged by ref. , who claimed that the final asymmetry is always fully determined by the initial conditions. While the analysis of ref.  is now discredited , it still remains true that no numerically reliable momentum dependent analysis has so far shown the existence of the chaotic region.333 Although ref.  found support for chaoticity, their conclusion is not firm because of the reported loss of the numerical accuracy of their methods when approaching the potentially chaotic region. It is important to settle the issue because, as mentioned above, large $`L_{\nu _e}`$ could significantly alter the helium-4 abundance. Indeed, if the sign of $`L_{\nu _e}`$ were found to be chaotic, then $`Y_{He}`$ could not be reliably determined from BBN .

In this paper we study the dynamics of active-sterile neutrino oscillations in the early universe using full momentum-dependent quantum-kinetic equations (in homogeneous space). Our main technical improvement over the previous works is the introduction of a novel dynamical discretization method for the momentum variable which enables us to model accurately the density matrices over the entire momentum range, including the extremely pronounced and narrow structures close to resonant momenta. On physical context, our results confirm the existence of a ”chaotic” region of mixing parameters, for which the final sign of the asymmetry is not deterministic. It also confirms the expectation , that the size of the chaotic region is somewhat smaller than indicated by the momentum averaged code . These results can be seen as giving more strength on the cosmological constraints on neutrino mixing parameters.

The paper is organized as follows. In section 2 we set up the quantum kinetic equations for the problem. In section 3 we introduce the novel dynamically adjusted discretization of the momentum variable and in section 4 we set up our kinetic equations with this parametrization. The method we develop allows us to have enough resolution to solve the problem with relatively small number of momentum bins (we use at most 400 bins). In section 5 we display our numerical results for the key variables driving the oscillation and discuss the numerical stability of our solutions. Finally, in section 6 we present our results on the asymmetry growth and oscillations as well as the interpretation of the physical consequences, and the section 7 contains a summary and outlook.

## 2 Quantum Kinetic Equations

In the early Universe the oscillating neutrinos experience frequent scatterings which interrupt the coherent evolution of the state and introduce collisional mixing between different momentum states. The mathematical formalism which can incorporate all these features has been developed elsewhere , and it takes the form of the quantum kinetic equations for the reduced density matrices for the neutrino and antineutrino ensembles. We parametrize the density matrices by the Bloch vector presentation :

$$\rho _\nu \frac{1}{2}f_0(P_0+𝐏\sigma ),\rho _{\overline{\nu }}\frac{1}{2}f_0(\overline{P}_0+\overline{𝐏}\sigma ),$$
(4)
where $`\sigma `$ are the Pauli spin matrices and $`f_0=(1+\mathrm{exp}(p/T))^1`$ is the Fermi-Dirac distribution function without chemical potential. Then to the the lowest order in the interaction , the QKE’s for the density matrix take the form

$`d_tP_0`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{f_0}}\left[f_{eq}\rho _{\alpha \alpha }\right]`$
$`d_tP_x`$ $`=`$ $`V_zP_yDP_x`$
$`d_tP_y`$ $`=`$ $`V_zP_xV_xP_zDP_y`$
$`d_tP_z`$ $`=`$ $`V_xP_y+{\displaystyle \frac{\mathrm{\Gamma }}{f_0}}\left[f_{eq}\rho _{\alpha \alpha }\right],`$ (5)
where $`d_t_tHp_p`$ and $`H`$ is the Hubble expansion factor. The rotation vector $`𝐕`$ has the following components ($`V_y0`$): in $`x`$-direction one has only the vacuum contribution

$$V_x=\frac{\delta m^2}{2p}\mathrm{sin}2\theta _0,$$
(6)
whereas in the $`z`$-direction also the matter induced effective potential contributes:

$$V_z=V_0+V_1+V_L,$$
(7)
where

$`V_0`$ $`=`$ $`{\displaystyle \frac{\delta m^2}{2p}}\mathrm{cos}2\theta _0,`$ (8)
$`V_1^e`$ $`=`$ $`14\sqrt{2}{\displaystyle \frac{\zeta (4)}{\zeta (3)}}{\displaystyle \frac{G_F}{M_W^2}}N_\gamma pT\left(1+\frac{1}{4}\mathrm{cos}^2\theta _W\left[n_{\nu _e}+n_{\overline{\nu }_e}\right]\right)`$ (9)
$`V_1^{\mu ,\tau }`$ $`=`$ $`{\displaystyle \frac{7\sqrt{2}}{2}}{\displaystyle \frac{\zeta (4)}{\zeta (3)}}{\displaystyle \frac{G_F}{M_Z^2}}N_\gamma pT\left[n_{\nu _{\mu ,\tau }}+n_{\overline{\nu }_{\mu ,\tau }}\right]`$ (10)
$`V_L`$ $`=`$ $`\sqrt{2}G_FN_\gamma L^{(\alpha )}.`$ (11)
Here $`N_\gamma `$ is the photon equilibrium number density, $`n_{\nu _{\alpha ,(\overline{\alpha })}}`$ is the normalized to unity (in equilibrium) actual neutrino (antineutrino) number density, and the effective neutrino asymmetries $`L^{(\alpha )}`$ are given by

$`L^{(e)}`$ $`=`$ $`\left(\frac{1}{2}+2\mathrm{sin}^2\theta _W\right)L_e+\left(\frac{1}{2}2\mathrm{sin}^2\theta _W\right)L_p\frac{1}{2}L_n+2L_{\nu _e}+L_{\nu _\mu }+L_{\nu _\tau }`$ (12)
$`L^{(\mu )}`$ $`=`$ $`L^{(e)}L_eL_{\nu _e}+L_{\nu _\mu }`$ (13)
$`L^{(\tau )}`$ $`=`$ $`L^{(e)}L_eL_{\nu _e}+L_{\nu _\tau }.`$ (14)
where $`L_f(n_fn_{\overline{f}})N_f/N_\gamma `$.

The repopulation term $`\mathrm{\Gamma }(f_{eq}\rho _{\alpha \alpha })`$ in the diagonal part of (5) is written in the relaxation time approximation . The distribution $`f_{eq}`$ in the repopulation term is the usual equilibrium Fermi-Dirac distribution

$$f_{eq}=\frac{1}{1+e^{\frac{p}{T}\frac{\mu }{T}}}.$$
(15)
and the reaction rates are

$$\mathrm{\Gamma }=C_\alpha G_F^2pT^4$$
(16)
with $`C_e1.27`$ and $`C_{\mu ,\tau }0.92`$. The damping terms $`DP_i`$ correspond to the decohering scatterings, which tend to destroy the coherence of the oscillation. In density matrix language such terms appear as relaxation terms suppressing the off-diagonals of the density matrix, or equivalently the $`P_{x,y}`$-components of the Bloch vectors. The magnitude of the decohering terms is just half of the corresponding scattering rate

$$D=\frac{1}{2}\mathrm{\Gamma }.$$
(17)
The equation of motion for anti-neutrinos can be found by substituting $`L^{(\alpha )}=L^{(\alpha )}`$ and $`\mu _\nu =\mu _\nu `$ to the above equations (the latter condition is true when neutrinos are in chemical equilibrium).

The repopulation term used above is only an approximation for the correct elastic collision integral, and unfortunately it breaks the lepton number conservation. This is not a serious problem however, and it can be circumvented by introducing an explicit lepton number conserving evolution equation for $`L^{(\alpha )}`$. The appropriate equation, introduced in , is

$$d_tL^{(\alpha )}=\frac{1}{8\zeta (3)}\underset{0}{\overset{\mathrm{}}{}}𝑑pp^2(V_xP_y\overline{V}_x\overline{P}_y).$$
(18)
Equation (18) can be derived from the original equations in the approximation where the $`L`$-violation due to the the approximate repopulation term is neglected. So, when using (18) to evolve $`L`$, it is guaranteed that the repopulation term does not directly affect the lepton asymmetry generation. The accuracy of the method can be monitored by comparing the value of $`L^{(\alpha )}`$ derived from (18) and from directly integrating the neutrino distribution functions. In practice the agreement was always to be good.

Neutrino and antineutrino sectors are extremely strongly coupled through the asymmetry $`L^{(\alpha )}`$ in equations (5). In order to compute $`L^{(\alpha )}`$ accurately enough despite the numerical round-off errors, it is then necessary to avoid computing $`L^{(\alpha )}`$ through differences of large quantities, such as naively appear in the equation (18). To this end we will define the “small” and “large” linear combinations of dynamical variables in particle and antiparticle sector

$$P_i^\pm =P_i\pm \overline{P}_i.$$
(19)
For convenience we also separate active and sterile sectors by defining

$`P_a^\pm `$ $`=P_0^\pm +P_z^\pm `$ $`=2\rho _{\alpha \alpha }^\pm /f_0`$ (20)
$`P_s^\pm `$ $`=P_0^\pm P_z^\pm `$ $`=2\rho _{ss}^\pm /f_0.`$ (21)
In terms of variables (20-21) the equations take the form

$`d_tP_a^\pm `$ $`=`$ $`V_xP_y^\pm +\mathrm{\Gamma }\left[2f_{eq}^\pm /f_0P_a^\pm \right],`$ (22)
$`d_tP_s^\pm `$ $`=`$ $`V_xP_y^\pm ,`$ (23)
$`d_tP_x^\pm `$ $`=`$ $`(V_0+V_1)P_y^\pm V_LP_y^{}DP_x^\pm ,`$ (24)
$`d_tP_y^\pm `$ $`=`$ $`(V_0+V_1)P_x^\pm +V_LP_x^{}{\displaystyle \frac{1}{2}}V_x(P_a^\pm P_s^\pm )DP_y^\pm ,`$ (25)
where we assumed $`\mathrm{\Gamma }=\overline{\mathrm{\Gamma }}`$ and defined

$$f_{eq}^\pm =f_{eq}(p,\mu )\pm f_{eq}(p,\mu )$$
(26)
Because $`\overline{V}_x=V_x`$, the equation for the asymmetry now reads simply

$$d_tL^{(\alpha )}=\frac{1}{8\zeta (3)}\underset{0}{\overset{\mathrm{}}{}}𝑑pp^2V_xP_y^{}.$$
(27)
Finally, due to the particular form of the differential operator $`d_t=_tHp_p`$, it is possible to reduce our set of partial differential equations to ordinary differential equations by changing variables $`pxp/T`$, and $`tT`$, and correspondingly using the time-temperature relation

$$\frac{dT}{dt}=HT.$$
(28)
The differential operator then becomes simply

$$d_t\frac{1}{HT}_T.$$
(29)
Let us stress that our introducing the small and large linear combinations is an absolutely necessary step for obtaining reliable solutions to QKE’s (22-25,27). This can best be appreciated from the figure (3) below, where we show the energy spectrum in neutrino and antineutrino sectors for a representative set of oscillation parameters; the asymmetry corresponds to the integral over the difference (not even visible!) of the neutrino and antineutrino distributions displayed.

## 3 Discretization of the momentum variable

While necessary, introducing parametrization (19) is unfortunately not sufficient to tackle the momentum dependent QKE’s. Additional difficulties arise due to structures in the momentum direction (in variable $`x`$) and in particular the very narrow ones introduced by the resonances. In this paragraph we explain step by step how we discretize the momentum variable such that these difficulties can be overcome.

First, the momentum variable ranges over the semi-infinite range $`x[0,\mathrm{}]`$, but in practice we will introduce cut-offs to both ends of spectrum. First, the ultra-relativistic neutrino approximation breaks down when $`x0`$, formally appearing as a singularity in at $`x=0`$ in our equations. Second, it is useful to cut out very large $`x`$, in the exponentially suppressed tail of the distribution. So, we restrict $`x`$ to a range $`x[x_{min},x_{max}]`$, where in practice it suffices to use444Note that even for $`x=10^4`$ still $`p1`$ keV throughout our computation, so that the ultrarelativistic approximation always remains valid. $`x_{min}=10^4`$ and $`x_{max}=100`$.

Second, the variable $`x`$ is not ideal for discretization, because most of the variation in the spectrum occurs at small $`x3`$. We therefore have mapped the variable $`x`$ to $`u(x)`$ belonging to range $`[0,1]`$ by a transformation

$$u(x)=\frac{x(1+ϵ_2)x_{\mathrm{ext}}ϵ_1}{x+x_{\mathrm{ext}}}.$$
(30)
where $`ϵ_1k_1(1+ϵ_2)`$ and $`ϵ_2(1+k_1)/(k_2k_1)`$ with $`k_1x_{\mathrm{min}}/x_{\mathrm{ext}}`$ and $`k_2x_{\mathrm{max}}/x_{\mathrm{ext}}`$. Apart from a small correction due to the cut-off parameters (30) is designed so that the extremum point of the thermal momentum distribution $`x_{ext}`$ gets mapped to $`u1/2`$. We will also need the inverse of this function, which is

$$x(u)=\frac{x_{\mathrm{ext}}(u+ϵ_1)}{1+ϵ_2u}.$$
(31)
While clearly a significant improvement for binning the momentum variable, (30) is not clever enough to solve the numerical problems arising from the the pronounced structures around the resonant momenta. Fortunately, the resonance positions $`x_{r_1}`$ and $`x_{r_2}`$ can be solved analytically from the conditions $`V_z=0`$ (neutrinos) and $`\overline{V}_z=0`$ (antineutrinos):

$`x_{r_1}`$ $`=`$ $`\sqrt{\phi ^2+\chi }\phi \overline{\nu }\mathrm{resonance}`$ (32)
$`x_{r_2}`$ $`=`$ $`\sqrt{\phi ^2+\chi }+\phi \nu \mathrm{resonance}`$ (33)
where in the case of $`\nu _\mu `$ and $`\nu _\tau `$ neutrinos

$`\phi `$ $`=`$ $`{\displaystyle \frac{\zeta (3)M_Z^2L^{(\alpha )}}{14\zeta (4)T^2}}6.60\times 10^8L^{(\alpha )}T_{\mathrm{MeV}}^2`$ (34)
$`\chi `$ $`=`$ $`{\displaystyle \frac{\zeta (3)M_Z^2|\delta m^2|\mathrm{cos}2\theta _0}{14\sqrt{2}\zeta (4)G_FN_\gamma T^3}}1.64\times 10^8|\delta m^2|_{\mathrm{eV}^2}\mathrm{cos}2\theta _0T_{\mathrm{MeV}}^6.`$ (35)
Using (30) one easily finds the corresponding resonant values in the $`u`$-variable.

The resonance widths can be estimated from the distance over which the matter mixing angle drops to a some fraction of its maximum value of unity. In this way one finds that the relative resonance width is proportional to the vacuum mixing angle:

$$\frac{\mathrm{\Delta }x}{x}\mathrm{tan}2\theta _0.$$
(36)
From (36) it is clear that for small mixing angles the structures near resonances can only be resolved by an extremely finely spaced momentum grid. For example555One is driven to such small values of vacuum mixing to find acceptable phenomenology for the asymmetry growth . with $`\mathrm{sin}^22\theta _0=10^6`$, equation (36) predicts structures in the scale $`\mathrm{\Delta }x/x10^3`$. Assuming then that at least 50 points are needed to describe the resonant region accurately, one sees that a linear binning of the $`x`$-variable with $`x_{max}=100`$ would require $`5\times 10^6`$ grid points! Solving such problem is clearly beyond capabilities of any computer we have today.

Our mapping (30) helps the situation by roughly a factor of hundred, but the remaining problem would still be way too large. We therefore need to find a second transformation $`uu(v)`$ that would redistribute the points such that more grid points are clustered around the resonances in the physical momentum variable. There are obviously many ways to implement such a transformation. We found it most convenient to construct it by using polynomials of the form $`a+b(vv_r)^3`$. Because we have two resonance points the mapping function is a little complicated:

$$u(v)=\alpha v+(a+b(vv_{r_1})^3)\theta (v_cv)+(c+d(vv_{r_2})^3)\theta (vv_c),$$
(37)
where the two separate mappings around resonant values $`v_{r_1}`$ and $`v_{r_2}`$ are joined smoothly at the geometric mean of the resonances $`v_c\frac{1}{2}(v_{r_1}+v_{r_2})`$. The linear term $`\alpha v`$ was added in order to have a control how steep is the distribution function of the grid points.

The mapping (37) has six free parameters $`a`$, $`b`$, $`c`$, $`d`$, $`v_{r_1}`$ and $`v_{r_2}`$, which will be determined from boundary conditions at $`v=0`$ and $`v=1`$ and continuity conditions at $`v=v_c`$. First, the constants $`a`$, $`b`$, $`c`$ and $`d`$ are fixed by the boundary conditions

$`u(0)`$ $`=`$ $`0,`$ (38)
$`u(v_{r_1})`$ $`=`$ $`u_{r_1},`$ (39)
$`u(v_{r_2})`$ $`=`$ $`u_{r_2},`$ (40)
$`u(1)`$ $`=`$ $`1,`$ (41)
which results in

$`a`$ $`=`$ $`u_{r_1}\alpha v_{r_1}`$ (42)
$`b`$ $`=`$ $`{\displaystyle \frac{u_{r_1}}{v_{r_1}^3}}{\displaystyle \frac{\alpha }{v_{r_1}^2}}`$ (43)
$`c`$ $`=`$ $`u_{r_2}\alpha v_{r_2}`$ (44)
$`d`$ $`=`$ $`{\displaystyle \frac{1u_{r_2}}{(1v_{r_2})^3}}{\displaystyle \frac{\alpha }{(1v_{r_2})^2}}.`$ (45)
Note that $`u_{r_i}`$ are not unknowns, but can be solved analytically from the resonance conditions (32-33) and the inverse mapping (31). Finally, the continuity conditions

$`\underset{vv_c+}{lim}`$ $`u(v)`$ $`=\underset{vv_c}{lim}u(v)`$ (46)
$`\underset{vv_c+}{lim}`$ $`_vu(v)`$ $`=\underset{vv_c}{lim}_vu(v).`$ (47)
provide implicit equations for the remaining two unknowns $`v_{r_i}`$:

$`a+{\displaystyle \frac{1}{8}}b(v_{r_2}v_{r_1})^3`$ $`=`$ $`c+{\displaystyle \frac{1}{8}}d(v_{r_1}v_{r_2})^3`$ (48)
$`b`$ $`=`$ $`d.`$ (49)
Unfortunately these are sixth order algebraic equations and can only be solved by iterative numerical algorithms.

The novel feature of the above change of variables is that it automatically follows the resonances, placing majority of the grid points dynamically in their vicinity. As an alternative to our approach, one could think of using library routines with automatic grids to take care of this distribution. However, the resonance positions evolve as a function of expansion scale of the universe and in particular as a function of $`L^{(\alpha )}`$. High point densities are thus needed at different momenta at different times, and this movement can be fast. None of the library routines we tried came close to being able to cope with these requirements and failed completely to solve the problem. With our parametrization on the other hand, the numerical program can be written for a fixed grid in the variable $`v`$, and the redistribution of points is provided by the dynamical equations themselves.

In Fig. (1) we show the transformations $`uv`$ when $`L^{(\alpha )}`$ is almost zero (dotted line), and when $`L^{(\alpha )}`$ has grown to a substantial value (solid line). In the former case one cannot distinguish the two resonances by eye, whereas in the second case the structures have become very distinct. Fig. (2) shows the density of points as a function of $`x`$ (defined as the inverse of separation of grid points in the $`x`$-variable).

## 4 Changes in the QKE’s due to the parametrization

The changes of variables (30) and (37) affect the partial derivatives appearing in Eqns. (22-25,27). These equations have been written assuming the evolution of $`\rho `$ and $`L`$ is along constant $`x`$ contours. However, we have now fixed our grid in the variable $`v`$, and hence need to find the evolution equations with $`T`$ and $`v`$ as the independent variables. These are easily found from the old ones by the chain rule partial differentiation

$`\left({\displaystyle \frac{\rho (T,x(T,v))}{T}}\right)_v`$ $`=`$ $`\left({\displaystyle \frac{\rho }{T}}\right)_x+\left({\displaystyle \frac{x}{T}}\right)_v\left({\displaystyle \frac{\rho }{x}}\right)_T`$ (50)
$`=`$ $`\left({\displaystyle \frac{\rho }{T}}\right)_x+\left({\displaystyle \frac{x}{T}}\right)_v\left({\displaystyle \frac{u}{x}}\right)_T\left({\displaystyle \frac{v}{u}}\right)_T\left({\displaystyle \frac{\rho }{v}}\right)_T`$
$`=`$ $`\left({\displaystyle \frac{\rho }{T}}\right)_x+\left({\displaystyle \frac{u}{T}}\right)_v\left({\displaystyle \frac{v}{u}}\right)_T\left({\displaystyle \frac{\rho }{v}}\right)_T.`$
In the last equality we used the fact that $`u=u(x)`$ only. The differentials $`(_T\rho )_x`$ are of course just given by the l.h.s. of equations (25), and the last coefficient $`(_v\rho )_T`$ is computed numerically as a part of the system of the partial differential equations. To be precise, we are using the central derivative formula

$$\left(\frac{\rho }{v}\right)_T=\frac{\rho (v+h)\rho (vh)}{2h}.$$
(51)
Of the remaining differentials $`(_vu)_T`$ is also easily solved from the parametrization Eq. (37:

$$\left(_vu\right)_T=\frac{1}{\alpha +3[b(vv_{r_1})^2\theta (v_cv)+d(vv_{r_2})^2\theta (vv_c)]}.$$
(52)
However, solving the remaining differential $`\left(_Tu\right)_v`$ is somewhat more complicated. Using equation (37) one can write

$`\left(_Tu\right)_v`$ $`=`$ $`\left[_Ta+_Tb(vv_{r_1})^33b(vv_{r_1})^2_Tv_{r_1}\right]\theta (v_cv)`$ (53)
$`+\left[_Tc+_Td(vv_{r_2})^33d(vv_{r_2})^2_Tv_{r_2}\right]\theta (vv_c),`$
where

$`_Ta`$ $`=`$ $`_Tu_{r_1}\alpha _Tv_{r_1}`$ (54)
$`_Tb`$ $`=`$ $`{\displaystyle \frac{1}{v_{r_1}^3}}\left(_Tu_{r_1}+\left(2\alpha {\displaystyle \frac{3u_{r_1}}{v_{r_1}}}\right)_Tv_{r_1}\right)`$ (55)
$`_Tc`$ $`=`$ $`_Tu_{r_2}\alpha _Tv_{r_2}`$ (56)
$`_Td`$ $`=`$ $`{\displaystyle \frac{1}{(1v_{r_2})^3}}\left(_Tu_{r_2}+\left(2\alpha 3{\displaystyle \frac{1u_{r_2}}{1v_{r_2}}}\right)_Tv_{r_2}\right).`$ (57)
Here the only unknowns are the temperature derivatives of the resonances in the $`v`$-variable $`_Tv_{r_i}`$. These can be solved from our original matching conditions (48) and (49), which can formally be written as

$`f_1(v_{r_1}(T),v_{r_2}(T),T)`$ $`=`$ $`0`$ (58)
$`f_2(v_{r_1}(T),v_{r_2}(T),T)`$ $`=`$ $`0.`$ (59)
Differentiating equations (58-59) with respect to $`T`$ results in a linear set of equations for $`_Tv_{r_i}`$, which after a little algebra can be brought to the form

$$_Tv_{r_i}=U_{ij}_Tu_{r_j},$$
(60)
where the matrix $`U_{ij}`$ depends both on $`u_{r_i}`$ and $`v_{r_i}`$. Note that despite the apparent complexity of our parametrization, we have been able to reduce everything except solving for the resonance parameters $`v_{r_i}`$ to a simple linear algebra. Solving $`v_{r_i}`$ from the matching equations at each time step would be very time consuming however. This is so in particular because we need to solve them with very high accuracy in order for not to lose the continuity of the matching at $`v=v_c`$ and to avoid inducing random numerical errors into the equations. Fortunately, the equation (60) itself offers the way out of this problem. Namely, we only need to solve $`v_{r_i}`$ from the matching conditions (48) and (49) once, in the beginning of the iteration, as their value can later be evolved by equations (60) as the solution proceeds. In this way even solving for the parametrization from $`v`$ to $`x`$ becomes part of the dynamical equations. The accumulation of numerical error can be easily traced by solving the matching conditions independently for some temperatures. In practice the errors do not accumulate at all and it turns out that using the differential equation is by far the superior method for solving $`v_{r_i}`$ in comparison to using the matching conditions at each temperature step.

## 5 Numerical Results

We have solved the evolution equations numerically for a range of representative parameter sets. In figure (3) we show the neutrino and antineutrino distribution functions in a relatively early stage in the evolution for $`\delta m^2=10\mathrm{eV}^2`$ and $`\mathrm{sin}^22\theta _0=10^{6.1}`$. The distributions show a very nice thermal shape, and no visible difference between particle and antiparticle sectors can be seen. This is so despite the fact that at the temperature at which the distributions were plotted, the momenta $`p1.8T`$ are currently resonant (indicated by the gray line in the plot); the amplitude of resonant structure is simply too small to be seen in the scale of the “large variables”.

The resonance does alter the asymmetry spectrum however. In figure (4) we show the distribution of the difference of neutrino and antineutrino densities $`\rho _{\alpha \alpha }\rho _{\overline{\alpha }\overline{\alpha }}`$ in the physical variable $`x`$. The prominent and extremely narrow resonance structure (one cannot separate the neutrino and antineutrino resonances here) is clearly visible. The same spectrum is shown in the integration variable $`v`$ in figure (5). Now the sharp kink-like resonant structure has broadened into a broad wave, which is easily represented by the linearly discretized function in variable $`v`$ (grid points are shown by the open circles). In this example we used just 400 grid points, but essentially identical results were found with only 200 points. Moreover, comparing the scales in figures (3-4), one sees that the scale of variation of $`L(x)`$ is ten orders of magnitude below unity. So, in order to follow the evolution of $`L^{(\alpha )}`$ accurately when using the original variables, one should be able to integrate over distributions with an accuracy much better than ten digits! This seems totally impossible, and clearly indicates that our use of variables $`P^\pm `$ is essential.

These examples prove the success of the key elements in our approach, namely that by separating the small and large variables, and by introducing the novel dynamical discretization of the momentum grid, we can model the evolution of the ensemble of oscillating neutrinos accurately. Moreover, we can do so using small enough lattices such that the code can be run even in a powerful modern workstation (although a single run may then take as long as a couple of days).

### 5.1 Regulating the sterile neutrino spectrum

As one lets the system evolve further the resonance $`x`$ moves to larger values while the integrated $`L^{(\alpha )}`$ decreases. This continues until one reaches roughly the point when the main bulk of the active neutrino distribution is getting through the resonance. At this point one expects that the total asymmetry begins to oscillate or grow, depending on the parameters. In our initial runs we did not get that far however, but the program crashed due to an accumulation of structures in sterile neutrino spectrum with momenta less than $`x_{res}T`$. How these structures arise is easy to understand. Indeed, the active sector distributions are kept smooth away from the resonances by diffusion induced by collisions with the background particles. Sterile states interact only through the effective interaction due to mixing however, and since this interaction is neglibly small away from the resonances, the diffusion does not work on sterile sector. As a result the moving resonance leaves the sterile neutrino momentum distribution quite oscillatory for $`x<x_{res}`$. While these developments in the sterile sector have little effect on variables in the active sector (because active and sterile states are essentially decoupled away from resonances), they still pose a numerical problem: as the resonance moves forward in $`x`$, the density of points reduces drastically in the region $`xx_{res}`$, until the point that the grid becomes too sparse to model the oscillatory structures in the unsmoothed sterile neutrino distribution. Eventually the amplitudes of these structures overflow causing numerical instability and the solution breaks down.

The above explanation of the instability also suggests an obvious cure for the problem. If one is not interested in studying the sterile neutrino asymmetry spectrum, all we need to do is to add a small regulatory interaction term to sterile neutrino sector, which will smooth their distribution away from resonances, and yet does not alter the evolution of the variables in the active sector. The latter requirement can be met by a suitable form of the interaction term and by making it small enough so that it does not change the physics in the resonance region. The challenge is to do this so that the regulatory term remains efficient enough to numerically stabilize the system. To this end we have added the following “sterile repopulation terms” onto equation (23):

$`R_{\nu _s}^+`$ $`=`$ $`r_s\mathrm{\Gamma }\left(n_{\nu _s}f_{eq}^+(x,\mu )\rho _{ss}^+\right)`$ (61)
$`R_{\nu _s}^{}`$ $`=`$ $`r_s\mathrm{\Gamma }\left(\rho _{aa}^{\mathrm{init}}\rho _{aa}^{}\rho _{ss}^{}\right).`$ (62)
Here $`\rho _{aa}^{\mathrm{init}}`$ is the initial density matrix in the active sector. We used the active sector interaction rate $`\mathrm{\Gamma }`$ above, and let the parameter $`r_s`$ control the size of the regulatory term. The regulators (61-62) differ from the repopulation terms in the active sector in (22) in that we have added normalization factors to make sure that the integrals over $`R_{\nu _s}^\pm `$ always vanish. (In this way we avoid spurious sterile state equilibration even with large control parameter $`r_s`$.)

We obviously need to show that there exists a range of values for which the results converge and become independent of $`r_s`$, while the regulator still continues to stabilize the solution. That we do reach this goal is evidenced in Fig. (6), where we show the evolution of the total asymmetry for $`r_s=0.1`$, 0.01, 0.001, $`10^4`$ and $`10^5`$, for $`\delta m^2=0.1\mathrm{eV}^2`$ and $`\mathrm{sin}^22\theta _0=10^{6.1}`$. Unsurprisingly, for large $`r_s`$ there is a significant effect, but for small $`r_s`$ the results do converge to the point that no visible difference can be seen between the two curves corresponding to two smallest $`r_s`$’s. Let us stress that the mixing parameters were here deliberately chosen such as to show a particularly large sensitivity for the regulator. For many other cases that we studied it would have been difficult to show any effect at all; for example the results displayed in figure (8) showed no appreciable dependence on the regulator up to $`r_s=10`$.

We have also studied the stability of our results against changes in various numerical error handling parameters in our code. The stepwise error tolerance for example can be relaxed by an order of magnitude from what we actually employed, before any visible deviation from the solution can be seen over the entire integration range. This is true also for our examples corresponding to the oscillatory solutions in or near the chaotic region. Of course, the numerical accuracy would be eventually compromised for parameters for which the asymmetry undergoes very large number of oscillations before settling to the growing curve, but as in , that does not affect our main conclusions. Also cut off parameters $`x_{min}`$ and $`x_{max}`$ and the tilt $`\alpha `$ which controls the density of points at the resonance can be varied within reasonable limits without any observable effects on the results. We therefore are confident that our observed pattern of asymmetry oscillations is indeed physical, and not of numerical origin.

### 5.2 Argument for chaoticity in $`L`$-growth

Having introduced a method that can handle the evolution of momentum dependent density matrices and in particular the asymmetry, and having proven its numerical stability and accuracy, we now pursue the question of the “chaotic” behaviour of the final sign of the integrated asymmetry. Note that while the use of word chaos is customary in this context, we strictly speaking mean only sensitivity of the final sign of the asymmetry on initial seed asymmetry and on oscillation parameters666 In fact the system does exhibit true chaoticity at a certain level, but the related information loss is very small in practice ..

In the context of our earlier work using the momentum averaged equations , we explained what was the physical origin for the sign sensitivity. Quite simply, changes in the boundary conditions (initial baryon asymmetry and the oscillation parameters) change the direction of the motion of the system, and the topography of the attractor solutions in the phase space at the onset of resonance. Roughly speaking, in the momentum averaged case one can draw an analogy between the system and a somewhat sticky ball initially rolling down a single valley floor that later branches to two, with a low maximum in between, at the resonance point. Changing baryon asymmetry would then be analogous to changing the speed and direction of the ball at the branching point and changing oscillation parameters to that of changing the shape of the valleys themselves. Sign sensitivity arises for those initial conditions and valley forms for which the ball makes many oscillations between the two degenerate valleys before stickiness (damping) eventually wins and makes it to settle into either one of them.

In the momentum dependent case these structures become smeared out, and to draw a mechanical analogy, one should rather think, instead of a ball, of a string of varying size beads (largest ones in the middle, corresponding to the peak in the thermal distribution) coupled by a very elastic cord, wiggling down the valleys described above. While this is obviously a much more complicated system with qualitatively new features, one would still expect to see oscillations in the weighted average position of beads at least when the largest ones roll past the branching point (in $`L`$ when the maximum of the distribution crosses the resonance). Moreover, if the average position of the string ($`L`$) makes many oscillations over the central hill before stickiness (damping) wins, one again should expect that the final valley that the string settles in (sign of $`L`$) strongly depends on initial and boundary conditions. While this analog is certainly bears only a crude resemblance to the true system, it should help convince the reader of the fact that there is nothing mysterious about the sign sensitivity in the $`L`$-evolution, but that it is rather something to be expected.

Obviously, all we need to do to prove the existence of chaoticity, is to find examples of almost identical sets of oscillation parameters for which the final asymmetry does show significantly different oscillation pattern, and leads to a different final sign. We show such a case in figure (7). All three curves in the figure correspond to $`\delta m^2=10\mathrm{eV}^2`$ but have slightly different vacuum mixing angles. The two curves which behave almost alike, have mixing angles $`\mathrm{sin}^22\theta _0=10^{6.1}`$ (solid curve) and $`\mathrm{sin}^22\theta _0=10^{6.099}`$ (dotted curve) respectively. The dash-dotted curve with a very different behaviour corresponds to $`\mathrm{sin}^22\theta _0=10^{6.09}`$. So, if such mixing were ever observed and the mixing angle was measured to a precision better than a tenth of a per cent, the final sign of the lepton asymmetry and its effect on nucleosynthesis could be predicted from our computation. However, if the mixing angle was measured with less than about one per cent accuracy, we would not be able to say what the final asymmetry will be, and therefore what is the precise BBN prediction for helium abundance.

The sign sensitivity does not occur for all oscillation parameters, however. In the figure (8) we show results from “stable” region of parameters (same mass difference as above, but much smaller mixing angle). In this case, while visible differences in the evolution become visible when mixing angle is changed by a factor or three, the asymmetry never becomes oscillatory, and the final sign is determined by the initial sign of the asymmetry.

Our first example was chosen from close to the edge of the chaotic region, and just as making $`\mathrm{sin}^22\theta _0`$ smaller made system more stable, increasing it has the opposite effect. The number of oscillations and also degree of sensitivity of the final sign to the oscillation parameters increases very rapidly with $`\mathrm{sin}^22\theta _0`$. Nevertheless, the degree of sensitivity to parameters is weaker and thus the chaotic region smaller than what was found in the momentum averaged treatment . Quantitatively our new results agree roughly with the boundaries of the chaotic region inferred in , but since the numerical work is, even with the improvements we have made, very time consuming, we will not attempt a precise mapping of the chaotic region here.

Let us finally to try to understand a little more deeply what precisely drives the oscillations of the total asymmetry in the momentum dependent case. A clue can be found from Figure (9) which shows the variable $`P_y^{}`$ in the region where the asymmetry oscillates significantly. (At earlier times, when the total $`L`$ does not oscillate, the structure of the $`P_y^{}`$-spectrum is similar to that of $`\rho _{\alpha \alpha }^{}`$ shown in Fig. (5).) The two new peaks seen outside the resonance grow slowly in time, while their sign oscillates rapidly. Because these structures are asymmetric around the resonance (in the physical variable $`x`$), the integral over $`P_y^\pm (x)`$ is nonzero and oscillates. As a result also $`L`$, which is solved from (27), becomes oscillatory.

What causes these structures, is a delicate interplay between the coherence length and the matter oscillation length. The former is simply given by the inverse of the damping rate $`\mathrm{}_{\mathrm{coh}}=1/D`$ and the latter is

$$\mathrm{}_\mathrm{m}=\frac{l_{\mathrm{vac}}}{\sqrt{12\chi \mathrm{cos}2\theta _0+\chi ^2}}$$
(63)
where $`l_{\mathrm{vac}}4\pi p/|\delta m^2|`$ and $`\chi 2p|V|/|\delta m^2|`$ where $`|V|`$ is the magnitude of the matter contribution to the $`V_z`$ in equation (7). Near the resonance the matter oscillation length is very long, and in particular at high temperatures it exceeds the coherence length by a large margin. In such a situation the resonance is strongly overdamped, and despite the maximal mixing angle oscillations do not have time to develop because the state is continuously projected to the active direction. This in part explains why the resonant features have so modest amplitudes at early times (cf. Fig. (4)).

In the upper panel of the figure (9) we have plotted the coherence length and $`\mathrm{}_\mathrm{m}/4`$, which corresponds to the distance over which $`P_y^\pm `$ first rises to maximum due to oscillation. It is clear that the new structures correspond to the case when the matter oscillations first time break the overdamping situation. Further away from the resonance the damping no more can block the oscillations, but since the matter mixing angle decreases very fast away from the resonance, the amplitude of the oscillation features dies off quickly. As the temperature decreases, the new structures move away from the resonance, and they also become broader (one can show that their separation scales like $`\delta xT^3`$). This is why their amplitude slowly increases in time (with decreasing $`T`$) causing larger and larger amplitude oscillations in the total asymmetry $`L`$. The lesson to bring home from these observations is that the dynamics of $`L`$-growth and oscillations is very dependent on the detailed fine structure of the various components of the density matrix. Such fine details on the hand can only be studied in a numerical approach using both separated small and large variables ($`P_i^\pm `$) and dynamically adjusting momentum grids.

## 6 Conclusions

Active-sterile mixing, while constrained by the present solar and atmospheric data, are nevertheless required if one wishes to incorporate the LSND-neutrino anomaly into the neutrino models. Moreover, the active-sterile neutrino mixing would have a rich phenomenology in the early universe, where it provides a unique theoretical challenge in the form of a system for which the quantum effects play an essential role in the kinetic equations.

Here we have studied the phenomenon of neutrino asymmetry growth in the early universe as a result of active sterile neutrino oscillations. Our main new contribution is the introduction of a novel method of discretizing the momentum variable such that the sharply pronounced structures near the resonances, which are here shown to drive the oscillation phenomena, can be treated numerically accurately. The only concession to rigor we made along the way to solution of the problem, was adding collision terms to regulate the sterile neutrino spectrum away from the resonances. However, we demonstrated that our regulators had no effect on the active neutrino evolution, and hence for the results presented here. Nevertheless, our treatment is obviously not adequate if the precise form of the sterile neutrino spectrum is important. This might be the case for the models where sterile neutrinos are invoked to provide a non-thermal component for the dark matter .

We have demonstrated that even tiny changes in the oscillation parameters may drastically alter the oscillation pattern, and even change the sign of the final asymmetry. This behaviour does not occur for all oscillation parameters, but instead the dependence of the sign of $`L`$ on oscillation parameters is weak in the “stable” region (in general small $`\mathrm{sin}^22\theta _0`$ and large negative $`\delta m^2`$) and strong in the chaotic region (in general large $`\mathrm{sin}^22\theta _0`$ and small negative $`\delta m^2`$). These results are in qualitative agreement with our earlier findings, which were based on a momentum averaged treatment . Obviously, if active-sterile mixing were to be observed with parameters residing in the chaotic region, it would be, due to inevitable errors in the experimentally measured parameters, impossible to reliably determine the sign of the neutrino asymmetry created by that mixing in the early universe. Moreover, since large neutrino asymmetries (and electron neutrino asymmetry in particular) directly affect, in a $`sign(L)`$-dependent way, the weak interaction rates that determine proton-to-neutron ratio in the early universe, this indeterminacy would undermine our ability to accurately compute the BBN prediction for light element abundances for chaotic oscillation parameters.

These conclusions hold in the spatially homogeneous calculation. However, the sensitivity of $`L`$-growth on initial conditions also lends to a speculation that in reality the initially inhomogeneous seed asymmetry might give rise to an inhomogeneous texture of domains of very large lepton asymmetries with oscillating sign. This phenomenon has been shown to occur in a one-dimensional, momentum averaged model , and could plausibly occur in a realistic three dimensional world. The implications of such a scenario for SBBN would obviously be very different, including the possibility of very efficient equilibration of sterile neutrinos via MSW-effect within the domain boundaries . If this was the case, then the asymmetry growth mechanism would lead to even stronger bounds on mixing than what is displayed in equations (2). It is beyond the scope of this work (and the reach of the present computers) to study this phenomenon quantitatively in the momentum dependent case, however.

Let us finally comment on the effects of large homogeneous lepton asymmetries for cosmic microwave background radiation (CMBR). The possibility of measuring $`L`$ by the Planck satellite data has been considered for example in  and in . However, the only effect of $`L`$ on CMBR comes through the associated fluctuations in the energy density, and correspondingly only the total asymmetry has relevance. In the oscillation scenarios, such as discussed in this paper, the total asymmetry (here the sum of the sterile and active sector asymmetries) is conserved however, and hence the oscillation-induced asymmetries, in contrast to the situation with SBBN, would have no direct effect on CMBR.

###### Acknowledgments.

We thank Kari Enqvist for useful conversations and a collaboration at earlier stages of this project. We also thank Steen Hannestad for discussions on the effects of $`L`$ on CMBR. AS wishes to thank Nordita for hospitality during several visits in the course of completing this project.


format: text

text:
# A superburst from GX 3+1

## 1 Introduction

Recently, six long X-ray flares lasting several hours have been identified in five low-mass X-ray binaries (Cornelisse et al. 2000, 2001, Strohmayer & Brown 2001; Wijnands 2001; Kuulkers et al. 2001). For one source, 4U 1636$``$53, two such events were reported which occurred about 4.7 years after each other (Wijnands 2001). So far, these flares have only been seen in X-ray bursters with persistent pre-flare luminosities of $``$0.1–0.3 times the Eddington luminosity (Wijnands 2001; Kuulkers et al. 2001).

The long X-ray flares share many of the characteristics of type I X-ray bursts111Type I X-ray bursts have light curves with a rise which is faster than the exponential-like decay; their emission is well described by black-body radiation with temperatures, $`kT`$, around 2 keV and apparent black-body radii around 10 km; they show X-ray spectral softening during the decay. They have durations of seconds to minutes. For a review, see Lewin et al. (1993). and are, therefore, attributed to thermonuclear runaway events on a neutron star (e.g. Cornelisse et al. 2000). The differences with type I bursts are their long duration (exponential decay times of a few hours), their large fluences (about 1042 erg), and their rarity. Because of the large fluences, the flares are referred to as ‘superbursts’. A likely fuel for the superbursts is carbon, left over from stable and unstable hydrogen and/or helium burning (Cumming & Bildsten 2001; Strohmayer & Brown 2001). Unstable electron capture by protons with subsequent capture of the resulting neutrons by heavy nuclei is another conceivable option (Kuulkers et al. 2001).

I here report on a superburst seen with the Rossi X-ray Timing Explorer (RXTE) All Sky Monitor (ASM) from GX 3+1 in June 1998. For a preliminary announcement of this event see Kuulkers (2001). The overall X-ray intensity of the GX 3+1 varies slowly on time scales of months to years by a factor of about 2 (e.g. Makishima et al. 1983; see Fig. 1a). Type I X-ray bursts in GX 3+1 were first discovered by Hakucho (Makishima et al. 1983). A type I X-ray burst with radius expansion (due to the burst luminosity reaching the Eddington limit) was observed by the Proportional Counter Array (PCA) onboard RXTE, enabling one to estimate the distance to the source to be in the range 4–6 kpc (Kuulkers & van der Klis 2000).

## 2 Observations and Analysis

The ASM (Levine et al. 1996) is one of the three instruments onboard RXTE. It consists of three Scanning Shadow Cameras (SSCs) mounted on a rotating drive such that the center of the field of view of one camera is perpendicular to that of the other two. The cameras are held stationary for 90-s intervals, called ‘dwells’, during which data are accumulated. The mount rotates the assembly through a 6$`\mathrm{°}`$ angle between dwells until a full rewind is necessary. In this manner about 80% of the sky is observed every 90-min orbit around the earth. Each dwell provides intensities in the 1.5–12 keV band for all known sources in the field of view of each camera, and results that meet a set of reliability criteria are saved in electronic tables that are available over the World Wide Web.

The ASM data are also available in three sub-bands: 1.5–3 keV, 3–5 keV, and 5–12 keV. The boundaries are, however, crude approximations. The actual channel boundaries change somewhat with time, anode, software changes, etcetera. Moreover, the amount of these changes differ between the three SSCs (A. Levine 2001, priv. comm.). A pre-flight instrument effective area matrix for a source at the center of the field of view of an SSC is available, providing the detector response between 1 and 20 keV on an energy grid with a resolution of 0.1 keV222Note that this matrix is also used by PIMMS (http://legacy.gsfc.nasa.gov/docs/software/tools/ pimms.html).. It was used for estimating data system properties and rough sensitivities during the ASM development stage, and was certainly not designed for doing X-ray spectral fits to actual observations in three energy bands. There were small changes in the instrument design after the matrix numbers were computed, but the numbers were never updated (A. Levine 2001, priv. comm.). However, the matrix can be still useful for rough calculations, as done here. Note that a simple spectral model such as a black-body with only two free parameters may still be constrained by data in three energy bands.

The Crab (pulsar and nebula) is one of the main bright and steady X-ray calibration sources. Its observed ASM count rate is about 75 cts s-1 SSC-1 (1.5–12 keV, when at the center of an SSC field of view; see e.g. Levine et al. 1996). The pre-flight instrument effective area matrix, however, predicts a Crab count rate of about 88 cts s-1 SSC-1 in the same energy band. (The Crab spectrum in the energy band of interest can be described by an absorbed power-law; for the parameter values see e.g. Schattenburg & Canizares (1986), and references therein.) For doing X-ray spectral analysis one, therefore, has to renormalize the source count rates. The observed source count rates in a certain energy band were multiplied by the ratio of the predicted Crab count rate and the average observed Crab count rate in the same energy band. The latter were determined from the individual dwells during an interval centered around the observations, i.e MJD 50900–51050 (1998 March 28 to August 25). The resulting source count rates were then input for my spectral fitting. The rms in the average Crab count rates was included in the uncertainties of the source count rates (4.5%, 13%, 5% and 5%, in the 1.5–12 keV, 1.5–3 keV, 3–5 keV and 5–12 keV bands, respectively). Note that the errors in the individual dwell measurements contain already, apart from counting statistics, a 3% systematic error; according to the ASM team this will likely be an underestimate of the error in many cases. The method outlayed above is in principle only applicable if the source spectrum is comparable to that of the Crab. However, since I am only interested in rough estimates my procedure is adequate.

In the spectral analysis the interstellar absorption column, NH, was fixed at 1.7$`\times `$1021 cm-2 (Christian & Swank 1997). Whenever the goodness of fits in the spectral or temporal analyses, expressed as reduced $`\chi ^2`$ ($`\chi _{\mathrm{red}}^2`$) for certain degrees of freedom (dof), exceeded 1, the error bars in the data points were arbitrarily increased so that $`\chi _{\mathrm{red}}^2=1`$. I then determined uncertainties in the fitted parameters using $`\mathrm{\Delta }\chi ^2=1`$.

## 3 Results

By visually inspecting the ASM 1.5–12 keV database of GX 3+1 a single flare was found on 1998 June 9 (MJD 50973). It had a maximum count rate of about 82 cts s-1 SSC-1, which was about 5 times brighter than the count rate seen before and after the event (Fig. 1b). It started between UT June 8 22:53:30 and June 9 00:50:40. The decay is exponential-like with a decay time of $``$1.6 hr (1.5–12 keV). About 4.4 hours after the peak of the flare the count rate (27$`\pm `$2 cts s-1 SSC-1) was still above the pre-flare persistent level ($``$17 cts s-1 SSC-1), suggesting the event lasted even longer. From the end of the last pre-flare measurement and the first measurement at the persistent count rate after the flare (June 9 15:05:04) I infer an upper limit on the duration of the flare of 16.2 hours.

In Table 1 the fitted exponential decay times of the long flare in the total energy band and the three sub-bands are listed. They were determined using the flare data plus the persistent data up to a day after the start of the flare. The decay time in the 5–12 keV band is significantly shorter than that in the 3–5 keV band, indicating spectral softening during the decay. This was verified by determining hardness ratios, using the following procedure. I first averaged the count rates in the three energy bands of data within 30 min bins. Then the ratio of the count rates in a hard to low energy band was computed (there are five possible independent combinations with the three available energy bands). I then investigated the flare hardness values during the decay and compared them with those of the persistent emission. It turned out that ratio of the count rates in the 5–12 keV and 1.5–3 keV energy bands gave the best indication for spectral cooling during the day, see Fig. 1c. At flare maximum the source is harder with respect to just before the flare; during the decay of the flare the source softens.

To estimate the observed maximum flux and the fluence of the flare I performed X-ray spectral fits to the ASM data (see Sect. 2). The emission of normal type I bursts can be well described by emission from a black-body (Swank et al. 1977; Lewin et al. 1993). This also holds for the superbursts (e.g. Cornelisse et al. 2000). The flare was divided into three intervals which are indicated by dotted lines in Fig. 1b. The average pre-flare persistent emission between MJD 50971.5 and 50973 was subtracted from the total flare emission per interval. A black-body emission model subjected to interstellar absorption was used; The results of the fits are given in Table 1. Since the peak luminosity is well below the Eddington luminosity (see Sect. 4) the radius of the emission region is expected to be constant during the decay. I, therefore, also fitted the spectra from the three intervals together, coupling the apparent black-body radius, $`R_{\mathrm{bb}}`$, to one value. The results of these fits are also given in Table 1; I get $`R_{\mathrm{bb}}`$=5.4$`\pm `$1.5 km. The values of $`R_{\mathrm{bb}}`$ obtained in the above spectral fits are more or less comparable to those found during the decay of the radius expansion burst seen from GX 3+1 (Kuulkers & van der Klis 2000: $``$7.6 km at 5 kpc). Similar results are derived for e.g. KS 1731$``$260 (Kuulkers et al. 2001). I, therefore, also performed spectral fits with $`R_{\mathrm{bb}}`$ fixed at 7.6 km, of which the results are also shown in Table 1. Note that some of the values for $`\chi _{\mathrm{red}}^2`$ are unrealistically low; this is most likely due to an overestimate of the errors in the data points by including the rms in the average Crab values. The spectral fits, however, indicate that the black-body temperature, $`kT_{\mathrm{bb}}`$, is in the range 1–2 keV, and decreases, as expected, during the decay of the flare.

To get the values for the observed maximum flux and fluence, the individual dwell net-flare (= total flare minus pre-flare) count rates were converted into fluxes, using the count rates and average fluxes in the three intervals as obtained in the spectral fits with $`R_{\mathrm{bb}}`$ left free. The decay in net-flare flux was fitted with an exponential and assumed it approached zero at infinity. The fluence was calculated by integrating the exponential from the start of the flare up to infinity. I note that this may underestimate the fluence, since the flare probably started somewhat earlier. Extrapolating the fitted exponential backwards in time and assuming the maximum flux can not exceed the Eddington limit (see Sect. 4) an estimate of the upper limit on the fluence, $`E_{\mathrm{b},\mathrm{max}}`$, may be obtained. The results are displayed in Table 1.

An estimate of the recurrence time of such flares is an important constraint on the physics of the flare. In my case this is very uncertain, given only one observed flare in 5.8 yrs. One could naively deduce that flares happen every 5.8 yrs, but this obviously assumes they occur periodically and that the source was continuously observed in that interval. There are, however, about 15600 dwells of 90 s in that timespan, implying an exposure time of only 16 days; they are also not equally spaced in time and the number of dwells per day varies wildly up to about 140 per day. Another flare could thus have happened during data gaps of longer than the flare duration, whereas one could rule out another flare when there are more than one dwell data point with the source at its persistent count rate in a time span shorter than the flare duration. I used the observed flare to estimate the total effective exposure time during which only one flare occurred, which may be used as a rough indicator on the recurrence time. The flare profile was characterized by an exponential in the 1.5–12 keV band, as determined by the fit to the decay, with a start point corresponding to the observed peak of 82 cts s-1 SSC-1. A window of 0.2 day was used (approximately the lower limit on the duration of the flare), which was moved throughout all dwell data points. I determined whether or not the flare could fit within the window by determining the $`\chi _{\mathrm{red}}^2`$ between the expected profile and the observed data points (note that I varied the expected persistent count rates and logged the smallest values of $`\chi _{\mathrm{red}}^2`$). Whenever $`\chi _{\mathrm{red}}^2`$ exceeded 2.5 (about the same value as for the exponential fit to the decay light curve), I concluded that within the window no such flare could have occurred; if it was less than 2.5 a flare could have occurred and only the individual dwell durations were taken into account. In this way I estimate a total effective exposure time of about 490 days. Note that this estimate is dependent on the used window length. If I use a window of 0.1 or 0.3 day I derive $``$320 days and $``$600 days, respectively.

## 4 Discussion

I found an X-ray flare from GX 3+1 in RXTE/ASM data spanning about 6 yrs; the flare had a decay time of 1.6 hr and a duration of longer than 4.4 hr, but shorter than 16.2 hr. During the exponential-like decay the flare spectrum softened. I conclude that the flare has its origin in unstable thermonuclear burning. The total fluence of the event is between about 5$`\times `$1041 erg and 2$`\times `$1042 erg for a distance of 5 kpc. The maximum net-burst flux reached during a normal type I X-ray burst seen by Kuulkers & van der Klis (2000) was 6.9$`\times `$10-8 erg cm-2 s-1. Since this burst was a radius-expansion event the corresponding luminosity at the neutron star surface reached the Eddington limit. The flare, therefore, had an observed maximum of about 0.5 times the Eddington value (but note that the actual peak of the flare was missed). The persistent ASM count rate before and after the flare was similar to that observed around the radius expansion burst, indicating the persistent flux near the X-ray flare was about 0.2 times the Eddington value.

The long X-ray flare bears all the characteristics of the superbursts discovered recently in five other X-ray burst sources: 4U 1735$``$44, Ser X-1, 4U 1636$``$53, 4U 1820$``$30 and KS 1731$``$260 (Cornelisse et al. 2000, 2001; Wijnands 2001; Strohmayer & Brown 2001; Kuulkers et al. 2001). They have been shown to be thermonuclear events on the surfaces of neutron stars, in which carbon (Cumming & Bildsten 2001; Strohmayer & Brown 2001) or electron capture (Kuulkers et al. 2001) plays the dominant role.

Assuming the persistent luminosity is a direct measure of the mass accretion rate, $`\dot{M}`$, onto the neutron star, one would infer $`\dot{M}`$$``$2$`\times `$1017 g s-1 in GX 3+1 just before the superburst. But since the persistent luminosity is variable by a factor of $``$2 and at the time of the superburst it was at a minimum as seen by the ASM in $``$6 yrs, I assume an average $`\dot{M}`$ of $``$3$`\times `$1017 g s-1. If the ignition of the superburst occurs in a carbon-rich layer, recurrence times of $``$13 yrs are then expected, unless the ignition is prematurely triggered (Strohmayer & Brown 2001). Note that this case is applicable to a neutron star accreting pure helium, such as 4U 1820$``$30. Of the other systems, only 4U 1735$``$44 and 4U 1636$``$53 are known to have hydrogen-rich donors (Augusteijn et al. 1998); this is not as yet clear for GX 3+1. For hydrogen/helium accreting neutron stars it has been shown that carbon may ignite earlier when in an ocean of heavy elements, leading to a smaller recurrence time of $``$2 yrs at the average $`\dot{M}`$ (Cumming & Bildsten 2001). If, on the other hand, the superbursts are due to unstable electron capture by protons with subsequent capture of the resulting neutrons by heavy nuclei, then recurrence times of less than 0.5 yr are expected (Kuulkers et al. 2001). The expected recurrence times in the latter two cases are more or less compatible with seeing one superburst from GX 3+1 in a total effective exposure time of about 1 to 1.5 yrs, as well as seeing two such events 4.7 yrs apart from 4U 1636$``$53 (Wijnands 2001). Effective exposure times can in principle be determined for the other superburst sources, using RXTE/ASM data together with data from other instruments (e.g. BeppoSAX/WFCs); however, this is outside the scope of this paper. Longer and more frequent monitoring of these sources will enable one to discriminate the origin of the superbursts.

###### Acknowledgements.

I thank Al Levine for providing the ASM on-source effective area matrix and discussions related to it, Jean in ’t Zand for comments on an earlier draft, Andrew Cumming and Lars Bildsten for a discussion on the recurrence times, and the referee for useful comments. I acknowledge the use of ASM data products provided by the ASM/RXTE team (http://xte.mit.edu; http://heasarc.gsfc.nasa.gov/docs/xte/asm\_products.html).


text:
# References

1. In this Letter, I shall discuss some phenomenological aspects of the ground state of spontaneously broken theories: the ‘Higgs condensate’. The name itself (as for the closely related gluon, chiral,..condensates) indicates that, in our view, this represents a kind of medium made up by the physical condensation process of some elementary quanta. If this were true, such a vacuum should support long-wavelength density fluctuations. In fact, the existence of density fluctuations in any known medium is a basic experimental fact depending on the coherent response of the elementary constituents to disturbances whose wavelength is much larger than their mean free path. This leads to an universal description, the ‘hydrodynamical regime’, that does not depend on the details of the underlying molecular dynamics. By accepting this argument, and quite independently of the Goldstone phenomenon, the energy spectrum of a Higgs condensate should terminate with an ‘acoustic branch’, say $`\stackrel{~}{E}(𝐩)=c_s|𝐩|`$ for $`𝐩0`$, as for the propagation of sound waves in ordinary media.

Some arguments suggest that, indeed, the vacuum of a ‘pro forma’ Lorentz-invariant quantum field theory may be such kind of medium. For instance, a fundamental phenomenon as the macroscopic occupation of the same quantum state (say $`𝐩=0`$ in some frame) may represent the operative construction of a ‘quantum aether’ . This would be quite distinct from the aether of classical physics, considered a truly preferred reference frame and whose constituents were assumed to follow definite space-time trajectories. However, it would also be different from the empty space-time of special relativity, assumed at the base of axiomatic quantum field theory to deduce the exact Lorentz-covariance of the energy spectrum.

In addition, one should take into account the approximate nature of locality in cutoff-dependent quantum field theories. In this picture, the elementary quanta are treated as ‘hard spheres’, as for the molecules of ordinary matter. Thus, the notion of the vacuum as a ‘condensate’ acquires an intuitive physical meaning. For the same reason, however, the simple idea that all deviations from Lorentz-covariance take place at the cutoff scale may be incorrect. In particular, non-perturbative vacuum condensation may give rise to a hierarchy of scales such that the region of Lorentz-covariance is sandwiched both in the high- and low-energy region.

In fact, in general, an ultraviolet cutoff induces vacuum-dependent reentrant violations of special relativity in the low-energy corner . In the simplest possible case, these extend over a small shell of momenta, say $`|𝐩|<\delta `$, where the energy spectrum $`\stackrel{~}{E}(𝐩)`$ deviates from a Lorentz-covariant form. Since Lorentz-covariance becomes an exact symmetry in the local limit, for very large $`\mathrm{\Lambda }`$, the scale $`\delta `$ is naturally infinitesimal in units of the energy scale associated with the Lorentz-covariant part of the energy spectrum, say $`M`$. By introducing dimensionless quantities, this means $`ϵ\frac{\delta }{M}0`$ when $`t\frac{\mathrm{\Lambda }}{M}\mathrm{}`$ so that the continuum limit can equivalently be defined either as $`t\mathrm{}`$ or $`ϵ0`$. Notice that, formally, $`𝒪(\frac{\delta }{M})`$ vacuum-dependent corrections represent $`𝒪(\frac{M}{\mathrm{\Lambda }})`$ effects which are always neglected when discussing  how Lorentz-covariance emerges at scales much smaller than the ultraviolet cutoff. In this sense, Lorentz-covariance is formally valid in the local limit but, for finite $`\mathrm{\Lambda }`$, there are infinitesimal deviations in an infinitesimal region of momenta that cannot be understood without exploring the physical properties of the vacuum.

For instance, for large but finite $`t`$, does the group velocity $`c_s=\frac{d\stackrel{~}{E}}{d|𝐩|}`$ remain below unity (in units of the light velocity $`c`$) when $`|𝐩|0`$ ? Of course, a group velocity $`c_s>1`$ does not correspond to any ‘classical’ velocity being faster than light. Moreover, when $`t\mathrm{}`$, the relevant wavelengths are $`𝒪(1/ϵ)`$ and thus infinitely long on the physical length scale defined by $`1/M`$. As a consequence, in a strict continuum limit there would be no way to form the sharp wave fronts needed to transfer any type of information. However, what about the cutoff theory ? For instance, if we just play with the numbers, and choose a mass unit $`M=𝒪(10^2)`$ GeV with values of $`t`$ as large as $`𝒪(10^{16})`$, we find conceivable to observe deviations from Lorentz-covariance at scales $`𝒪(10^{33})`$ cm. On the other hand, if such deviations were reentrant with an $`ϵ1/t`$, what about vacuum fluctuations whose wavelengths were larger than a few millimeters ? Should we consider them as ‘infinitely’ long ? In this sense, a non-trivial structure of the vacuum raises delicate issues for which no general answers exist ‘a priori’.

2. After this preliminary discussion, our analysis will start after having understood that, quite independently of the Goldstone phenomenon, there is a gap-less mode of the singlet Higgs field in the spontaneously broken phase of $`\lambda \mathrm{\Phi }^4`$ theories.

The presence of such a gap-less mode reflects the quantum nature of the scalar condensate that cannot be treated as a purely classical c-number field. In fact, either considering the re-summation of the one-particle reducible zero-momentum tadpole graphs  in a given background field or performing explicitely the last functional integration over the strength of the zero-momentum mode of the singlet Higgs field , one finds two different solutions for the inverse zero-4-momentum propagator in the spontaneously broken phase of a (one component) $`\lambda \mathrm{\Phi }^4`$ theory: a) $`G_a^1(p_\mu =0)=M_h^2`$ and b) $`G_b^1(p_\mu =0)=0`$. For the convenience of the reader, we shall briefly repeat the main argument of ref..

When discussing spontaneous symmetry breaking, the starting point is the separation of the scalar field into a constant background and a shifted fluctuation field, namely

$$\mathrm{\Phi }(x)=\varphi +h(x)$$
(1)
In order Eq.(1) to be unambiguous, $`\varphi `$ denotes the spatial average in a large 4-volume $`\mathrm{\Omega }`$

$$\varphi =\frac{1}{\mathrm{\Omega }}d^4x\mathrm{\Phi }(x)$$
(2)
and the limit $`\mathrm{\Omega }\mathrm{}`$ has to be taken at the end.

In this way, the full functional measure can be expressed as

$$[d\mathrm{\Phi }(x)]\mathrm{}=_{\mathrm{}}^+\mathrm{}𝑑\varphi [dh(x)]\mathrm{}$$
(3)
and the functional integration on the r.h.s. of Eq.(3) is over all quantum modes for $`p_\mu 0`$.

After integrating out all non-zero quantum modes, the generating functional in the presence of a space-time constant source $`J`$ is given by

$$Z(J)=_{\mathrm{}}^+\mathrm{}𝑑\varphi \mathrm{exp}[\mathrm{\Omega }(V_{\mathrm{NC}}(\varphi )J\varphi )]$$
(4)
and $`V_{\mathrm{NC}}(\varphi )`$ denotes the usual non-convex (‘NC’) effective potential obtained order by order in the loop expansion. Finally, by introducing the generating functional for connected Green’s functions $`w(J)`$ through

$$\mathrm{\Omega }w(J)=\mathrm{ln}\frac{Z(J)}{Z(0)}$$
(5)
one can compute the field expectation value

$$\phi (J)=\frac{dw}{dJ}$$
(6)
and the zero-momentum propagator

$$G_J(p_\mu =0)=\frac{d^2w}{dJ^2}$$
(7)
In this framework, spontaneous symmetry breaking corresponds to non-zero values of Eq.(6) in the double limit $`J\pm 0`$ and $`\mathrm{\Omega }\mathrm{}`$.

Now, by denoting $`\pm v`$ the absolute minima of $`V_{\mathrm{NC}}`$ and $`M_h^2=V_{\mathrm{NC}}^{\prime \prime }`$ its quadratic shape there, one usually assumes

$$\underset{\mathrm{\Omega }\mathrm{}}{lim}\underset{J\pm 0}{lim}\phi (J)=\pm v$$
(8)
and

$$\underset{\mathrm{\Omega }\mathrm{}}{lim}\underset{J\pm 0}{lim}G_J(p_\mu =0)=\frac{1}{M_h^2}$$
(9)
In this case, the excitations in the broken phase would be massive particles (the conventional Higgs bosons) whose mass $`M_h`$ is determined by the positive curvature of $`V_{\mathrm{NC}}`$ at its absolute minima.

The main result of  is that this conclusion is not true. In fact, at $`\phi =\pm v`$, besides the value $`\frac{1}{M_h^2}`$, one also finds

$$\underset{\mathrm{\Omega }\mathrm{}}{lim}\underset{J\pm 0}{lim}G_J(p_\mu =0)=+\mathrm{}$$
(10)
a result that has no counterpart in perturbation theory. As discussed in ref., the existence of such divergent behaviour admits a simple geometric interpretation in terms of the Legendre transform of $`w(J)`$. Differently from $`V_{\mathrm{NC}}`$, this other definition of effective potential does not provide an infinitely differentiable function in the presence of spontaneous symmetry breaking . Therefore, its left- and right- second derivatives at $`\phi =\pm v`$ do not coincide. In this sense, the existence of a singular zero-4-momentum propagator in the broken phase is a genuine quantum-field theoretical effect, quite independent of any physical interpretation. This will be discussed below.

3. ‘A priori’, the existence of two possible values for the zero-4-momentum propagator implies two possible types of excitations with the same quantum numbers but different energies when the 3-momentum $`𝐩0`$: a massive one, with $`\stackrel{~}{E}_a(𝐩)M_h`$, and a gap-less one with $`\stackrel{~}{E}_b(𝐩)0`$. However, the latter dominates the exponential decay $`e^{\stackrel{~}{E}_b(𝐩)T}`$ of the connected euclidean correlator for $`𝐩0`$. In this sense, the massive excitation is unphysical in the infrared region. Therefore, differently from the simplest perturbative indications, in a (one-component) spontaneously broken $`\lambda \mathrm{\Phi }^4`$ theory there is no energy-gap associated with the ‘Higgs mass’ $`M_h`$, as it would be for a genuine massive single-particle spectrum where the relation

$$\stackrel{~}{E}_a(𝐩)=\sqrt{𝐩^2+M_h^2}$$
(11)
remains true for $`𝐩0`$. Rather, the infrared region is dominated by gap-less collective excitations

$$\stackrel{~}{E}_b(𝐩)c_s|𝐩|$$
(12)
depending on an unknown parameter $`c_s`$ that controls the slope of the spectrum for $`𝐩0`$ and represents the ‘sound velocity’ for the density fluctuations of the scalar condensate. The massive branch, however, can become important at higher momenta. To understand this point, let us explore the analogy with superfluid 4He. This analogy is based on the observation  that, as for the interatomic 4He-4He potential, the low-energy limit of cutoff $`\lambda \mathrm{\Phi }^4`$ is also a theory of quanta with a short-range repulsive core and a long-range attractive tail .

The essential point is that, for superfluid 4He, the existence of two types of excitations was first deduced theoretically by Landau on the base of very general arguments . According to this original idea, there are phonons with energy $`E_{\mathrm{ph}}(𝐩)=v_s|𝐩|`$ and rotons with energy $`E_{\mathrm{rot}}(𝐩)=\mathrm{\Delta }+\frac{𝐩^2}{2\mu }`$. Only later, it was experimentally discovered that there is a single energy spectrum $`E(𝐩)`$ which is made up by a continuous matching of these two different parts. This unique spectrum agrees with the phonon branch for $`𝐩0`$ and agrees with the roton branch at higher momenta.

Although these results are well established, all details of the energy spectrum of superfluid 4He in the matching region are not yet completely understood . For this reason, in our case of the spontaneously broken phase, we shall just extract the main conclusion: the existence of a single energy spectrum $`\stackrel{~}{E}(𝐩)`$ that tends to $`\stackrel{~}{E}_b(𝐩)`$ for $`𝐩0`$ and approaches $`\stackrel{~}{E}_a(𝐩)`$ for larger values of $`|𝐩|`$.

However, for $`\mathrm{\Delta }=\mu M_h`$, the matching between Eqs.(11) and (12) is only possible for $`c_s>1`$, in units of the light velocity $`c`$ (this can easily be checked looking for a possible intersection between Eqs.(11) and (12)). Only in this case, for sufficiently high momenta where $`\stackrel{~}{E}_a(𝐩)<\stackrel{~}{E}_b(𝐩)`$, the gap-less collective modes will become unphysical and the lowest excitations of the vacuum will correspond to the familiar, massive Higgs boson.

Independently of this argument, the idea that $`c_s`$ is actually infinitely larger than the light velocity is also supported by a semi-classical argument due to Stevenson  that we shall briefly report. Stevenson’s argument starts from a perfect-fluid treatment of the Higgs condensate. In this approximation, energy-momentum conservation is equivalent to wave propagation with a squared velocity given by ( $`c`$ is the light velocity)

$$c_s^2=c^2(\frac{𝒫}{})$$
(13)
where $`𝒫`$ is the pressure and $``$ the energy density. Introducing the condensate density $`n`$, and using the energy-pressure relation

$$𝒫=+n\frac{}{n}$$
(14)
we obtain

$$c_s^2=c^2(\frac{𝒫}{n})(\frac{}{n})^1=c^2(n\frac{^2}{n^2})(\frac{}{n})^1$$
(15)
For a non-relativistic Bose condensate of neutral particles with mass $`m`$ and scattering length $`a`$, where $`\frac{na\mathrm{}^2}{m^2c^2}1`$ (in this case we explicitely introduce $`\mathrm{}`$ and $`c`$) one finds

$$=nmc^2+n^2\frac{2\pi a\mathrm{}^2}{m}$$
(16)
so that

$$c_s^2=\frac{4\pi na\mathrm{}^2}{m^2}$$
(17)
which is the well known result for the sound velocity in a dilute hard sphere Bose gas .

On the other hand, in a fully relativistic case, the additional terms in Eq.(16) are such that the scalar condensate is spontaneously generated from the ‘empty’ vacuum where $`n=0`$ for that particular equilibrium density where

$$\frac{}{n}=0$$
(18)
Therefore, in this approximation, approaching the equilibrium density one finds

$$c_s^2\mathrm{}$$
(19)
thus implying that long-wavelength density fluctuations would propagate instantaneously in the spontaneously broken vacuum.

As Stevenson points out , Eq.(19) neglects all possible corrections to the perfect-fluid approximation. These require to introduce the effects of a mean free path $`R_{\mathrm{mfp}}`$ for the elementary constituents. Due to its finite value, in fact, sound waves will stop to propagate at a typical momentum $`|𝐩|\delta \frac{1}{R_{\mathrm{mfp}}}`$ . At this value, the collective modes Eq.(12) become unphysical and for higher momenta the spectrum Eq.(11) applies. The transition corresponds to

$$\sqrt{\delta ^2+M_h^2}c_s\delta $$
(20)
so that for $`c_s\mathrm{}`$, $`\frac{\delta }{M_h}\frac{1}{c_s}0`$. Therefore, $`c_s`$ represents the inverse of the parameter $`ϵ`$ introduced before to characterize the continuum limit. For this reason, we shall replace the condition $`ϵ0`$ with $`c_s\mathrm{}`$.

As anticipated, in the continuum limit, superluminal wave propagation is restricted to the region $`𝐩0`$ and , therefore, one cannot use these waves to form a sharp wave front. As such, there would be no possibility to transfer informations giving rise to violations of causality . These require a group velocity $`\frac{d\stackrel{~}{E}}{d|𝐩|}>1`$ at finite $`|𝐩|`$ that cannot occur due to the change of the energy spectrum from Eq.(12) to Eq.(11).

In this sense, the perfect-fluid value $`c_s=\mathrm{}`$ simulates an exact Lorentz-covariant limit where the energy spectrum maintains its massive form for $`𝐩0`$. Yet, this is not entirely true due to the subtleties associated with the zero-measure set $`𝐩=0`$. This set, in fact, belongs to the range of Eq.(12) and therefore the right $`c_s=\mathrm{}`$ limit is always $`\stackrel{~}{E}(𝐩=0)=0`$ and not $`\stackrel{~}{E}(𝐩=0)=M_h`$. Just for this reason, the correct procedure is to include both branches of the spectrum in the spectral representation of the fluctuation field $`h(x)`$. It may be convenient, however, to separate out the long-wavelength modes Eq.(12), say $`\stackrel{~}{h}(x)`$, from the more conventional massive part of Eq.(11). The observable effects due to $`\stackrel{~}{h}(x)`$ depend crucially on the value of $`c_s`$ and will be discussed below.

4. Let us ignore, for the moment, the previous indication in Eq.(19) and just explore the phenomenological implications of long-wavelength modes in the spectrum as in Eq.(12). Whatever the value of $`c_s`$, these dominate the infrared region so that a general yukawa coupling of the Higgs field to fermions will give rise to a long-range attractive potential between any pair of fermion masses $`m_i`$ amd $`m_j`$

$$U_{\mathrm{}}(r)=\frac{1}{4\pi c_s^2\mathrm{\Phi }^2}\frac{m_im_j}{r}$$
(21)
The above result would have a considerable impact for the Standard Model if we take the value $`\mathrm{\Phi }246`$ GeV related to the Fermi constant. Unless $`c_s`$ be an extremely large number (in units of $`c`$) one is faced with strong long-range forces coupled to the inertial masses of the known elementary fermions that have never been observed. Just to have an idea, for $`c_s=1`$ the long-range interaction between two electrons in Eq.(21) is $`𝒪(10^{33})`$ larger than their purely gravitational attraction. On the other hand, invoking a phenomenologically viable strength, as if $`c_s\mathrm{\Phi }`$ were of the order of the Planck scale, is equivalent to re-obtain nearly instantaneous interactions transmitted by the scalar condensate as in Eq.(19).

Independently of phenomenology, the identification $`c_s\mathrm{\Phi }=M_{\mathrm{Planck}}`$ is also natural  noticing that for $`c_s\mathrm{}`$ the energy-spectrum becomes Lorentz-covariant (with the exception of $`𝐩=0`$). Therefore, in a picture where the ‘true’ dynamical origin of gravity is searched into long-wavelength deviations from exact Lorentz-covariance , it would be natural to relate the limit of a vanishing gravitational strength, $`M_{\mathrm{Planck}}\mathrm{}`$, to the limit of an exact Lorentz-covariant spectrum, $`c_s\mathrm{}`$ .

Returning to more phenomenological aspects, we observe that the potential in Eq.(21) can also be derived as a static limit from the effective lagrangian ($`\stackrel{~}{\sigma }\frac{\stackrel{~}{h}}{\mathrm{\Phi }}`$)

$$_{\mathrm{eff}}(\stackrel{~}{\sigma })=\frac{\mathrm{\Phi }^2}{2}\stackrel{~}{\sigma }[c_s^2\mathrm{\Delta }\frac{^2}{t^2}]\stackrel{~}{\sigma }\stackrel{~}{\sigma }\underset{f}{}m_f\overline{\psi }_f\psi _f$$
(22)
where the free part takes into account the peculiar nature of the energy spectrum Eq.(12). Eq.(22) is useful to represent the effects of $`\stackrel{~}{\sigma }`$ over macroscopic scales as required by its long-range nature. The key-ingredient is the replacement of $`m\overline{\psi }\psi `$ with $`T_\mu ^\mu (x)`$, the trace of the energy-momentum tensor of ordinary matter, a result embodied into the well known relation

$$f|T_\mu ^\mu |f=m_f\overline{\psi }_f\psi _f$$
(23)
This relation allows for an intuitive transition from the quantum to the classical theory. In fact, by introducing a wave-packet corresponding to a particle of momentum $`𝐩`$ and normalization $`d^3𝐱\overline{\psi }\psi =\frac{m}{E(𝐩)}`$ we obtain

$$md^4x\overline{\psi }\psi =m𝑑s$$
(24)
where $`ds=dt\sqrt{1𝐯^2}`$ denotes the infinitesimal element of proper time for a classical particle with velocity $`𝐯`$. Therefore, using the relation

$$\underset{n}{}m_n𝑑s_n=d^4xT_\mu ^\mu (x)$$
(25)
where

$$T_\mu ^\mu (x)\underset{n}{}\frac{E_n^2𝐩_n𝐩_n}{E_n}\delta ^3(𝐱𝐱_n(t))$$
(26)
Eq.(22) is finally replaced by

$$_{\mathrm{eff}}(\stackrel{~}{\sigma })=\frac{\mathrm{\Phi }^2}{2}\stackrel{~}{\sigma }[c_s^2\mathrm{\Delta }\frac{^2}{t^2}]\stackrel{~}{\sigma }\stackrel{~}{\sigma }T_\mu ^\mu $$
(27)
In this way, we get the equation of motion

$$[c_s^2\mathrm{\Delta }\frac{^2}{t^2}]\stackrel{~}{\sigma }=\frac{T_\mu ^\mu }{\mathrm{\Phi }^2}$$
(28)
and it is clear that, for those very large values $`c_s\mathrm{}`$ suggested by the properties of the vacuum, the $`\stackrel{~}{\sigma }`$ Green’s function has practically no retardation effects. In this limit, Eq.(28) reduces to an instantaneous interaction

$$\mathrm{\Delta }\stackrel{~}{\sigma }=\frac{T_\mu ^\mu }{c_s^2\mathrm{\Phi }^2}$$
(29)
of vanishingly small strength. Finally for very slow motions, when the trace of the energy-momentum tensor reduces to the mass density

$$\rho (x)\underset{n}{}m_n\delta ^3(𝐱𝐱_n(t))$$
(30)
one gets, formally, a Poisson equation where, however, the Newton constant $`G_N=M_{\mathrm{Planck}}^2`$ is replaced by the product $`c_s^2\mathrm{\Phi }^2`$

$$\mathrm{\Delta }\stackrel{~}{\sigma }=\frac{1}{c_s^2\mathrm{\Phi }^2}\rho (x)$$
(31)
For this reason, again, one may be tempted to exploit the identification $`M_{\mathrm{Planck}}=c_s\mathrm{\Phi }`$ .

Before concluding, we comment on the momentum $`\delta `$ Eq.(20) associated with the transition between the two branches of the spectrum Eqs.(12) and (11). For $`rR_{\mathrm{mfp}}=\delta ^1`$ the interparticle potential is not a simple $`1/r`$, as for asymptotic distances, but has to be computed from the Fourier transform of the $`h`$field propagator

$$D(r)=\frac{d^3𝐩}{(2\pi )^3}\frac{e^{i𝐩𝐫}}{\stackrel{~}{E}^2(𝐩)}$$
(32)
and depends on the detailed form of the spectrum that interpolates between Eqs.(12) and (11). However, in the Standard Model, for $`M_h=𝒪(\mathrm{\Phi })`$ and $`c_s\mathrm{\Phi }=M_{\mathrm{Planck}}`$, one would predict in any case a length scale $`R_{\mathrm{mfp}}=\delta ^1`$ in the millimeter range. In this framework, the tight infrared-ultraviolet connection embodied in the relation $`\delta \frac{\mathrm{\Phi }^2}{M_{\mathrm{Planck}}}`$ would formally be identical to that occurring in models  with extra space-time dimensions compactified at a size $`R_\mathrm{c}=R_{\mathrm{mfp}}`$.

5. Summarizing: the conventional description of the singlet Higgs field as a purely massive field has to be modified to take into account that the energy spectrum has the form $`c_s|𝐩|`$ for $`𝐩0`$. In this sense, the Higgs condensate is a truly physical medium that can support long-wavelength density fluctuations. However, their velocity $`c_s`$ becomes infinitely large and their wavelengths become infinitely long when approaching the continuum limit. Physically, this corresponds to treat the Higgs condensate in a perfect-fluid approximation.

In a cutoff theory, $`c_s`$ is finite and this amounts to introduce a finite mean free path $`R_{\mathrm{mfp}}\delta ^1`$ for the elementary constituents in the condensate. In this situation, one finds a weak, attractive long-range force proportional to $`\frac{1}{c_s^2}`$ and the energy spectrum acquires its single-particle form $`\sqrt{𝐩^2+M_h^2}`$ at very small momenta $`|𝐩|\delta `$ with $`\delta =𝒪(\frac{M_h}{c_s})`$. There are some arguments for the identification $`c_s\mathrm{\Phi }=M_{\mathrm{Planck}}`$, that suggest the possible relevance of our picture in connection with the problem of gravity . In this case, the length $`R_{\mathrm{mfp}}`$ would be a fundamental scale placed in the millimeter range. Acknowledgements I thank P. M. Stevenson for useful discussions and collaboration.