7.6. Linear Stability Analysis¶

7.6.3.1. Linearize the EoM¶

To linear the EoM we start from a state where almost all neutrinos are in one flavor,

Suppose we have a Hamiltonian in flavor basis of the form

the commutator of Hamiltonian and density matrix is

(* Define the density matrices and Hamiltonians *)
rho={
{1,epsilon},
{Conjugate@epsilon,0}
};
rhodiag={
{1/2,epsilon},
{Conjugate@epsilon,-1/2}
};
hamil={
{h1,h},
{Conjugate@h,h2}
};
hamil2={
{-h1,h},
{Conjugate@h,h2}
};

(* Calculate the Commutators *)

hamil.rho-rho.hamil//Simplify

{{h Conjugate[epsilon] - epsilon Conjugate[h], -h + epsilon (h1 - h2)}, {(-h1 + h2) Conjugate[epsilon] + Conjugate[h], -h Conjugate[epsilon] + epsilon Conjugate[h]}}

hamil2.rho-rho.hamil2//Simplify

{{h Conjugate[epsilon] - epsilon Conjugate[h], -h - epsilon (h1 + h2)}, {(h1 + h2) Conjugate[epsilon] + Conjugate[h], -h Conjugate[epsilon] + epsilon Conjugate[h]}}

hamil.rhodiag-rhodiag.hamil//Simplify

{{h Conjugate[epsilon] - epsilon Conjugate[h], -h +epsilon (h1 - h2)}, {(-h1 + h2) Conjugate[epsilon] + Conjugate[h], -h Conjugate[epsilon] + epsilon Conjugate[h]}}

We linearize the equation by keeping only the first order terms of \(\delta\). For this purpose, we need to calculate the neutrino self-interaction \(H_{\nu\nu}\).

However, from the general form of \(H_{\nu\nu}\), which is an integral or convolution of \(\rho\), we would expect that the off diagonal element of the Hamiltonian \(h\), is of first order, if we start from a density matrix that has first order, which is what we do. Thus we expect \(h \delta^*\) is second order effect, which we will neglect.

Finally, we obtain one equation for each beam, which can either be the (1,2) element or the (2,1) element.

7.6.3.2. Four-Beam Line Model¶

For any line model of finite beams, we can specify each beam by three parameters,

where \(\rho\) is the density matrix of the beam, \(\theta\) is the angle of the beam defined by some rule, \(\alpha\) is the ratio of the particle number density to the neutrino number density. If we are talking about a neutrino beam instead of an anti-neutrino beam, we have \(\alpha=1\).

In the four-beam case, we define the system using the following three lists of parameters,

where the \(\epsilon\)’s are used to construct the perturbed density matrix,

So we can now write down the equation of motion for the system with this perturbed density.

With all the definitions and conventions specified, we can write down the equation of motion without trouble, IN PRINCIPLE.

First of all, we find the Hamiltonian,

The neutrino self interaction term requires some elabration on it. We take the left neutrino beam as an example. It experiences interactions with three beams, \(\{\bar\rho^L, \bar\theta^L, \alpha\}\), \(\{\bar\rho^R, \bar\theta^R, \alpha\}\), as well as \(\{ \rho^R, \theta^R, 1\}\). So \(H_{\nu\nu}^L\) should have three terms,

This procedure works for all other beams. Or we can use the power of the ALMIGHTY Mathematica.

The equation of motion is reduced to one equation about \(\delta\)’s for each beam.

where \(h_1\), \(h_2\) are determined by both \(H_v\), \(H_m\), and the self-interaction term \(H_{\nu\nu}\). \(h\) is determined by the expression of \(H_{\nu\nu}\). Then we rewrite the equation into the form

where \(M\) is the coefficient matrix that generates the equations we previously derived. This procedure can be done by Mathematica.

Notice that we have the equation with r as the variable, which is not very convenient. Even we solve the equation, it is very hard to interpretate the solutions since r is different at the same height z. So we have to rewrite the equation into one with vertical height z as the variable using \(i\partial_r = i \sin \theta \partial_z + i \cos \theta \partial_x\). Be very careful with the sign of \(+ i \cos \theta \partial_x\). In the four beam case, we have

The equation for the perturbations becomes

transMatrix[eta_, mk_, mu_, alpha_, theta1_, theta2_, lambda_, omegav_: 0] := Module[{},

   (* theta1 for neutrinos, theta2 for antineutrinos *)

   {
   {lambda + mu - 2 alpha mu - eta omegav + mu Cos[2 theta1] +
       2 alpha mu Sin[theta1] Sin[theta2],
      alpha mu - alpha mu Cos[theta1 - theta2], -mu -
       mu Cos[2 theta1], alpha mu + alpha mu Cos[theta1 + theta2]}/
     Sin[theta1] + mk Cot[theta1] {1, 0, 0, 0},
   {-mu + mu Cos[theta1] Cos[theta2] + mu Sin[theta1] Sin[theta2],
      lambda + 2 mu - alpha mu + eta omegav -
       alpha mu Cos[2 theta2] - 2 mu Sin[theta1] Sin[theta2], -mu -
       mu Cos[theta1] Cos[theta2] + mu Sin[theta1] Sin[theta2],
      alpha mu + alpha mu Cos[2 theta2]}/Sin[theta2] +
   mk Cot[theta2] {0, 1, 0, 0},
   {-mu - mu Cos[2 theta1],
      alpha mu + alpha mu Cos[theta1] Cos[theta2] -
       alpha mu Sin[theta1] Sin[theta2],
      lambda + mu - 2 alpha mu - eta omegav + mu Cos[2 theta1] +
       2 alpha mu Sin[theta1] Sin[theta2],
      alpha mu - alpha mu Cos[theta1] Cos[theta2] -
       alpha mu Sin[theta1] Sin[theta2]}/Sin[theta1] -
   mk Cot[theta1] {0, 0, 1, 0},
   {-mu - mu Cos[theta1] Cos[theta2] + mu Sin[theta1] Sin[theta2],
      alpha mu + alpha mu Cos[2 theta2], -mu +
       mu Cos[theta1] Cos[theta2] + mu Sin[theta1] Sin[theta2],
      lambda + 2 mu - alpha mu + eta omegav -
       alpha mu Cos[2 theta2] - 2 mu Sin[theta1] Sin[theta2]}/
     Sin[theta2] - mk Cot[theta2] {0, 0, 0, 1}
   }
]

If we are using a model that is homogeneous in x direction, the derivative is gone. We assume the solution is of the form \(\epsilon = \epsilon_0 e^{i\Omega z}\). By put the assumption back into the equation we obtain

Linear stability analysis basically becomes finding the eigenvalues of matrix \(M\). A negative imaginary part in \(\Omega\) means the solution can grow exponentially.

For this four-beam model, the eigenvalues can be found analytically by Mathematica, eventhough the solution is a bit tedious. We work out the example using unit of \(\omega_v\), i.e., \(\hat \lambda=\lambda/\omega_v\) and \(\hat\mu = \mu/\omega_v\).

7.6.3.3. Neutrino Line Model with Fourier Analysis¶

Neutrino line model is discussed in [duan2015]. We’ll follow the paper.

The equation of motion is

where Hamiltonian

We discuss the equilibrium case so that the time dependent part vanishes.

For the line model, we have only two directions \(x\) and \(z\), thus the density matrix depends on these two directions, i.e., \(\rho(x,z)\). Since all the neutrinos emitted from a line located at \(z=0\), we can Fourier decompose the density matrix \(\rho(x,z)\) in the x direction

where \(k_0\) is determined by the size of the line \(L\),

Equivalently, we can linearize the equation first then Fourier transform the perturbations. Both methods works.

7.6.3.3.2. Fourier Transform Perturbations¶

The other idea is to linearize the equations first then Fourier transform only the perturbations. The result for the equations of perturbations can be obtained directly from Eq. (7.6),

\[i v_z \partial_z \rho^i(z) = \left[ -\beta^i \eta \omega_v \sigma_3/2, \rho^i(z) \right] + \left[ \lambda \sigma_3/2, \rho^i(z) \right] + \left[\sum_j \mu \alpha^j (1 - \hat v^{i}\cdot \hat v^{j}) \beta^{j}\rho^j(z),\rho^i(z) \right].\]

We linearize the equation first before Fourier decomposition is applied. The linearized equation is basically (7.4), with different notations. Then we Fourier transform the equations,

\[\begin{split}iv_z \partial_z \epsilon^L_m =& \left( m k_0 v_x + (\lambda - \eta \omega_v) + (1 - \alpha) \chi \right) \epsilon^L_m - \chi \epsilon^R_m + \alpha \chi \bar\epsilon^R_m\\ iv_z \partial_z \bar\epsilon^L_m =& \left( m k_0 v_x + (\lambda + \eta \omega_v) + (1-\alpha)\chi \right) \bar\epsilon^L_m - \chi \epsilon^R_m + \alpha \chi \bar\epsilon^R\\ i v_z \partial_z \epsilon^R_m =& \left( -m k_0 v_x + (\lambda - \eta \omega_v ) + (1-\alpha)\chi \right) \epsilon^R_m - \chi \epsilon^L_m + \alpha \chi \bar \epsilon^L_m \\ iv_z \partial_z \bar\epsilon^R_m =& \left( -m k_0 v_x + (\lambda + \eta \omega_v) +(1-\alpha)\chi \right) \bar\epsilon^R_m - \chi \epsilon^L + \alpha \chi \bar\epsilon^L_m.\end{split}\]

The reason we have no coupling between different modes is that linearized equation of motion is linear to all the perturbations.

Construct a vector

(7.7)¶\[\begin{split}\begin{pmatrix} \epsilon^L_m \\ \bar\epsilon^L_m\\ \epsilon^R_m \\ \bar\epsilon^R_m \end{pmatrix},\end{split}\]

from which we develop the matrix equation

\[\begin{split}i\partial_z\begin{pmatrix} \epsilon^L_m \\ \bar\epsilon^L_m\\ \epsilon^R_m \\ \bar\epsilon^R_m \end{pmatrix} = \frac{1}{v_z} \begin{pmatrix} m k_0 v_x + (\lambda - \eta \omega_v) + (1 - \alpha) \chi & 0 & -\chi & \alpha \chi \\ 0 & m k_0 v_x + (\lambda + \eta \omega_v) + (1-\alpha)\chi & -\chi & \alpha \chi \\ -\chi & \alpha \chi & -m k_0 v_x + (\lambda - \eta \omega_v ) + (1-\alpha)\chi & 0 \\ -\chi & \alpha \chi & 0 & -m k_0 v_x + (\lambda + \eta \omega_v) +(1-\alpha)\chi \end{pmatrix}\begin{pmatrix} \epsilon^L_m \\ \bar\epsilon^L_m\\ \epsilon^R_m \\ \bar\epsilon^R_m \end{pmatrix}.\end{split}\]

We define

\[\begin{split}\Upsilon_m = \frac{1}{v_z} \begin{pmatrix} m k_0 v_x + (\lambda - \eta \omega_v) + (1 - \alpha) \chi & 0 & -\chi & \alpha \chi \\ 0 & m k_0 v_x + (\lambda + \eta \omega_v ) + (1-\alpha)\chi & -\chi & \alpha \chi \\ -\chi & \alpha \chi & -m k_0 v_x + (\lambda - \eta \omega_v ) + (1-\alpha)\chi & 0 \\ -\chi & \alpha \chi & 0 & -m k_0 v_x + (\lambda + \eta \omega_v) +(1-\alpha)\chi \end{pmatrix}.\end{split}\]

This non-Hermitian ‘evolution’ matrix \(\Upsilon_m\) introduces many new features in the evolutions of the perturbations since the eigenvalues of it are not garanteed to be real. Any imaginary part of the eigenvalues of it will give us exponential increase.

Plus-Minus Modes

In the paper [duan2015] the authors introduced the definition

\[\begin{split}D^{\pm}_m =& \frac{1}{2} (\epsilon^L_m \pm \epsilon^R_m) - \frac{\alpha}{2} (\bar\epsilon^L_m\pm \bar\epsilon^R_m)\\ S^{\pm}_m =& \frac{1}{2} (\epsilon^L_m \pm \epsilon^R_m) + \frac{\alpha}{2} (\bar\epsilon^L_m\pm \bar\epsilon^R_m).\end{split}\]

The vector of functions to be solve is

\[\begin{split}\begin{pmatrix} D^+_m\\ S^+_m \\ D^-_m\\ S^-_m \end{pmatrix}.\end{split}\]

This is simply a transformation of the vector we have, i.e., Eq. (7.7). The transformation matrix is

\[\begin{split}\mathcal R=\frac{1}{2}\begin{pmatrix} 1 & -\alpha & 1 & -\alpha \\ 1 & \alpha & 1 & \alpha \\ 1 & -\alpha & -1 & \alpha \\ 1 & \alpha & -1 & -\alpha \end{pmatrix},\end{split}\]

so that

\[\begin{split}\begin{pmatrix} D^+_m\\ S^+_m \\ D^-_m\\ S^-_m \end{pmatrix} = \mathcal R \begin{pmatrix} \epsilon^L_m \\ \bar\epsilon^L_m\\ \epsilon^R_m \\ \bar\epsilon^R_m \end{pmatrix}.\end{split}\]

We can find the corresponding ‘Hamiltonian’ matrix for the new vector by applying

\[\Lambda_m = \mathcal R \Upsilon_m \mathcal R^{-1}.\]

What I get is

\[\begin{split}\Lambda_m = \mathcal R \Upsilon_m \mathcal R^{-1}=\frac{1}{v_z}\begin{pmatrix} \lambda & - \eta \omega_v & m k_0 v_x & 0 \\ - (1+\alpha)\chi -\eta \omega_v & (1-\alpha)\chi + \lambda & 0 & m k_0 v_x \\ m k_0 v_x & 0 & 2(1-\alpha)\chi + \lambda & -\eta\omega_v\\ 0 & m k_0 v_x & (1+\alpha)\chi - \eta \omega_v & (1-\alpha)\chi + \lambda \end{pmatrix},\end{split}\]

which is the same as the form derived in the paper [duan2015].

We can easily find the eigenvalues for the matrix \(Uplison_m\) or \(\Lambda_m\). Any imaginary part (not real part because we have an extra i in the equation) of the eigenvalue will lead to exponential growth of the perturbations.

Eigenvalues of \(\Upsilon_m\) and \(\Lambda_m\) are the same

The reason is that the determinant of the following two matrice are the same,

\[\lvert \Upsilon_m \rvert = \lvert \mathcal R \Upsilon_m \mathcal R^{-1} \rvert,\]

since determinant has cyclic permutation symmetric.

We also notice that matter has no effect on this phenomenon because we can remove matter effect by minus a matrix \(\frac{1}{v_z}\lambda \mathbf I\) from matrix \(\Upsilon_m\), while the eigenvalues is not changed. What’s more important, the evolution of the perturbation doesn’t change under such a manipulation.

duan2015(1,2,3)

Duan, H., & Shalgar, S. (2015). Flavor instabilities in the neutrino line model. Physics Letters, Section B: Nuclear, Elementary Particle and High-Energy Physics, 747, 139–143.