text:
# Magnetism in the dilute Kondo lattice model

(April 14, 2003)

## Abstract

The one dimensional dilute Kondo lattice model is investigated by means of bosonization for different dilution patterns of the array of impurity spins. The physical picture is very different if a commensurate or incommensurate doping of the impurity spins is considered. For the commensurate case, the obtained phase diagram is verified using a non-Abelian density-matrix renormalization-group algorithm. The paramagnetic phase widens at the expense of the ferromagnetic phase as the $`f`$-spins are diluted. For the incommensurate case, antiferromagnetism is found at low doping, which distinguishes the dilute Kondo lattice model from the standard Kondo lattice model.

Heavy fermion systems have been of great theoretical interest since their discovery some twenty years ago.hf  The central problem posed by heavy fermion materials is to understand the interaction between an array of localized moments (generally $`f`$ \- electrons in lanthanide or actinide ions) and conduction electrons (generally $`s`$ \- or $`d`$ \- band). This situation is well described by an antiferromagnetically coupled Kondo type model.

The solution of Kondo type models is well understood in two limiting cases; the single impurity limitkondo  which can be reduced to a one dimensional problem and solved via Bethe Ansatz, and secondly the Kondo lattice model (KLM), which was solved via bosonizationboso ; graeme  and numerous numerical approachesian ; numerics  in one dimension for half - and partial conduction band - filling. For half - filling the results indicate the existence of a finite spin and charge gap. Accordingly in this case the Kondo lattice model is an insulator with well defined massive solitonic excitations of the spin sector.

For partial conduction band filling, the conduction electrons form a Luttinger liquid, with spin and charge separation.boso ; numerics  The localized spins, however, exhibit ferromagnetism, due to an effective double exchange coupling.graeme ; ian  The double exchange is driving the system toward ferromagnetism, while the fluctuations generated by Kondo singlets compete against this tendency. As a consequence, the paramagnetic to ferromagnetic phase transition is of the quantum order - disorder type, typical to models with an effective random field.graeme  However, for small Kondo coupling and close to half filling an RKKY liquid state and polaronic regime are always present.ian

Beyond these two solvable limits, no rigorous results exist for the intermediate cases, where the number of impurities are neither one, nor equal to the number of sites. This is the focus of our study. We concentrate on the one dimensional case, and start from the Kondo lattice limit introducing impurity spin holes, that is we will be dealing with a dilute Kondo lattice model (DKLM):

$$H=t\underset{i=1,\sigma }{\overset{L1}{}}(c_{i,\sigma }^{}c_{i+1,\sigma }+\mathrm{h}.\mathrm{c}.)+J\underset{i=1}{\overset{L}{}}𝒫𝐒_i^c𝐒_i𝒫,$$
(1)
where $`L`$ is the number of sites and $`t>0`$ is the conduction electron hopping. We measure the Kondo coupling $`J`$ in units of the hopping $`t`$. We denote by $`N_f`$ ($`n_f=N_f/L`$) the number (concentration) of impurities and $`N_c`$ ($`n_c=N_c/L`$) the number (concentration) of conduction electrons. The constraint $`N_fL`$ is imposed by $`𝒫`$, which is an operator that projects out a pre-determined set of $`f`$-spins. $`𝐒_i`$ are spin $`1/2`$ operators for the localized spins, e.g. $`f`$, and $`𝐒_i^c=\frac{1}{2}_{\sigma ,\sigma ^{}}c_{i,\sigma }^{}𝝈_{\sigma ,\sigma ^{}}c_{i,\sigma ^{}}`$ with $`𝝈`$ the Pauli spin matrices and $`c_{i,\sigma }^{}`$, $`c_{i,\sigma }`$ the electron creation and annihilation site operators.

We investigate the behavior of the DKLM both by an analytical approach, based on a standard bosonization scheme, and by numerical calculations. The latter were performed using the newly developed non-Abelian density-matrix renormalization-group (DMRG) algorithm,dmrg  which preserves the total spin and pseudospin symmetry. This choice of basis first of all greatly facilitates the observation of magnetic phases and secondly it gives a dramatic performance improvement compared to the standard DMRG basis.

The bosonization we use takes the standard approachhaldane  by first decomposing the on-site operators into Dirac fields, with spinor components $`\tau =\pm `$ (otherwise known as the right, $`\tau =+`$, and left, $`\tau =`$, movers): $`c_{x,\sigma }_\tau c_{\tau ,x,\sigma }_\tau e^{ik_Fx}\mathrm{\Psi }_{\tau ,\sigma }(x)`$, where $`k_F=\pi n_c/2`$ and we consider the lattice spacing to be unity. Next we bosonize the Dirac fields with $`\mathrm{\Psi }_{\tau ,\sigma }=\mathrm{exp}(i\mathrm{\Phi }_{\tau ,\sigma })/\sqrt{2\pi \lambda }`$, where $`1/\lambda `$ is the ultraviolet cutoff. For the scalar Bose fields, $`\mathrm{\Phi }_{\tau ,\sigma }(x)`$ and its canonical conjugate momenta, $`\mathrm{\Pi }_{\tau ,\sigma }(x)`$, $`\mathrm{\Phi }_{\tau ,\sigma }(x)=_{\mathrm{}}^x𝑑x^{}\mathrm{\Pi }_{\tau ,\sigma }(x^{})`$, we use the standard Mandelstam representation,bosorep  which introduces a momentum cutoff function $`\mathrm{\Lambda }(k)=\mathrm{exp}(\lambda |k|/2)`$ via the Fourier transforms. Thus, the electron field can be represented in terms of collective density operators which satisfy Bose commutation relations:

$$c_{\tau ,x,\sigma }\mathrm{exp}(i\tau k_Fx)\mathrm{exp}i\{\theta _\rho (x)+\tau \varphi _\rho (x)+\sigma [\theta _\sigma (x)+\tau \varphi _\sigma (x)]\}/2,$$
(2)
where the Bose fields (for both $`\nu =\rho ,\sigma `$) are defined by $`\varphi _\nu /\theta _\nu =i(\pi /N)_{k0}e^{ikx}[\nu _+(k)\pm \nu _{}(k)]\mathrm{\Lambda }(k)/k`$, where $`\varphi _\nu `$ are the number fields and $`\theta _\nu `$ the current fields. The charge (holon) and spin (spinon) number fluctuations are defined as $`\rho _\tau (k)=_\sigma \rho _{\tau ,\sigma }(k)`$, and $`\sigma _\tau (k)=_\sigma \sigma \rho _{\tau ,\sigma }(k)`$. This type of Bose representation provides a non-perturbative description of the conduction electrons in terms ofhaldane  holons and spinons. We will neglect for the moment all the rapidly oscillating (umklapp) terms. These will give a contribution only at half filling, i.e. $`n_c=n_f`$, and will be analyzed later on. Thus, the bosonized form of DKLM is:

$`H`$ $`=`$ $`{\displaystyle \frac{v_F}{4\pi }}{\displaystyle \underset{j,\nu }{}}\left\{\mathrm{\Pi }_\nu ^2(j)+[_x\varphi _\nu (j)]^2\right\}+{\displaystyle \frac{J}{2\pi }}{\displaystyle \underset{j}{}}[_x\varphi _\sigma (j)]S_j^z`$ (3)
$`+`$ $`{\displaystyle \frac{J}{4\pi \lambda }}{\displaystyle \underset{j}{}}\{\mathrm{cos}[\varphi _\sigma (j)]+\mathrm{cos}[2k_Fj+\varphi _\rho (j)]\}(e^{i\theta _\sigma (j)}S_j^++\mathrm{h}.\mathrm{c}.)`$
$``$ $`{\displaystyle \frac{J}{2\pi \lambda }}{\displaystyle \underset{j}{}}\mathrm{sin}[\varphi _\sigma (j)]\mathrm{sin}[2k_Fj+\varphi _\rho (j)]S_j^z.`$
This equation has the same form as for a standard Kondo latticeboso ; graeme  except i) we have to keep in mind that the impurity spin, i.e. terms containing $`S_j^z`$, $`S_j^+`$ and $`S_j^{}`$, contributes only if there is an $`f`$ spin at site $`j`$, and ii) the even cutoff function $`\mathrm{\Lambda }(k)`$, defined in Eq. (2), satisfying $`\mathrm{\Lambda }(k)1`$ for $`|k|<1/\lambda `$ and $`\mathrm{\Lambda }(k)0`$ otherwise, is needed in the Bose fields to ensure that delocalized conduction electrons are described. Delocalization is essential to describe ferromagnetism,graeme ; ian  since ferromagnetism in the Kondo lattice models is due to the double exchange, which only requires that $`N_c<N_f`$.zener

In this situation each electron has on average more than one localized spin to screen, and since hopping between localized spins is energetically most favorable for electrons which preserve their spin as they hop, this tends to align the underlying localized spins.zener  This also means that double exchange will vanish if the distance between impurity spins is larger than $`\lambda `$. At lengths beyond $`\lambda `$ the electrons will behave as collective density fluctuations, as usual in one dimensional systems.haldane  Hence, $`\lambda `$ measures the effective range of the double exchange, and in principle it is a function of $`J`$, $`n_c`$ and $`n_f`$.

The most straightforward method to determine an ordering of the localized spins is by applying a unitary transformation $`\stackrel{~}{H}=e^{\widehat{\mathrm{S}}}He^{\widehat{\mathrm{S}}}`$. We choose the transformation which changes to a basis of states in which the conduction electron spin degrees of freedom are coupled directly to the localized spins: $`\widehat{\mathrm{S}}=i\frac{J}{2\pi v_F}_j\theta _\sigma (j)S_j^z`$. We perform the unitary transformation up to infinite order, so there is not any artificial truncation error generated (for details see Ref. graeme, ). In the new transformed basis the double exchange interaction leading to ferromagnetism is clearly exhibited and we obtain the effective Hamiltonian for the localized spins:

$`H_{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \frac{J^2}{2\pi ^2v_F}}{\displaystyle \underset{i,j}{}}{\displaystyle \frac{\lambda }{\lambda ^2+(ij)^2}}S_i^zS_j^z`$ (4)
$`+`$ $`{\displaystyle \frac{J}{2\pi \lambda }}{\displaystyle \underset{i}{}}\{\mathrm{cos}[K(i)]+\mathrm{cos}[2k_Fi]\}S_i^x`$
$``$ $`{\displaystyle \frac{J}{2\pi \lambda }}{\displaystyle \underset{i}{}}\mathrm{sin}[K(i)]\mathrm{sin}[2k_Fi]S_i^z.`$
The $`K(i)`$ is a term introduced by the unitary transformation, which counts all the $`S_i^z`$’s to the right of the site $`i`$ and subtracts from those to the left of $`i`$: $`K(i)=(J/2v_F)_{l=1}^{\mathrm{}}(S_{i+l}^zS_{il}^z)`$. This term gives the crucial difference between KLM and DKLM, as will be explained later on. The most important term in Eq. (4) is the first one, which clearly shows that a ferromagnetic coupling emerges for DKLM. This coupling is non-negligible for $`N_c<N_f`$ and $`ij\lambda `$ and its strength will decrease with increasing distance between impurity spins.

For $`N_c<N_f`$, the physical picture given by Eq. (4) will be crucially different if the lattice of impurity spins contains commensurate or incommensurate array of holes. Hence, we analyze these two cases separately. If we have a commensurate doping of the impurity spins, then we can approximate the ferromagnetic term in the usual way by taking $`1/n_f`$ for the shortest average distance between $`f`$ spins: $`\{J^2n_f^2\lambda /[2\pi ^2v_F(1+n_f^2\lambda ^2)]\}_iS_i^zS_{i+1/n_f}^z`$. Lattice sites which are not occupied by $`f`$ spins are inert and do not contribute to the ferromagnetic phase. This was verified by DMRG: the calculated $`ff`$ spin correlation functions behave similarly as those of the normal KLM. The empty sites are inert. The $`f`$-structure factor has the usual peak at $`k/\pi =N_c/N_f`$ for low $`J`$, hence in the commensurate case the DKLM behaves similarly to the standard KLM model.

To understand the behavior of the second and third term from Eq. (4), we notice that $`K(i)`$ is vanishingly small for the commensurate case, as the number of $`f`$ impurity spins to the left and to the right of a given site $`i`$ is the same. So the effective Hamiltonian will reduce to the random transverse field Ising model, as in the KLM.graeme  The randomness is generated by $`(1+\mathrm{cos}[2k_Fi])`$ at large distances and it is driven by a cosine distribution, similarly to spin-glasses.abrikosov  To determine the phase transition we need the dependence $`\lambda =\lambda (J,n_c,n_f)`$, which however is very difficult to determine and as such we use the low density value $`\lambda \sqrt{2/J}`$ close to criticality, similarly to previous works.graeme ; shimoi  In this way we obtain the critical phase transition line to be: $`J=\pi \mathrm{sin}(\pi n_c/2)/[1\pi \mathrm{sin}(\pi n_c/2)/n_f^2]`$, which represents a quantum order - disorder transition with variable exponents.graeme  However, this ferromagnetic phase disappears for larger distances between impurities because, as mentioned earlier, the double exchange interaction vanishes if the average distance between impurity spins, $`1/n_f`$, is larger than $`\lambda `$. This is very important because it ensures that the single impurity limit, $`n_f0`$, is free of ferromagnetism, as it should be.

The incommensurate case is more difficult than the commensurate case. The reason is that in the low concentration limit the properties of DKLM will be very much dependent on the random distribution of $`f`$ spins. We may observe phase separation or clusterization processes in this case. In this limit, where the average distance between impurities is very large, then the single impuritykondo  approximation seems natural. However, if we look at small doping of $`f`$ electrons only, then the main difference compared to the commensurate limit studied previously is that the $`K(i)`$ term, in Eq. (4), is not negligible anymore. The impurity $`f`$ spins are no longer equally distributed to the left and right of a given site $`i`$, hence $`K(i)(1)^i(J/2v_F)`$, which gives rise to a staggered field. The properties of Eq. (4) are then given by the staggered field Ising model, which gives an antiferromagnetic ordering. This antiferromagnetic ordering of the impurity spins represents a new element in DKLM compared to the Kondo lattice. This corresponds to the soliton lattice obtained by Schlottmann in a dynamical mean-field treatment of the three dimensional dilute Kondo lattice.schlottmann

Similar behavior also occurs above half-filling, i.e. $`N_c>N_f`$, where double exchange (as shown previously) does not appear. But bosonization still works: the effective Hamiltonian reduces to the second and third terms of Eq. (4), from which the most dominant term, for low doping of impurity spins, as in the case described previously, is a staggered $`S_i^z`$ field. As the first term in Eq. (4) is missing in this case, the only fluctuation which can destroy a locked staggered order is $`S_i^x`$. For large $`J`$ ($`4J`$) the staggered order wins. While for smaller values of $`J`$ the systems will be disordered.

As we approach half filling from both sides, the bosonization approach breaks down as the strongly oscillating (umklapp) fields start dominating. The DKLM will undergo a metal-insulator transition as in a standard quantum sine-Gordon modeloldpaper  by dynamical mass generation. A Spin gap will also appear. This can be understood easily, because the half filled DKLM is equivalent to the quarter filled periodic Anderson model, which has a antiferromagnetic order.shimoi  The only difference from the Kondo lattice is that the massive solitons obtained for DKLM are of Su, Schrieffer and Heeger type.SSH

To confirm the previously obtained magnetic phases, we performed non-Abelian DMRG analysis of the DKLM model in the commensurate case, for both the $`N_c<N_f`$ and $`N_c>N_f`$, where ferro- and antiferromagnetism, respectively exist. For the $`N_c<N_f`$ case, we have investigated several commensurate dilution patterns of the form: “0 0 $`f`$ 0 0 $`f`$ 0 0 $`f`$ $`\mathrm{}`$”; “0 $`f`$ 0 $`f`$ 0 $`f`$ 0 $`f`$ $`\mathrm{}`$”; and “0 $`f`$ $`f`$ 0 $`f`$ $`f`$ 0 $`f`$ $`f`$ $`\mathrm{}`$” in a 64-site long chain. While in the $`N_c>N_f`$ limit we have studied the $`n_f=0.8`$ localized spin filling (to be in accordance with the bosonization requirement of low dilution) on a 80-site long DKLM chain, i.e., we investigated the pattern: “0 $`f`$ $`f`$ $`f`$ $`f`$ 0 $`f`$ $`f`$ $`f`$ $`f`$ $`f`$ $`\mathrm{}`$”.

These patterns were selected in such a way that the chain middle-point reflection symmetry was preserved to accelerate the calculations. There was only one exception: the “0 $`f`$ 0 $`f`$ 0 $`f`$” chain has an impurity in the middle: “0 $`f`$ 0 $`f`$ $`f`$ 0 $`f`$ 0”, but its effect is rather small compared to our final errors and so it was neglected.

In the DMRG calculations careful error and convergence analysis were used, and we extrapolated the energy linearly to zero truncation error (we saw no quadratic terms large enough to affect the fit). For each given dilution pattern, filling $`N_c/N_f`$, and interaction constant $`J`$, we used several DMRG sweeps of between 200 and 400 $`SO(4)`$-symmetric states.

For $`N_c<N_f`$ ferromagnetism appears at large $`J`$ values, see Fig. 1. A point on the phase diagram shown in Fig. 1 is judged to be ferromagnetic if the extrapolated energy of the spin $`S_{\mathrm{max}}`$ run is lower than the spin zero energy. This energy difference to the spin singlet excited state can be calculated directly using the $`SO(4)`$ basis set of the non-Abelian DMRG. The phase transition line can be determined with high accuracy if one plots the energy gap between the spin $`S_{\mathrm{max}}`$ and spin zero states as a function of $`J`$ \- the gap rapidly decreases as the transition is approached. It should be mentioned that, as well as in the standard Kondo case,ian  there also exists additional ferromagnetic phases inside the paramagnetic region for dilute Kondo chains.

In the opposite limit, i.e., $`N_c>N_f`$, we have confirmed the existence of antiferromagnetism by calculating the spin structure factor $`S(k)`$. As can be seen in Fig. 2, for small $`J`$ and away from half filling the peak in $`S(k)`$ is at $`2k_F`$ (similarly to KLM ian ), i.e., the DKLM is a disordered paramagnet. The dominant $`2k_F`$ backscattering processes are manifest of a system of free localized spins embedded in effective fields determined by conduction electron scattering. The half filled case, i.e., $`n_c=1`$ as in the KLM model halffilled  is always an antiferromagentic insulator. The tendency of Kondo lattice models to form an antiferromagnetic insulator at half filling for any nonzero $`J`$ is a consequence of the Kondo effect.halffilled  As we increase $`J`$, the distinctive antiferromagnetic peak at $`k=\pi `$ appears away from half filling (see Fig. 2). As we further increase $`J`$ the $`2k_F`$ peak broadens and decreases in strength while the antiferromagentism becomes more robust, at $`J=5`$ it appers already at $`n_c0.925`$.

In conclusion, we have studied the dilute Kondo lattice model in one dimension both numerically, using DMRG, and analytically, with a standard bosonization approach. We have derived an effective Hamiltonian for the $`f`$ spins, which accounts for the appearance of a ferromagnetic phase seen with a commensurate dilution pattern of the impurity spin array. The paramagnetic-ferromagnetic phase transition shifts to higher coupling $`J`$ values as the $`f`$-spins of the chain are diluted, in agreement with the numerical DMRG calculation. We have also shown that for incommensurate dilution or above half filling, i.e., $`N_c>N_f`$ for low doping of impurity spins, antiferromagnetism is found. This distinguishes the dilute Kondo lattice model from the standard Kondo lattice model.

Work in Australia was supported by the Australian Research Council and Department of Industry, Science and Resources. In Sweden by The Swedish Natural Science Research Council, The Swedish Foundation for International Cooperation in Research and Higher Education (STINT), and The Göran Gustafsson foundation. Some of the numerical calculations were performed at National Facility of the Australian Partnership for Advanced Computing, via an award from the ANU Supercomputer Time Allocation Committee.


format: text

text:
# Analysis of energy- and time-dependence of supernova shock effects on neutrino crossing probabilities

## I Introduction

The realization that the supernova shock propagation can affect neutrino flavor transitions Schi  a few seconds after core bounce Wils  is gaining increasing attention in the recent literature on supernova neutrinos Raff ; Beac ; Lead ; Reli ; Spin ; Shoc ; Chou ; Luna . Indeed, time-dependent variations of the neutrino potential in matter Matt  along the supernova shock profile can leave an interesting imprint on the energy and time structure of the (anti)neutrino signal at the Earth Schi ; Shoc ; Luna . Observation of such possible effects, although very challenging from an experimental viewpoint, would open a unique opportunity to study some aspects of the supernova shock dynamics in “real time.” Moreover, the strong sensitivity of such effects to the neutrino oscillation parameters Schi ; Shoc ; Luna  might, in principle, provide us with additional constraints (or hints, at least) on the neutrino mass spectrum and mixing angles.

From a conservative viewpoint, it should be stressed that current simulations of core-collapse supernovae, despite remarkable efforts, may still require substantial physics improvements Jank . It is not excluded that more refined simulations might significantly modify, e.g., the main features of shock profile discussed in Schi , and especially its propagation velocity and density gradient, which govern the structure of the neutrino crossing probabilities Schi ; Shoc ; Luna . Moreover, the time and energy dependence of the source (anti)neutrino fluxes might be more complicated and uncertain Unce ; Flux  than is customarily assumed, making it more difficult to identify shock-related signals. All these potential systematic uncertainties in the basic physics ingredients might be large enough to weaken the significance of shock-related signals, even with hypothetically large experimental statistics. However, despite all these caveats, the stakes are so high that further investigations on possible shock effects on supernova neutrinos are widely justified, in our opinion.

The purpose of this work is to investigate in some detail the energy- and time-dependence of the neutrino crossing probabilities along the supernova shock. We assume the active $`3\nu `$ oscillation scenario (Sec. II), and construct an empirical parameterization for the shock profile (Sec. III). Analytical calculations are shown to reproduce well (and to be much more convenient than) numerical approaches to the overall crossing probabilities, and are used to break down the multiple transition structure along the shock density profile, both in the energy domain (Sec. IV) and in the time domain (Sec. V). Finally, applications are given in terms of antineutrino event spectra observable at the Earth through inverse beta decay (Sec. VI). Conclusions and prospects for further work are given in Sec. VII.

The bottom line of our work is that simple, analytical calculations of the crossing probabilities along the supernova shock profile appear to be adequate to study the structure of neutrino flavor transitions in the energy and time domain, for all practical purposes. The specific probability values estimated in this work, however, should not be taken too literally, being based on a simplified description of the matter density profile, which is subject to change as more detailed supernova shock simulations are performed and become publicly available.

## II $`3\nu `$ notation and analytical approximations

Throughout this work, we consider only flavor transitions among the three known active neutrinos $`\nu _e`$, $`\nu _\mu `$, and $`\nu _\tau `$. In the following, we set the notation and briefly describe the analytical expressions for the crossing probabilities.

### II.1 Notation for $`3\nu `$ kinematics and dynamics

The $`3\nu `$ squared mass spectrum is parameterized as

$$(m_1^2,m_2^2,m_3^2)=(\frac{\delta m^2}{2},+\frac{\delta m^2}{2},\pm \mathrm{\Delta }m^2),$$
(1)
where the sign before $`\mathrm{\Delta }m^2`$ distinguishes the cases of normal ($`+`$) and inverted ($``$) hierarchy. The mixing matrix $`U`$ is defined in terms of three rotation angles $`\theta _{ij}`$, ordered as for the quark mixing matrix Hagi ,

$$U=U(\theta _{12},\theta _{13},\theta _{23}).$$
(2)
When needed, we take

$`\delta m^2/\mathrm{eV}^2`$ $`=`$ $`7.3\times 10^5,`$ (3)
$`\mathrm{sin}^2\theta _{12}`$ $`=`$ $`0.315,`$ (4)
corresponding to the so-called LMA-I best fit to the solar and reactor neutrino data LMAI . We also take

$$\mathrm{\Delta }m^2=3\times 10^3\mathrm{eV}^2,$$
(5)
close to the best fit to the atmospheric and accelerator neutrino data detected in Super-Kamiokande Deco . The value of $`\mathrm{sin}^2\theta _{23}`$ is basically irrelevant for supernova neutrinos, since it corresponds to an unobservable rotation in the $`(\nu _\mu ,\nu _\tau )`$ flavor subspace. Concerning the mixing parameter $`\mathrm{sin}^2\theta _{13}`$, we will pick up several representative values, well below the current $`3\sigma `$ upper limit ($`\mathrm{sin}^2\theta _{13}<0.05`$ Deco ).

The mass gap hierarchy ($`\delta m^2/\mathrm{\Delta }m^21`$) and the smallness of $`\mathrm{sin}^2\theta _{13}`$ guarantee, to a very good approximation, the factorization of the $`3\nu `$ neutrino dynamics into a “low” ($`L`$) and a “high” ($`H`$) $`2\nu `$ subsystem, whose dominant oscillation parameters are $`(\delta m^2,\theta _{12})`$ and $`(\mathrm{\Delta }m^2,\theta _{13})`$, respectively (see KuPa ; Ours  and references therein). The corresponding neutrino wavenumbers are defined as

$`k_L`$ $`=`$ $`\delta m^2/2E,`$ (6)
$`k_H`$ $`=`$ $`\mathrm{\Delta }m^2/2E,`$ (7)
where $`E`$ is the neutrino energy.

Potentially large matter effects Matt  are expected when either $`k_L`$ or $`k_H`$ take values close to the neutrino potential $`V`$ in matter,

$$V=\sqrt{2}G_FN_e(x),$$
(8)
where $`N_e`$ is the electron density at the supernova radius $`x`$. In appropriate units,

$$\frac{V(x)}{\mathrm{eV}^2/\mathrm{MeV}}=7.6\times 10^8Y_e\frac{\rho (x)}{\mathrm{g}/\mathrm{cm}^3},$$
(9)
where $`Y_e`$ is the electron/nucleon number fraction (that we assume equal to 1/2), while $`\rho (x)`$ is the radial mass density profile. We conventionally take $`k_L`$, $`k_H`$, and $`V`$ to be positive (as in Ref. Ours ), irrespective of the neutrino type ($`\nu `$ or $`\overline{\nu }`$) and hierarchy (normal or inverted).

The $`V(x)`$ profile in supernovae is often approximated by a static power law, $`V(x)x^n`$ with $`n3`$. In the presence of a shock wave, however, such static profile can be significantly modified Schi ; Wils . In particular the shock wave, while propagating outwards at supersonic speed, leaves behind a rarefaction zone, followed by a high-density region and by a sharp drop of density (down to the static value) at

$$x_s=\mathrm{shock}\mathrm{front}\mathrm{radius}.$$
(10)
This abrupt change of density occurs over a microscopic length scale governed by the mean free path of ions and electrons, which, for our purposes, can be effectively taken as zero.111The authors of Schi  correctly recognized that the apparent smoothness of numerical shock front profiles Wils ; Shoc  is actually an artifact of hydrodynamical simulations, which tend to widen the true (step-like) profile over several numerical grid zones. Correspondingly, at $`x=x_s`$ the matter potential $`V(x)`$ drops from

$$V^+=\underset{xx_s^{}}{lim}V(x)$$
(11)
to

$$V^{}=\underset{xx_s^+}{lim}V(x),$$
(12)
with a typical ratio Schi

$$\xi =\frac{V^+}{V^{}}10.$$
(13)
Given the peculiar step-like features of the matter potential at the shock front, in the following we discuss the analytical expressions of the crossing probabilities for $`xx_s`$ and $`x=x_s`$ separately. We then study the general case of multiple crossings along the neutrino trajectory.

### II.2 Crossing probability at $`xx_s`$

As previously mentioned, potentially large matter effects are expected when $`k_{L,H}V`$. In particular, if the condition

$$V(x_H)=k_H$$
(14)
occurs at some point $`x_H(x_s)`$, there is a finite probability $`P_H`$ of crossing between effective mass eigenstates in the $`H`$ subsystem, which is given to a good approximation by the so-called double exponential formula Petc

$$P_H=\frac{\mathrm{exp}(2\pi r_Hk_H\mathrm{cos}^2\theta _{13})1}{\mathrm{exp}(2\pi r_Hk_H)1}$$
(15)
for both neutrinos and antineutrinos (see Ours  and related bibliography therein), with

$$r_H=\left|\frac{d\mathrm{ln}V(x)}{dx}\right|_{x=x_H}^1.$$
(16)
For all the $`V(x)`$ profiles and $`\mathrm{sin}^2\theta _{13}`$ values that will be used in the following, the above formula reduces to the well-known Landau-Zener (LZ) limit

$$P_H\mathrm{exp}(2\pi r_Hk_H\mathrm{sin}^2\theta _{13}),$$
(17)
at any $`\nu `$ energy of practical interest.

In principle, the occurrence of the condition $`k_L=V(x_L)`$ at some radius $`x_L`$ could also induce neutrino level crossing in the $`L`$ subsystem with probability $`P_L`$ (generally different for neutrinos and antineutrinos; see, e.g., Digh ; Ours ). However, for $`\mathrm{sin}^2\theta _{12}`$ as large as in Eq. (4) it turns out that the $`L`$-transition is adiabatic,

$$P_L(\nu )0P_L(\overline{\nu }),$$
(18)
at all points $`xx_s`$.

### II.3 Crossing probability at $`x=x_s`$

The analytical approximations in Eqs. (17) and (18) break down at $`x=x_s`$, where the density gradient explodes (and $`r_H0`$). In this case, the crossing probability $`P_s`$ is determined by the conservation of flavor across the $`V(x)`$ discontinuity, with the well-known result (see, e.g., KuPa ):

$$P_s=\mathrm{sin}^2(\theta _m^+\theta _m^{}),$$
(19)
where $`\theta _m^+`$ and $`\theta _m^{}`$ are the effective mixing angles in matter immediately before $`(V=V^+)`$ and after $`(V=V^{})`$ the shock front position $`x_s`$, respectively. For the $`H`$ transition, such angles are defined by

$$\mathrm{cos}2\theta _m^\pm =\frac{\mathrm{cos}2\theta _{13}V^\pm /k_H}{\sqrt{(\mathrm{cos}2\theta _{13}V^\pm /k_H)^2+(\mathrm{sin}2\theta _{13})^2}},$$
(20)
for both neutrinos and antineutrinos.222We remind the reader that, as proven in Ours , $`P_H(\nu )=P_H(\overline{\nu })`$, independently of the functional form of $`P_H`$. In the limit of very small $`\theta _{13}`$, the above equations imply $`P_s1`$ for $`k_H[V^{},V^+]`$, and $`P_s0`$ otherwise. Notice that the extremely nonadiabatic limit $`P_s1`$ would also be obtained from Eq. (15) for $`\theta _{13}0`$ and $`r_H0`$. However, the tempting approximation $`P_s=1`$ (for $`k_H[V^{},V^+]`$) does not work well for $`\mathrm{sin}^2\theta _{13}`$ as large as $`O(10^3)`$$`O(10^2)`$, the “top-hat” behavior for $`P_s`$ being significantly smeared out \[i.e., a tail $`P_s0`$ is developed for $`k_H`$ outside the range $`[V^{},V^+]`$, especially for $`k_H\stackrel{<}{}V^{}`$, as easily understandable on the basis of Eq. (20)\]. Therefore, we adopt the exact Eqs. (19) and (20) to calculate the $`H`$ crossing probability at the shock front, without further approximations.

Concerning the analogous crossing probability in the $`L`$ sector, similar expressions apply, modulo the substitutions $`\theta _{13}\theta _{12}`$ and $`k_Hk_L`$ in Eqs. (19) and (20).333Contrary to the $`H`$-transition case, it is $`P_L(\nu )P_L(\overline{\nu })`$ in general. $`P_L(\overline{\nu })`$ can be obtained from $`P_L(\nu )`$ through the further substitution $`V^\pm /k_LV^\pm /k_L`$ Ours . The strong nonadiabatic character of the transition at the shock front can now lead to $`P_L0`$ in some corners of the parameter space \[in contrast to Eq. (18)\]. However, as we shall discuss at the end of Sec. IV C, such parameter values happen to have little phenomenological relevance. Therefore, unless otherwise noted, we shall formally take $`P_L0`$ not only for $`xx_s`$ \[Eq. (18)\], but also for $`x=x_s`$.

### II.4 Analytical expressions for multiple crossings

For $`P_L0`$, the formal expression for the $`\nu _e`$ survival probability $`P_{ee}`$ Digh ; Ours  is exceedingly simple (up to Earth matter effects):444Examples of Earth matter effects will be illustrated in Sec. VI

$$P_{ee}\{\begin{array}{cc}\mathrm{sin}^2\theta _{12}P_H\hfill & (\nu ,\mathrm{normal}),\hfill \\ \mathrm{cos}^2\theta _{12}\hfill & (\overline{\nu },\mathrm{normal}),\hfill \\ \mathrm{sin}^2\theta _{12}\hfill & (\nu ,\mathrm{inverted}),\hfill \\ \mathrm{cos}^2\theta _{12}P_H\hfill & (\overline{\nu },\mathrm{inverted}),\hfill \end{array}$$
(21)
where “normal” and “inverted” refer to the hierarchy type. In the above equation, we have neglected small additive terms of $`O(\mathrm{sin}^2\theta _{13})`$, which are not relevant for our discussion. From the above equations, it appears that $`P_H`$ can modulate the (otherwise constant) survival probability of $`\nu _e`$ in normal hierarchy and of $`\overline{\nu }_e`$ in inverted hierarchy, thus providing an important handle to solve the current hierarchy ambiguity.

In the presence of a propagating shock wave, the calculation of $`P_H`$ is not entirely trivial. For the shock profiles studied in Schi ; Shoc , the condition in Eq. (14) can occur at up to three points $`x_ix_s(i=1,2,3)`$:

$$V(x_i)=k_H(i=1,2,3),$$
(22)
In the most general case, two such points belong to the rarefaction zone ($`x_1<x_2<x_s`$) and one to the static region below the shock front $`(x_3>x_s)`$. The corresponding crossing probabilities are denoted as $`P_1`$, $`P_2`$, and $`P_3`$. The $`P_i`$’s can be analitically calculated as in Eq. (17), with the inverse log-derivative $`r_i`$ \[Eq. (16)\] evaluated at the point $`x_i`$. In addition, one should consider the crossing probability $`P_s`$ at $`x=x_s`$, as given by Eqs. (19) and (20). Therefore, the global $`P_H`$ must be generally constructed in terms of four crossing probabilities $`P_1`$, $`P_2`$, $`P_s`$, and $`P_3`$, occurring at the points

$$x_1<x_2<x_s<x_3,$$
(23)
respectively.

If we assume that all relative neutrino phases can be averaged away, and that the four transitions can be exactly factorized, the overall crossing probability $`P_H`$ can be defined in terms of a simple matrix equation KuPa

$$\left(\begin{array}{cc}1P_H& P_H\\ P_H& 1P_H\end{array}\right)=\underset{i=1,2,s,3}{}\left(\begin{array}{cc}1P_i& P_i\\ P_i& 1P_i\end{array}\right),$$
(24)
whose solution is

$`P_H`$ $`=`$ $`P_1+P_2+P_s+P_32(P_1P_2+P_1P_s+P_1P_3+P_2P_s+P_2P_3+P_sP_3)`$ (25)
$`+4(P_1P_2P_s+P_1P_2P_3+P_1P_sP_3+P_2P_sP_3)8P_1P_2P_sP_3.`$
We will refer to this simple equation for our analytical calculations of $`P_H`$.

Although the smallness of $`\theta _{13}`$ implies that each of the four transitions $`P_{1,2,s,3}`$ is (in general) well localized and can thus be separated from the others, the intrinsic accuracy of the factorization in Eq. (24) is not obvious a priori in the whole relevant parameter space. There are zones of the shock profile where two transition points can become very close and eventually merge, or where the local density shape can be strongly different from the static, quasi-cubic power law which was used in Ours  to test the analytical recipe in of Eq. (15). Moreover, in the presence of multiple (interfering) transition amplitudes, it is worth checking that reasonable phase averaging effectively reproduces the incoherent factorization implicit in Eq. (24).

For such reasons, we think it is useful to compare the analytical calculations of $`P_H`$ \[based on Eqs. (16) and (17) for $`x=x_{1,2,3}`$, and on Eqs. (19) and (20) for $`x=x_s`$\] with the results of a numerical (Runge-Kutta) evolution of the neutrino flavor propagation equations along representative shock density profiles. To our knowledge, such reassuring comparison has not been explicitly performed in the available literature.555The authors of Refs. Shoc  and Schi  seem to have used a fully numerical and a mixed numerical+analytical approach to the neutrino evolution equations, respectively. The authors of Ref. Luna  discuss only qualitatively the analytic form of $`P_3`$, in the hypothesis that $`P_10P_2`$. In order to do so in both the time and energy domain (Secs. IV and V, respectively), we need a parameterization of the shock wave profile which is (almost everywhere) continuous in both the radius $`x`$ and in post-bounce time $`t`$. Unfortunately, a similar parameterization (or, equivalently, a dense numerical table) is not publicly available; only some representative profiles are graphically reported in Schi  and Shoc .666Notice that the bilogarithmic scales used in Schi ; Shoc  make it difficult to recover potentially interesting shock details by graphical reduction. Moreover, only the profiles in Schi  exhibit (by construction) the correct step-like profile at $`x=x_s`$. In the following section, we thus introduce a simplified parametrization of the supernova density profile which is continuous in the $`(x,t)`$ variables \[excepting the shock-front discontinuity point\]. Such empirical profile captures the main qualitative features of the time and radial dependence of the (supposed spherically symmetric) shock propagation, allowing us to perform a meaningful comparison of numerical and analytical results, and to break down explicitly the analytical calculation of $`P_H`$ into its four components $`P_1`$, $`P_2`$, $`P_s`$, and $`P_3`$. Needless to say, such a simplified profile can and must be improved when more informative or refined supernova shock simulations become available.

A final remark is in order. In this work we are mainly concerned with the factorization of transition within the $`H`$ subsystem, since the factorization between the $`L`$ and $`H`$ subsystems is guaranteed for $`\delta m^2`$ as low as in Eq. (3). However, small corrections to the $`L`$-$`H`$ factorization might arise for values of $`\delta m^2\mathrm{few}\times 10^4`$ eV2, not considered here (see, e.g., KuPa  for earlier discussions and Luna ; LimS  for more recent considerations). The analysis of these possible small corrections is beyond the scope of this work.

## III Empirical parametrization of the Shock density profile

In this section we introduce a simplified, empirical parametrization of the shock density profile $`\rho `$, which reproduces the main features of the graphical profile in Schi  and, at same time, is continuous in both $`x`$ and $`t`$ \[excepting the case $`x(t)=x_s(t)`$\]. In this way, numerical and analytical calculations of $`P_H`$ can continuously cover the relevant parameter space.

For post-bounce times $`t\stackrel{<}{}1\mathrm{s}`$, shock effects Schi ; Shoc  take place typically at $`V>k_H`$ and do not significantly affect the $`H`$ subsystem. The density profile can thus be effectively approximated by its static limit $`\rho _0`$ as taken from Schi :

$$t\stackrel{<}{}1\mathrm{s}\frac{\rho _0(x)}{\mathrm{g}/\mathrm{cm}^3}10^{14}\left(\frac{x}{\mathrm{km}}\right)^{2.4}.$$
(26)
For post-bounce times $`t\stackrel{>}{}1\mathrm{s}`$, the condition $`Vk_H`$ can be fulfilled, and the shock profile can thus modulate $`P_H`$. We characterize the profile in terms of the shock front position $`x_s`$ and of its shape variation $`f(x)`$ (with respect to $`\rho _0`$) for $`xx_s`$. Formally, one can write

$$t\stackrel{>}{}1\mathrm{s}\rho (x)=\rho _0(x)\{\begin{array}{cc}\xi f(x)\hfill & ,xx_s,\hfill \\ 1\hfill & ,x>x_s,\hfill \end{array}$$
(27)
where $`\xi `$ (the matter density enhancement at $`x=x_s`$) has been defined in Eq. (13) (see also Schi ).

In Eq. (27), the function $`f(x)`$ parametrizes the rarefaction zone (“hot bubble”) above the shock front, characterized by a drop of density over more than one decade in $`x`$, and by an asymptotic density increase for $`xx_s`$. After some trials, we have chosen the following (purely empirical) parametrization for $`f(x)`$, which reproduces the main features of the graphical profiles in Schi :

$$\mathrm{ln}f(x)=\left[0.280.69\mathrm{ln}(x_s/\mathrm{km})\right]\left[\mathrm{arcsin}(1x/x_s)\right]^{1.1},$$
(28)
In the above equations we assume (following Schi ) a slightly accelerating shock-front position $`x_s`$,

$$x_s(t)=x_s^0+v_st+\frac{1}{2}a_st^2,$$
(29)
with parameters approximately given by

$`x_s^0`$ $``$ $`4.6\times 10^3\mathrm{km},`$ (30)
$`v_s`$ $``$ $`11.3\times 10^3\mathrm{km}/\mathrm{s},`$ (31)
$`a_s`$ $``$ $`0.2\times 10^3\mathrm{km}/\mathrm{s}^2.`$ (32)
Notice that: ($`i`$) The time dependence of $`\rho (x)`$ in Eq. (27) is implicit, through the function $`x_s=x_s(t)`$; and ($`ii`$) at small radii, it is

$$xx_s\rho (x)\rho _0(x)\times 16(x_s/\mathrm{km})^{1.13}.$$
(33)
Figure 1 shows the neutrino potential $`V(x)`$ as derived from the above parametrization \[Eqs.(26)–(32)\] for representative post-bounce times $`1`$ s, as well as for the static profile $`(t0)`$. These curves reasonably reproduce the main features of the shock profiles from Schi . The horizontal bands represent the ranges where the neutrino potential equals the two wavenumbers ($`k_H=V`$ and $`k_L=V`$), for a representative energy interval $`E[4,\mathrm{\hspace{0.17em}70}]`$ MeV.

As previously mentioned, a line at constant $`k_H`$ can intersect the $`V(x)`$ profile in Fig. 1 in at most three points $`x_i`$, ordered as $`x_1x_2(x_s)x_3`$, and corresponding to level crossing probabilities $`P_1`$, $`P_2`$, and $`P_3`$. A further crossing probability $`P_s`$ is related to the step-like feature of the profile at the shock front $`x=x_s`$. Two of the four critical points ($`x_1,x_2,x_s,x_3`$) can merge in the following cases: (1) at the bottom of the rarefaction zone (denoted as $`x_r`$), where $`x_1=x_2`$; (2) at the shock front $`(x_s)`$, where it is $`x_2=x_s`$ for $`V^+=k_H`$ or $`x_3=x_s`$ for $`V^{}=k_H`$. Single transitions can instead occur at early times or high energies ($`V(x)=k_H`$ at $`x=x_3`$ only), as well as at late times or low energies ($`V(x)=k_H`$ at $`x=x_1`$ only). Therefore, the three critical cases $`V(x_r)=k_H`$ ($`x_1=x_2`$) and $`V^\pm =k_H`$ ($`x_{2,3}=x_s`$) are expected to mark significant changes in the behavior of $`P_H`$.

We conclude this section with a few cautionary remarks. At present, the detailed shape of $`V(x)`$ at the shock front is poorly known, since the physical requirement of a density discontinuity at $`x=x_s`$ is implemented “by hand” Schi , as a necessary correction to the artificially smooth profile from simulations Shoc .777In Ref. Schi , the authors explicitly say to have steepened the shock front because “… in hydrodynamics calculations \[it\] may be softened by numerical techniques.” A side effect of the shock front steepening is the appearance of a local “cusp” at $`x=x_s`$ (as apparent in Ref. Schi  and even more in our Fig. 1), which might well be unphysical. However, we do not smooth it out, since it actually provides us with a useful, nontrivial check of numerical versus analytical calculations of $`P_H`$ in an “extreme” condition (i.e., around a sudden change in the sign and value of the density gradient). Finally, we observe that some “leftover” turbulence might be expected in the rarefaction zone behind the shock front, inducing fuzzy variations of the local density. Such variations might lead to more than three solutions of Eq. (22) in some cases, and thus to a more complicated and “random” structure for $`P_H`$. This possibility (not considered in this work) should be kept in mind when more refined supernova shock simulations will become available.

## IV Analysis of crossing probabilities in the energy domain

In this Section we study in detail the energy dependence of $`P_H`$ at fixed time. We pay particular attention to the discussion of multiple transitions along the shock density profile.

### IV.1 Analytical vs numerical calculations of $`P_H(E)`$

Figure 2 shows our calculation of $`P_H(E)`$ at fixed post-bounce time $`t=4\mathrm{s}`$, for five representative values of $`\mathrm{sin}^2\theta _{13}`$, ranging from $`10^6`$ (top) to $`10^2`$ (bottom). In the left panels, a direct comparison is made between analytical calculations (solid curves) and Runge-Kutta numerical calculations (dots), performed at $`5\times 10^3`$ equally-spaced points in the interval $`E[0,80]`$ MeV. In general, up to four energy regimes can be identified, with significant changes at the three critical energies $`E5`$, 50, and 67 MeV. These energies fulfill (for the $`V(x)`$ profile at $`t=4`$ s) the three critical conditions $`k_H=V^+,V^{},V(x_r)`$ which, as discussed in the previous section, signal the merging of two of the four possible level crossing points.888Note that these critical conditions are independent of $`\theta _{13}`$.

In the left panel of Fig. 2, the results of the numerical calculations appear to be generally scattered around the analytical curve. The reason is that Runge-Kutta calculations keep track of the phase(s) of the level crossing amplitude(s), while our analytical approximations average out this information from the beginning. For a single numerical amplitude, the associated phase becomes ineffective after taking the squared modulus. For multiple (interfering) amplitudes, instead, the relative phases do provide an oscillating structure in $`P_H`$, which is apparent in the numerical results.999In the left panels of Fig. 2, for the sake of graphical clearness, we show the spread of numerical results as scatter plots, rather than through rapidly oscillating curves.

In the right panels of Fig. 2, the oscillating structure of the numerical results is averaged out by a convolution with a Gaussian function (with one-sigma width of $`\pm 1`$ MeV), which simulates a generic “smearing” process (e.g., experimental energy resolution). The same convolution is applied to the analytical results. It can be seen that there is very good agreement between the two different (numerical and analytical) approaches after smearing. The results of Fig. 2 (and of other checks that we have performed for different post-bounce times) show that analytical calculations of $`P_H`$ coincide with phase-averaged numerical results with very good accuracy, even close to the “critical” points where two transitions collapse. Moreover, analytical calculations are actually much more convenient than the numerical ones, which require a very high sampling rate in order to perform efficient phase-averaging and to prevent numerical artifacts.

Let us now discuss in more detail the behavior of the analytical $`P_H`$ in the left panels of Fig. 2. We remind that $`P_s`$ \[the crossing probability at the shock front, Eq. (19)\] has almost a top-hat behavior: it rapidly drops from $`1`$ to zero for $`k_H`$ outside the range $`[V^{},V^+]`$, i.e., for $`E`$ outside the range $`[5,50]`$ MeV (for the specific profile used in Fig. 2). This behavior depends only mildly on $`\mathrm{sin}^2\theta _{13}`$ \[as far as it is $`\stackrel{<}{}O(10^2)`$\]. The LZ probabilities $`P_{1,2,3}`$ are instead exponentially suppressed as $`\mathrm{sin}^2\theta _{13}`$ increases \[see Eq. 15\]. In particular, in the bottom panels of Fig. 2 ($`\mathrm{sin}^2\theta _{13}=10^3`$ and $`10^2`$), $`P_H`$ is dominated by $`P_s`$ (with its characteristic top-hat behavior), with very little contributions from $`P_{1,2,3}`$. The dominance of a single crossing amplitude in $`P_H`$ explains the suppression of the numerical oscillations in the same panels.

Conversely, for $`\mathrm{sin}^2\theta _{13}`$ as small as $`10^6`$ eV2 (top panel in Fig. 2), the probabilities $`P_{1,2,3}`$ (as well as $`P_s`$) can all become as large as $`1`$, rendering $`P_H1`$. In particular, for $`E\stackrel{<}{}5`$ MeV, only $`P_11`$ is active (single crossing, well above the shock front). For $`5\stackrel{<}{}E\stackrel{<}{}50`$ MeV, besides $`P_s1`$ it is also $`P_{1,2}1`$ (two crossings in the rarefaction zone). For $`50\stackrel{<}{}E\stackrel{<}{}67`$ MeV, $`P_{1,2}`$ are still active, $`P_31`$ is switched on, and $`P_s`$ drops to zero. Finally, for $`E67`$ MeV, only $`P_3`$ (level crossing well below the shock front) survives. In all such cases (characterized by either single or triple crossing with strong nonadiabatic character) one derives from Eq. (25) that $`P_H1`$.

The intermediate panels in Fig. 2, corresponding to $`\mathrm{sin}^2\theta _{13}=10^5`$ and $`10^4`$, are not as easily understood as the previous ones, since at least one of the three probabilities $`P_{1,2,3}`$ is definitely $`<1`$. In such cases, $`P_H`$ is better understood by explicitly separating its components $`P_{1,2,s,3}`$, as done in the next subsection.

### IV.2 $`P_H(E)`$ decomposition

Figure 3 shows the decomposition of $`P_H(E)`$ (analytical, unsmeared) into the four crossing probabilities $`P_{1,2,s,3}`$ for $`\mathrm{sin}^2\theta _{13}=10^5`$ (top panel) and $`\mathrm{sin}^2\theta _{13}=10^4`$ (bottom panel). In both cases, the shock wave profile at $`t=4`$ s (see Fig. 1) is assumed. The various transition probabilities are switched on and off at the three critical energies ($`E5`$, 50, and 67 MeV) associated to such profile, as previously discussed.

In the top panel of Fig. 3, both $`P_1`$ and $`P_3`$ are strongly nonadiabatic ($`1`$), and thus $`P_H=P_11`$ and $`P_H=P_31`$ in the single-transition regimes at low and high energy, respectively. At intermediate energies both $`P_s`$ and $`P_2`$ are nonzero, but the nonadiabatic character of $`P_2`$ is less pronounced, since at $`x_2`$ (outer part of the rarefaction zone), the density gradient is smaller \[and the scale factor in Eq. (16) is larger\] than at $`x_1`$ or $`x_3`$—at least in our simplified description of the density profile. Therefore, in the intermediate regime $`5\stackrel{<}{}E\stackrel{<}{}50`$ MeV, the strong nonadiabaticity of $`P_1`$ and $`P_s`$ implies $`P_HP_2`$ through Eq. (25). In summary, in the top panel of Fig. 3, the overall crossing probability $`P_H`$ initially takes the value $`P_1`$, then $`P_2`$, and eventually $`P_3`$, as the energy increases and passes through its critical values. All relevant parts of the density profile (static profile, shock front, and both sides of the rarefaction zone) can thus contribute to the $`P_H`$ energy profile.

In the bottom panel of Fig. 3, the higher value of $`\mathrm{sin}^2\theta _{13}`$ leads to an overall decrease of $`P_{1,2,3}`$, particularly for the second transition ($`P_20`$), while $`P_s`$ remains strongly nonadiabatic. Therefore, the behavior of $`P_H`$ is not obvious, except at very low energy (where only $`P_1>0`$) or at very high energy (where only $`P_3>0`$). We can only say that the smallness of $`P_2`$ implies that the outer part of the rarefaction zone (just behind the shock front) are not relevant in this case.

Summarizing, variations of $`\mathrm{sin}^2\theta _{13}`$ for fixed density profile can modulate significantly the $`P_{1,2,3}`$ contributions to the crossing probability $`P_H`$. The contribution of $`P_s`$ remains instead stable and strongly nonadiabatic. As a consequence, different parts of the neutrino profile can play leading or subleading roles, according to the chosen value of $`\theta _{13}`$. When $`\theta _{13}`$ will be known, one shall be able to gauge the relative importance of such contributions, and to determine which parts of the shock profile should be most intensively studied. In the meantime, all parts of the shock profile should be considered as potentially important and worth further study.

### IV.3 $`P_H(E)`$ at different times, and comments on $`P_L`$

We conclude the analysis in the energy domain by briefly discussing the variations of $`P_H(E)`$ for different post-bounce times $`t`$. We also comment on the contribution of the shock-front discontinuity to $`P_L`$.

Figure 4 shows the function $`P_H(E)`$ at three different times ($`t=0`$, 4, and 8 s) for $`\mathrm{sin}^2\theta _{13}=10^5`$ (top panel) and $`10^4`$ (bottom panel). The curves at $`t=4\mathrm{s}`$ have already been discussed. At earlier times $`t\stackrel{<}{}1\mathrm{s}`$ ($`t=0\mathrm{s}`$ in Fig. 4) the behavior of $`P_H`$ becomes simpler, being basically dictated by a single LZ crossing along the static profile. At later times $`(t=8\mathrm{s})`$, the critical energies at which $`k_H=V^\pm `$ and $`k_H=V(x_r)`$ move forward and partly fall outside the energy range relevant for supernova neutrinos. As a consequence, the nonadiabatic character of the innermost transition at $`x_1`$ makes $`P_H1`$ in almost the first half of the energy range of Fig. 4, with a rapid drop to smaller values as the more adiabatic transition $`P_2`$ becomes relevant. The drift of the zones where $`P_H`$ is high or low has thus an interesting time structure, which is analyzed in more detail in the next section.

We conclude this section with a few comments on the crossing probability $`P_L`$ in the $`L`$ subsystem. As mentioned in Sec. II C, the step-like nature of the shock front profile can make $`P_LP_s0`$ for $`k_L`$ roughly within the range $`[V^{},V^+]`$. In practice, however, this condition occurs only in the low-energy tail of the (anti)neutrino spectrum $`(E\stackrel{<}{}\mathrm{few}\mathrm{MeV})`$ and at relatively late times $`(t\stackrel{>}{}10\mathrm{s})`$. Low-energy effects are suppressed by the lower cross section and by experimental detection thresholds, while late-time effects are intrinsically suppressed by the exponential decrease of the supernova neutrino luminosity on a timescale of a few seconds (see also Sec. VI). Moreover, the $`L`$-transition is never strongly nonadiabatic; in fact, from Eq. (19) one can prove that $`P_L(\nu )`$ and $`P_L(\overline{\nu })`$ are always significantly smaller than $`\mathrm{cos}^2\theta _{12}`$ and $`\mathrm{sin}^2\theta _{12}`$, respectively. For such reasons, effects related to $`P_L0`$ at the shock front are hardly observable, and have been neglected throughout this work.

## V Analysis of crossing probabilities in the time domain

In this Section we study in detail the time dependence of $`P_H`$ at fixed energy. As in Sec. IV, we start with a comparison between numerical and analytical calculations of $`P_H`$, and then we study the decomposition of $`P_H`$ in terms of $`P_{1,2,s,3}`$.

### V.1 Analytical vs numerical calculations of $`P_H(t)`$

Figure 5 shows our calculation of $`P_H(t)`$ at fixed neutrino energy $`E=30`$ MeV and for five representative values of $`\mathrm{sin}^2\theta _{13}`$, ranging from $`10^6`$ (top) to $`10^2`$ (bottom). In the left panels, a direct comparison is made between analytical calculations (solid curves) and Runge-Kutta numerical calculations (dots), performed at $`5\times 10^3`$ equally-spaced points in the shown time interval. Once again, the interference among multiple $`\nu `$ transition amplitudes appears to generate fast oscillations, which scatter the numerical results around the analytical curves. In the right panels of Fig. 5, phase averaging (smearing) is instead enforced through a convolution with a “top-hat” time resolution function with $`\pm 0.5\mathrm{s}`$ width. It can be seen that such phase averaging brings the numerical and analytical calculations in very good agreement at any $`t`$, up to small residual artifacts in the smeared numerical results, which would disappear by sampling the abscissa more densely (not shown). From the results of Figs. 2 and 5 we conclude that simple, analytical calculations provide the correct phase-averaged value of $`P_H`$ as a function of both $`E`$ and $`t`$. Therefore, time-consuming numerical estimates of neutrino transitions along the shock profile Schi ; Shoc  can be accurately replaced by much faster and elementary calculations. This is one of the main results of our work.

Let us now discuss in more detail the behavior of the analytical $`P_H`$ in the left panels of Fig. 5. In analogy with the previous discussion in the energy domain, the structure of $`P_H`$ appears to be characterized by (up to) four different regimes also in the time domain. These regimes are separated by three “critical times,” namely, $`t3.2\mathrm{s}`$ (when $`k_H`$ equals $`V`$ at the bottom of the rarefaction zone $`x_r`$), $`t3.3\mathrm{s}`$ (when $`k_H`$ equals $`V^{}`$ at the bottom of the shock front), and $`t7.7\mathrm{s}`$ (when $`k_H`$ equals $`V^+`$ at the top of the shock front). \[Of course, for $`E30`$ MeV such critical times will be different.\] The probabilities $`P_{1,2,s,3}`$ play different roles in forming the global $`P_H(t)`$ in these different regimes. For $`t\stackrel{<}{}3.2`$ s, the relevant density profile is essentially static, and $`P_HP_3\mathrm{const}^{\mathrm{sin}^2\theta _{13}}`$. Similarly, for $`t\stackrel{>}{}7.7\mathrm{s}`$, it is $`P_HP_1\mathrm{const}^{\mathrm{sin}^2\theta _{13}}`$. The single-transition values of $`P_H`$ at early and late times are thus reduced by successive powers of ten as $`\mathrm{sin}^2\theta _{13}`$ is reduced by factors of ten from top to bottom in Fig. 5. At intermediate times the situation is instead less obvious. For $`t3.2`$–3.3 s, the crossing probabilities $`P_2`$ and $`P_s`$ are switched on in succession, leading to a complex local structure in $`P_H(t)`$ (especially at small $`\mathrm{sin}^2\theta _{13}`$); in this situation, neutrino matter effects are probing the lower part of both the rarefaction zone and of the shock front (see Fig. 1). In the interval $`3.3\stackrel{<}{}t\stackrel{<}{}7.7`$ s, $`P_s`$ is always strongly nonadiabatic, $`P_3`$ is negligible, while $`P_{1,2}`$ play a role only at small $`\mathrm{sin}^2\theta _{13}`$ (being exponentially suppressed for increasing $`\theta _{13}`$). For $`\mathrm{sin}^2\theta _{13}=10^2`$ (bottom panels), only $`P_s`$ is important in forming the global $`P_H`$.

The most complicated cases in Fig. 5 thus emerge at “intermediate times” (i.e., when both the rarefaction zone and the shock front are being probed) and at “intermediate $`\theta _{13}`$ mixing” (i.e., when $`P_{1,2,3}`$ are moderately nonadiabatic). In such cases, a decomposition of $`P_H`$ into $`P_{1,2,s,3}`$ is necessary to grasp the main features of $`P_H(t)`$.

### V.2 $`P_H(t)`$ decomposition

Figure 6 shows the decomposition of $`P_H(t)`$ (analytical, unsmeared) into the four crossing probabilities $`P_{1,2,s,3}`$ at $`E=30`$ MeV, for $`\mathrm{sin}^2\theta _{13}=10^5`$ (top panel) and $`\mathrm{sin}^2\theta _{13}=10^4`$ (bottom panel). The various transition probabilities are switched on and off at the three critical times $`t3.2`$, 3.3, and 7.7 s. In both panels, $`P_3`$ and $`P_1`$ dominate at early and late times, respectively. At intermediate times, the structure of $`P_H`$ depends sensitively \[through Eq. (25)\] upon the specific values of the crossing probabilities $`(P_{1,2})`$ in the rarefaction zone ($`P_s`$ being almost constant and close to $`1`$). The results in this figure confirm that variations of $`\theta _{13}`$ can lead to dramatic variations in the role of each $`P_i`$ (and thus of the different parts of the shock profile) in building the final crossing probability $`P_H`$ as a function of both energy and time.

## VI Impact on observable positron spectra

As an application of the analytical calculation of $`P_H`$ described in the previous section, we study the effect of the shock propagation on the energy (and time) spectra of positrons detectable at the Earth through the inverse beta-decay reaction

$$\overline{\nu }_e+pn+e^+.$$
(34)
We assume a “standard” supernova explosion at $`D=10`$ kpc, releasing a total energy $`E_B=3\times 10^{53}`$ erg, equally shared among the (anti)neutrino flavors $`\nu _\alpha `$, and distributed in time as

$$L(t)=\frac{E_B}{6\tau }e^{t/\tau },$$
(35)
where $`L`$ is the luminosity for each $`\nu _\alpha `$ flux, and the decay time is taken as $`\tau =3\mathrm{s}`$ Schi . For the sake of simplicity, we assume unpinched (normalized) Fermi-Dirac $`\nu _\alpha `$ spectra with time-independent temperatures $`T_\alpha `$,

$$f_\alpha (E)=\frac{0.176}{T_\alpha ^4}\frac{L(t)}{4\pi D^2}\frac{E^2}{e^{E/T_\alpha }+1},$$
(36)
where we take $`T_{\overline{e}}=4.5`$ MeV and $`T_{\overline{x}}=6.5`$ MeV ($`x=\mu ,\tau `$). 101010For recent discussions of supernova neutrino energy spectra, see Unce ; Flux .

In the presence of oscillations, the supernova $`\overline{\nu }_e`$ spectrum is given by

$$f(E)=f_{\overline{e}}(E)P_{ee}+f_{\overline{x}}(E)(1P_{ee}),$$
(37)
where $`P_{ee}`$ is taken from Eq. (21) (antineutrino case). The convolution of $`f(E)`$ with the differential cross section $`d\sigma (E,E_{\mathrm{pos}})/dE_{\mathrm{pos}}`$ Stru  provides the event spectrum $`dN/dE_{\mathrm{pos}}`$ in terms of the positron energy $`E_{\mathrm{pos}}`$. We assume perfect energy resolution and zero threshold, in order to show the shock effects in the most favorable conditions. It is understood that such effects will be somewhat degraded in realistic cases, depending on detector details. Finally, just to fix the overall scale, we assume a detector volume corresponding to 32 kton of water. However, for the main purposes of our discussion, the event rate units could be taken as arbitrary in the following figures.

Figure 7 shows our calculated positron spectra for no oscillations ($`P_{ee}=1`$, dashed curve) and hypothetical full conversion ($`P_{ee}=0`$, dot-dashed curve). For partial conversion, Eq. (37) implies that the final positron spectrum is basically a linear combination of the previous two spectra. Examples are given in Fig. 7 for $`\mathrm{sin}^2\theta _{13}=10^5`$, for both normal hierarchy (dotted curve) and inverted hierarchy (solid curve). For normal hierarchy, the coefficients of the combination are energy-independent, namely, $`P_{ee}=\mathrm{cos}^2\theta _{12}`$ and $`1P_{ee}=\mathrm{sin}^2\theta _{12}`$ \[see Eq. 21\]. Conversely, for inverse hierarchy, such coefficients become energy- and time-dependent through $`P_H`$ \[$`P_{ee}=\mathrm{cos}^2\theta _{12}P_H(E)`$\]. In the specific example of Fig. 7, the time is fixed through integration over a representative bin \[$`t=6\pm 0.25\mathrm{s}`$\]. The energy profile of $`P_H(E)`$ then governs the relative weight of the spectra $`f_{\overline{e}}`$ and $`f_{\overline{x}}`$ in the combination. As $`P_H`$ decreases from its low-energy value at the first critical energy (see Figs. 3 and 4 and related comments), the weight of $`f_{\overline{x}}`$ increases and the positron spectrum is shifted to higher energies. The “shoulder” around the peak of the solid curve is thus an imprint of the shock passage. This shoulder is expected to drift for time bins different from $`t=6\pm 0.25\mathrm{s}`$.

Figure 8 shows a set of successive positron energy spectra, binned in time intervals of 0.5 s, from about 6 to 9 seconds after bounce. The overall decrease of the spectra for increasing $`t`$ is due to the decrease in luminosity. The top panel does not include Earth effects, while the bottom panels includes representative Earth mantle crossing effects, calculated in the same conditions as in Ours . As in Fig. 7, the cases of normal and inverted hierarchy are distinguished through dotted and solid curves, respectively, and $`\mathrm{sin}^2\theta _{13}`$ is fixed at $`10^5`$. The top panel shows how the shock passage can leave an observable, time-dependent spectral deformation, in the case of inverse hierarchy. This possibility is very exciting, although one might realistically hope to see no more than variations in the first spectral moments Shoc ; Luna .

From the comparison of the top and bottom panel in Fig. 8, it turns out that the spectral modulation due to the shock might be partly hindered, at late times, by the additional modulation generated by Earth effects. In other words, the “wiggles” produced by oscillations in the Earth might make it difficult to identify the position of the genuine, shock-dependent spectral deformations at different times. Therefore, in some cases Earth effects may not necessarily represent an additional handle Luna  to study the shock dynamics via neutrinos.

Finally, Fig. 9 shows the drop in the positron event rate (arbitrary units) as a function of time and for four representative values of $`\mathrm{sin}^2\theta _{13}`$, ranging from $`10^5`$ to $`10^2`$. Four rates are displayed in each panel, integrated over consecutive 10 MeV energy bins. For normal hierarchy (dotted curves), the rates decrease according to the assumed exponential law in Eq. (35). For inverted hierarchy (solid curves), the characteristic time structure of $`P_H(t)`$ (see Fig. 5) is instead reflected in a characteristic deviation of the rate decrease with respect to a pure exponential drop: The lower $`P_H(t)`$, the higher the rate (enriched in $`\overline{\nu }_x`$ events) for inverse hierarchy, as compared with normal hierarchy. The two hierarchies become instead indistinguishable when $`P_H(t)1`$ (see Fig. 5), since $`P_{ee}\mathrm{cos}^2\theta _{12}`$ for both normal and inverted neutrino mass spectra \[Eq. (21)\]. Roughly speaking, the time structure of the event rates in Fig. 9 reflects (although “up-side-down”) the $`P_H(t)`$ pattern in Fig. 5.

Shock-wave signatures in the time domain have been emphasized in Schi . It goes without saying that very high statistics would be needed to discriminate such time structures in real experiments.111111A more complete and quantitative study of the observability of shock-induced structures in supernova (anti)neutrino signals will be performed elsewhere, in the context of future large-volume, high-statistics detectors. Nevertheless, the reward could be very high. Let us consider, e.g., the “bathtub” pattern of the solid curves in the lower right panel of Fig. 9 ($`\mathrm{sin}^2\theta _{13}=10^2`$). The (hypothetical) experimental detection of such rather distinctive pattern would, at the same time: ($`i`$) provide a snapshot (if not a “movie”) of the shock wave propagation; ($`ii`$) prove that the neutrino mass hierarchy is inverted; and ($`iii`$) put a significant lower bound on $`\mathrm{sin}^2\theta _{13}`$. Each of these results would have a dramatic impact on our understanding of both supernova and neutrino properties.

## VII Summary and prospects

Supernova shock propagation can produce observable effects in the energy and time structure of the neutrino signal, as suggested in Schi  and also investigated in Shoc ; Luna . In the calculation of such effects for $`3\nu `$ oscillations, the neutrino crossing probability $`P_H`$ in the 1-3 neutrino subsystem plays an important role. The evaluation of $`P_H`$ can be done either numerically or analytically, for given shock profiles in space and time. By assuming a simplified description of the time-dependent density profile, we have explicitly shown that simple analytical calculations accurately reproduce the phase-averaged numerical results, also in the nonobvious case of multiple transitions along the shock profile. The analytical approach has then been used to study some relevant characteristics of the structure of $`P_H`$ in energy and time. This structure can in part be reflected in observable signals, as we have shown through some simple and selected examples.

The simplifications used in this work, and especially the empirical parameterization of the density profile, do not alter our main conclusions about the validity and usefulness of a simple analytical approach to $`P_H`$, as compared with brute-force numerical calculations. However, such simplifications may affect the specific values of $`P_H`$ that we have estimated at given energy and time. These values (and the corresponding observable effects) should not be taken too literally, since they depend on features of the shock profile which are beyond our control. In particular, the shape of the rarefaction zone left behind by the shock wave should be constrained, if possible, by dedicated supernova simulations with higher resolution in space and time. When such issues will be clarified, more realistic calculations of $`P_H`$, and of its imprint on observable neutrino signals at the Earth, will become possible.

###### Acknowledgements.

This work was in part supported by the Italian Ministero dell’Istruzione, Università e Ricerca (MIUR) and Istituto Nazionale di Fisica Nucleare (INFN) within the “Astroparticle Physics” research project. We thank C. Lunardini and T. Janka for very useful comments. One of us (E.L.) acknowledges kind hospitality at the Institute for Advanced Study (Princeton, New Jersey) where this work was completed.