1.2.4. Dispersion Relation¶

Reference to Notes

Izaguirre, I., Raffelt, G., & Tamborra, I. (2016). Fast Pairwise Conversion of Supernova Neutrinos: Dispersion-Relation Approach, 21101(January), 1–6. https://doi.org/10.1103/PhysRevLett.118.021101

1.2.4.1. Polarization Tensor¶

In Raffelt’s paper, they defined the polarization tensor as

\[\Pi^{\mu\nu} = \eta^{\mu\nu} + \int \frac{d\Omega}{4\pi} G(\theta,\phi) \frac{v^\mu v^\nu}{k^\mu v_\mu}.\]

For numerical calculations, we lower the second index and multiply on both side \(\omega\),

\[\omega\Pi^\mu_\nu = \omega\delta^\mu_\nu + \int \frac{d\Omega}{4\pi} G(\theta,\phi) \frac{v^\mu v_\nu}{ 1- \frac{k}{\omega} \hat k\cdot \mathbf v },\]

where

\[\begin{split}v^\mu &= \begin{pmatrix} 1 & \sin\theta \cos\phi & \sin\theta \sin\phi & \cos\theta \end{pmatrix}\\ k^\mu &= \begin{pmatrix} \omega & k \sin\theta_k \cos\phi_k & k\sin\theta_k \sin\phi_k & k\cos\theta_k \end{pmatrix}.\end{split}\]

1.2.4.1.1. Parametrization of Polarization Tensor¶

Raffelt et al parametrize \(k^\mu(n)\) where \(n=k/\omega\). Then the polarization tensor is decomposed into two parts,

\[\eta^{\mu\nu}\]

and

\[\frac{1}{\omega}N^{\mu\nu},\]

where

\[N^{\mu\nu} = \int d\Gamma G_{\mathbf v} \frac{v^\mu v^\nu}{ 1-n \hat{\mathbf k} \cdot \mathbf v },\]

with \(k^\mu=(\omega,\mathbf k)\).

Note to self

The actual wave vector that determines the instability is \(K^\mu\) which is related to \(k^\mu\),

\[k^\mu= K^\mu - (\Lambda^\mu + \Phi^\mu).\]

Since \(\Lambda^\mu\) and \(\Phi^\mu\) are real, imaginary part of \(\omega\) (\(\mathbf k\)) equal imaginary part of \(\Omega\) (\(\mathbf K\)). Thus we only discuss the dispersion relation of \(k^\mu\).

Density matrix is written as

\[\begin{split}\rho = \begin{pmatrix} 1 & \epsilon \\ \epsilon^* & -1 \end{pmatrix}.\end{split}\]

The perturbation \(\epsilon\) is assumed to have the form

\[\epsilon = Q(\Omega, \mathbf K) e^{ -i( \Omega t - \mathbf K \cdot \mathbf x ) }.\]

This assumption indicates that even though we find instabilities, a proper initial condition/boundary condition is required to stimulate this instability.

The polarization tensor is in fact

\[\Pi^{\mu\nu} = \eta + \frac{1}{\omega}N^{\mu\nu}.\]

The equation of motion becomes

\[v_\mu \Pi^{\mu\nu} a_\nu =0 ,\]

where

\[a_\nu = - \int d\Gamma v_\nu G_{\mathbf v} Q_{\mathbf v}.\]

Since \(v_\mu\) is (component) of a null dual vector, we require \(\Pi^{\mu\nu} a_\nu\) to be (component) proportional to \(v^\mu\). Since we have a lot of directions, different \(v^\mu\) are independent of each other. So we require \(\Pi^{\mu\nu} a_\nu=0\).

Then we need to find the solution to

\[\mathrm{Det}(\Pi^{\mu\nu})=0,\]

which is simplified to

\[\mathrm{Det}(\omega \eta^{\mu\nu} + N^{\mu\nu}) = 0.\]

We can also use the polarization tensor \(\Pi^\mu_\nu\)

\[\mathrm{Det}(\omega \delta^\mu_\nu + N^{\mu}_\nu) =0\]

is the determinant of a matrix

\[\omega I + \mathbf N^\mu_\nu.\]

Equivalently, we only need to find the eigenvalues of \(\mathbf N^\mu_\nu\) then multiply on each of them by negative sign.

1.2.4.1.2. Examples of Parametrization¶

Raffelt et al proposed that can now solve the dispersion relation by finding the value of \(k^\mu(n)\) for each n. We make the plot \(\omega\) vs \(\mathbf k\).

Here is an example that I calculated.

The axial symmetric system can be calculated easily using this method. The paper gave an example of two polar angle beams with axial symmetry.

../../_images/listpltOmegan1.png — Fig. 1.62 \(\omega(n)\) for \(G=0.5 \delta(\cos\theta- 0.8) + 0.5 \delta(\cos\theta+0.2)\).¶

../../_images/listpltDispersionRelationDecompose1.png — Fig. 1.63 Dispersion relation.¶

We can check what happens for multibeams. I can plot the dispersion relation for similar configuration but with different number of beams.

../../_images/listanimi1.png — Fig. 1.64 Animition of dispersion relation.¶

dataPltNBeamsPlt[Join[Table[1/beams, {n, 1, beams/2}], Table[-1/beams, {n, 1, beams/2}]], Table[Pi/3 + n Pi/2/(beams - 1), {n, 0, beams - 1}], {-10, 10}, 0.049, {{-10, 10}, {-10, 10}}]

I plot the \(\omega(n)\) relation for different number of beams

../../_images/listpltOmegan12List-10.png

Similar to the previous example, confining the range of \(n\) leads to only a partial patch of the dispersion relation.

../../_images/pltDiffBeamsConfined-n-in--1-to-1-beams-10.png — Fig. 1.65 The code for it¶

pltDiffBeamsConfined[beams_] := dataPltNBeamsPlt[ Join[Table[1/beams, {n, 1, beams/2}], Table[-1/beams, {n, 1, beams/2}]], Table[Pi/3 + n Pi/2/(beams - 1), {n, 0, beams - 1}], {-1, 1}, 0.049, {{-10, 10}, {-10, 10}}]

This should be the continuous limit?

As a comparison, we can plot the dispersion relation in a larger range of n for 10 beams.

On the other hand, we can calculate the continuous limit for the same angle range.

../../_images/compare-continuous-and-10-beams-within-n-range--1-to-1.png — Fig. 1.67 Dispersion relation for 10 beams (\(n\in [-1,1]\)), and continuous limit.¶

MEH

1.2.4.2. Analyze the Symmetries in Polarization Tensor¶

Vectors Using Spherical Harmonics

Four velocity can be expressed in terms of spherical harmonics.

\[v^\mu = \sqrt{\pi}\begin{pmatrix} 2 & \sqrt{2/3} (Y_1^{-1} - Y_1^1) & i \sqrt{2/3} (Y_1^{-1} + Y_1^1) & 2\sqrt{1/3} Y_1^0 \end{pmatrix}.\]

In principle, solving the dispersion relation is not easy. Neverthless, symmetries would significantly simplify the problem.

Axial symmetry indicates that the integrals of first orders of \(\sin\phi\cos\phi\), \(\sin\phi\), and \(\cos\phi\) are 0 on the range \(\phi\in [0,2\Pi]\).

We denote the integral

\[\int \frac{d\Omega}{4\pi} G(\theta,\phi) \frac{v^\mu v_\nu}{\omega- k \hat{\mathbf k}\cdot \mathbf v}\]

as \(P^\mu_\nu\). The polarization tensor becomes

\[\Pi^\mu_{\phantom{\mu}\nu} = I + P^\mu_{\phantom{\mu}\nu}.\]

For axial symmetric emission, only terms \(P^0_{\phantom{0}0}, P^0_{\phantom{0}3}, P^3_{\phantom{3}0}, P^3_{\phantom{3}3}, P^1_{\phantom{1}1}, P^2_{\phantom{2}2}\) are nonzero, given \(\mathbf k\) in z direction, i.e., \(\phi_k=\theta_k=0\).

To simplify the calcualtion, we denote \(n=\frac{\lvert \mathbf k\rvert}{\omega}\). We will NOT solve \(\omega(n)\). Instead we write down the form of the eigenvalues of

\[N^\mu_{\phantom{\mu}\nu} = \omega P^\mu_{\phantom{\mu}\nu},\]

which shows is an analytical expression of \(\omega\). We do not solve this relation. Instead, we plugin the definition \(n=\frac{\lvert \mathbf k\rvert}{\omega}\) and find out the relation between \(\omega\) and \(k=\lvert \mathbf k \rvert\)

The four velocity is

\[v^\mu \to ( 1, \sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta )^{\mathrm T}.\]

We define

\[I_n(\theta)=\int_{\cos\theta_2}^{\cos\theta_1} d\cos\theta G(\theta) \frac{\cos^n\theta}{1 - n \cos\theta },\]

where \(\theta_1\) and \(\theta_2\) are

Since

\[\begin{split}\int_0^{2\pi} d\phi&=2\pi \\ \int_0^{2\pi} d\phi \cos^2\phi &=\int_0^{2\phi} d\phi \sin^2\phi = \pi,\end{split}\]

the matrix \(N^\mu_{\phantom{\mu}\nu}\) is simplified,

\[\begin{split}N^\mu_{\phantom{\mu}\nu} = \omega P^\mu_{\phantom{\mu}\nu}\to \begin{pmatrix} \frac{1}{2} I_0 & 0 & 0 & -\frac{1}{2}I_1\\ 0 & -\frac{1}{4}(I_0-I_2) & 0 & 0\\ 0 & 0 & -\frac{1}{4}(I_0-I_2) & 0 \\ \frac{1}{2}I_1 & 0 & 0 & -\frac{1}{2}I_2 \end{pmatrix}.\end{split}\]

The equation we are solving is

\[\operatorname{Det}(\omega + N^{\nu}_{\phantom{\mu}\nu}) = 0.\]

We find the solutions to omega,

(1.35)¶\[\omega = \frac{1}{4}(I_0-I_2), \quad -\frac{1}{4}\left(I_0-I_2\pm \sqrt{ (I_0-2I_1+I_2)(I_0+2I_1+I_2) }\right).\]

Mathematica Code

The Mathematica code for solving \(\omega\) as a function of \(n\) is

Det[{
      {omega + i0/2, 0, 0, -i1/2},
      {0, omega - (i0 - i2)/4, 0, 0},
      {0, 0, omega - (i0 - i2)/4, 0},
      {i1/2, 0, 0, omega - i2/2}
   }]

We plug in the definition \(n=k/\omega\) then solve dispersion relation from each of the solutions in Eq. (1.35).

The questions are

What does each of the solutions mean?

1.2.4.2.1. Eigenvalues and Axial Symmetry¶

By definition, the meaning of polarization tensor, \(\Pi^\mu_\nu a^\nu=0\) implies that \(a^1\) and \(a^2\) are the \(\phi\) angle dependent components. To prove this, we rewrite \(Q\),

\[Q = \frac{a^\mu v_\mu}{k^\nu v_\nu},\]

which clearly shows that the 1, and 2 component of \(a^\mu\) is related to the phi dependence of \(Q\). \(a^1=a^2=0\) indicates that \(Q\) has no \(\phi\) dependence.

Is this related to eigenvectors?

The eigenvalue \(\frac{1}{4}(I_0-I_2)\) of matrix \(N^\mu_{\phantom{\mu}\nu}\) corresponds to eigenvectors \((0,0,1,0)\) and \((0,1,0,0)\).

I don’t think it is related to eigenvalues. However, eigenvalues set limit on the actual solution. When we write down the solution to \(a^\mu\), the coefficients are related to each other because we have determinant of coefficient matrix being 0. There are degeneracies.

That is to say, the part

\[\begin{split}\begin{pmatrix} -\frac{1}{4}(I_0-I_2) & 0 \\ 0 & -\frac{1}{4}(I_0-I_2) \end{pmatrix}\end{split}\]

are the only elements that determines the whether we have a \(\phi\) dependence in \(Q\), since this is the only part that needs to be acted on in Gaussian elimination method. It is obvious that we have

\[a^1=a^2=0.\]

In turn, it determines the angle dependence of \(Q\),

\[Q = \frac{a^0 - a^1\sin\theta\cos\phi - a^2\sin\theta\sin\phi - a^3 \cos\theta}{k^\mu v_\mu} = Q_0 + Q_3(\theta).\]

We have no \(\phi\) dependence in \(Q\) if we foce the emission to be axially symmetric.

1.2.4.2.2. Continuous Emission within Angle Range¶

In this case we have to calculate \(I_n\) specifically for the angle range, then plug in the expression \(n=k/\omega\) to find the dispersion relation.

1.2.4.2.3. Discrete Emission Beams¶

For discrete emission \(G(\theta)=\sum_i G_i \delta(\cos\theta-\cos\theta_i)\), we can define new quantities

\[\tilde I_n(\theta)= \sum_i G_i \frac{\cos^n\theta_i}{1 - n \cos\theta_i }.\]

Thus

\[\begin{split}&\omega = \frac{1}{4}(I_0-I_2) \\ \Rightarrow &\omega = \frac{1}{4} \sum_i G_i \frac{1-\cos^2\theta_i}{1 - n \cos\theta_i }.\end{split}\]

For two sets of beams, we have

\[4 = G_1 \frac{1-\cos^2\theta_1}{\omega - k\cos\theta_1 } + G_2 \frac{1-\cos^2\theta_2}{\omega - k \cos\theta_2 },\]

which is a conic section. We have already used \(n=k/\omega\).

Hyperbola

For an quadratic equation [HyperbolaWikipedia]

\[A_{xx} x^2 + 2 A_{xy} xy + A_{yy} y^2 + 2 B_x x + 2 B_y y + C = 0,\]

it is hyperbola if

\[\begin{split}D := \begin{vmatrix} A_{xx} & A_{xy}\\A_{xy} & A_{yy} \end{vmatrix} < 0.\end{split}\]

Center of the hyperbola \((x_c,y_c)\) is

\[\begin{split}x_c &= -\frac{1}{D} \begin{vmatrix} B_x & A_{xy} \\ B_y & A_{yy} \end{vmatrix}\\ y_c &= -\frac 1 D \begin{vmatrix} A_{xx} & B_x \\A_{xy} & B_y \end{vmatrix}.\end{split}\]

Principal axis is tilted away from x axis by angle \(\beta\)

\[\tan 2\beta = \frac{2A_{xy}}{A_{xx}-A_{yy}}.\]

HyperbolaWikipedia: Hyperbola @ Wikipedia

We can prove that this is a hyperbola. Simplify the equation to standard form of conic sections

\[4\cos\theta_1\cos\theta_2 k^2 - 2 \times 2(\cos \theta_1+\cos\theta_2) \omega k +4\omega^2 + 2 \times \frac{ G_1(1-\cos^2\theta_1)\cos\theta_2 + G_2(1-\cos^2\theta_2)\cos\theta_1 }{2} k - 2 \times\frac{ G_1(1-\cos^2\theta_1) + G_2(1-\cos^2\theta_2) }{2} \omega = 0.\]

The \(\omega~k\) plane

We use \(\omega~k\) plane, so that we map \(k\) to \(x\) and \(\omega\) to \(y\).

So the coefficients are defined as

\[\begin{split}A_{kk} &= 4\cos\theta_1\cos\theta_2\\ A_{k\omega} &= -2(\cos\theta_1+\cos\theta_2)\\ A_{\omega\omega} & = 4 \\ B_k & = \frac{ G_1(1-\cos^2\theta_1)\cos\theta_2 + G_2(1-\cos^2\theta_2)\cos\theta_1 }{2} = \frac{\tilde G_1 \cos\theta_2 + \tilde G_2 \cos\theta_1}{2}\\ B_\omega & = - \frac{ G_1(1-\cos^2\theta_1) + G_2(1-\cos^2\theta_2) }{2} = -\frac{\tilde G_1 + \tilde G_2}{2}.\end{split}\]

The condition for it to be hyperbola is \(D<0\), where

\[\begin{split}D = \begin{vmatrix} 4\cos\theta_1\cos\theta_2 & -2(\cos\theta_1+\cos\theta_2) \\ -2(\cos\theta_1+\cos\theta_2) & 4 \end{vmatrix} = - 4 (\cos\theta_1-\cos\theta_2)^2.\end{split}\]

As long as we have different angles, \(D\) is always less than 0. We always have a hyperbola. The center of the hyperbola is

\[\begin{split}k_c &= -\frac{1}{D} \begin{vmatrix} B_k & A_{k\omega} \\ B_\omega & A_{\omega\omega} \end{vmatrix} = \frac{1}{4 (\cos\theta_1-\cos\theta_2)^2} \begin{vmatrix} \frac{ \tilde G_1\cos\theta_2 + \tilde G_2\cos\theta_1 }{2} & -2(\cos\theta_1+\cos\theta_2)\\ - \frac{ \tilde G_1 + \tilde G_2 }{2} & 4 \end{vmatrix} \\ \omega_c &= -\frac{1}{D} \begin{vmatrix} A_{kk} & B_k \\ A_{k\omega} & B_\omega \end{vmatrix} = \frac{1}{4 (\cos\theta_1-\cos\theta_2)^2} \begin{vmatrix} 4\cos\theta_1\cos\theta_2 & \frac{ \tilde G_1\cos\theta_2 + \tilde G_2\cos\theta_1 }{2}\\ -2(\cos\theta_1+\cos\theta_2) & - \frac{ \tilde G_1 + \tilde G_2 }{2} \end{vmatrix} .\end{split}\]

A special case for it is \(G_1=\pm G_2\), however, the expression for the center doesn’t really simplify that much.

We are interested in gaps, so the asymptotic lines are the lines that we are interested in.

First of all, we need to find out the principal axis. The angle between the principal axis and x axis is defined to be \(\beta\),

\[\tan2 \beta = \frac{2A_{k\omega}}{A_{kk}-A_{\omega\omega}} = \frac{\cos\theta_1+\cos\theta_2}{1-\cos\theta_1\cos\theta_2}.\]

Suppose we have angles \(\theta'_i=\arctan ( \cos\theta_i )\),

\[\beta = \frac{\theta'_1+\theta'_2}{2}.\]

This indicates that the angle \(\beta\) is always within range \(\beta \in [\pi/4,\pi/4]\).

The other solutions

For the solutions

\[\omega = -\frac{1}{4}\left(I_0-I_2\pm \sqrt{ (I_0-2I_1+I_2)(I_0+2I_1+I_2) }\right),\]

it becomes much more complicated.

1.2.4.3. Why is it called Polarization Tensor¶

Why is \(a^\mu\) called polarization vector?

We kind of see why \(a^\mu\) is some kind of polarization given the definition

\[a^\mu = -\int d\Gamma v^\mu G(\theta,\phi) Q(\theta,\phi).\]

In some sense it is a weighted average of \(Q\). And \(S = Q e^{-i k^\mu x_\mu}\) is the “field” we are insterested in.

Comparing to electrodynamics, where we have the field \(A^\mu\) which tells us about the polarization,

Polarization in Electrodynamics

We can assume that the four vector field is

\[A^\mu = \epsilon^\mu a e^{-ikx} + {\epsilon^\mu}^* a^* e^{ikx} \qquad \text{with}\quad k_\mu k^\mu =0.\]

We make sense of it by interpretating \(\epsilon^\mu\) as the polarization vector and \(a\) as the amplitude of the field strength. This solution is only for one particular case. We use it as an example because it is simple.

To calculate the electric field, \(E^i\) where \(i=1,2,3\), we apply the definition of it

\[E^i = - F^{0i}.\]

By working it out, we find \(\epsilon^i\), which is the spatial part of \(\epsilon^\mu\), indeed plays a role in the direction of field.

1.2.4.4. Solving N Beams¶

For consistancy check, we parametrize \(\omega(n)\). For MAA solution, we have

\[\omega(n) = \frac{1}{4} \sum_i G_i \frac{1-\cos^n\theta_i}{1 - n \cos\theta_i }.\]

../../_images/DR-omega-n-2-beams-lr.png — Fig. 1.68 \(\omega(n)\) with two angles \(\cos\theta = 0.9,0.3\) and \(G_1 = -G_2 = -2\pi/2\). The grid lines are the values of n for singularity, \(n=1/0.3,1/0.9\).¶

../../_images/DR-omega-k-2-beams-lr.png — Fig. 1.69 Same parameters as above. The orange dashed lines are the singularity lines.¶

Fig. 1.68 shows that the limit of \(n\to 1/0.9,1/0.3\) leads to \(\omega\) infinities. These infinities also indicates that \(\omega (k)\) at large \(\lvert\omega\rvert\) are straight lines with slope \(n=1/0.9,1/0.3\), as shown in Fig. 1.69.

Similar plots are made for 4,6,8 beams.

../../_images/DR-omega-n-2-4-6-8-beams-lr.png — Fig. 1.70 2, 4, 6, 8 beams with equal division of cosine of emission angle. For example 4 beams are emission at \(\cos\theta= 0.9, 0.7, 0.5, 0.3\).¶

../../_images/DR-omega-k-2-4-6-8-beams-lr.png — Fig. 1.71 Same parameters as above.¶

We also calculated the homogeneous emssion.

../../_images/DR-omega-n-2-4-6-8-beams-homogeneous.png — Fig. 1.72 2, 4, 6, 8 beams for homogeneous \(G\).¶

../../_images/DR-omega-k-2-4-6-8-beams-homogeneous.png — Fig. 1.73 Same parameters as above.¶

Gap

Whenever a gap in \(omega\) appears, we might find in the \(\omega(n)\) plot that we have less solutions to \(n\) for some given \(\omega\) values. It seems that this can only happen for ELN spectrum with crossing.

To illustrate this idea, I shaded the region that \(n(\omega)\) has no solution for the case of two emission angle and crossing in Fig. 1.74.

../../_images/DR-omega-k-2-beams-lr-shaded-region.png — Fig. 1.74 The shaded regions are the region that \(n(\omega)\) has no solution. It’s hard to analytically solve \(n(\omega)\). But we can get an idea of this using \(\omega(n)\) plot.¶

This provides an method to determine whether we have a gap in \(\omega\) or \(k\). For gap in \(k\), we plot out \(k(n)\).

../../_images/DR-k-n-2-beams-lr-shaded-region.png — Fig. 1.75 The shaded region are the region where \(n(k)\) has no solution.¶

The behavor of \(\omega(n)\) and \(k(n)\) are qualititatively the same.

../../_images/DR-k-n-2-beams-lr-and-homogeneous.png — Fig. 1.76 \(k(n)\) for non-crossing and crossing. The are qualitatively the same as \(\omega(n)\). That being said, we always have the same upside down U-shape for spectrum with crossing for both \(k(n)\) and \(\omega(n)\), if the region of emission angles \(\theta\) satisfies \(\cos\theta>0\).¶

1.2.4.5. Solving Continuous Emission¶

Suppose neutrinos are emitted within a angle range \([\theta_1,\theta_2]\). Using Mathematica, we find the three important integrals

\[\begin{split}I_0 &= \int_{c_2 }^{c_1 } d x \frac{1}{1-k\cos\theta/\omega} \\ &= \frac{\omega}{k} \ln \left( \frac{\omega-c_2 k}{\omega-c_1 k} \right) \\ I_1 &= \int_{c_2 }^{c_1 } d x \frac{x}{1-k\cos\theta/\omega} \\ & = \frac{\omega}{k} \left( c_2 -c_1 + \frac{ \omega }{k} \ln \left( \frac{\omega-c_2 k}{\omega-c_1 k} \right) \right) \\ I_2 &= \int_{c_2 }^{c_1 } d x \frac{x^2}{1-k\cos\theta/\omega} \\ &= \frac{\omega}{2 k} (c_2 -c_1 ) \left(2 \frac{\omega}{k} + c_1 +c_2 \right) + \left(\frac{\omega}{k}\right)^3 \ln \left( \frac{\omega-c_2 k}{\omega-c_1 k} \right) .\end{split}\]

where \(c_1=\cos\theta_1\) and \(c_2=\cos\theta_2\).

1.2.4.5.1. Homogeneous Emission¶

Assuming \(G=1\), the MAA solution is

\[\omega = \frac{g}{4}(I_0 - I_2),\]

which becomes

(1.36)¶\[\omega=\frac{(c_2-c_1)(1+(c_1+c_2)n/2)}{4n^3} + \frac{ 1-n^2 }{4n^3}\ln\left( \frac{1-c_2 n}{1-c_1 n} \right)\]

Meanwhile we could write down the MZA/bimodal solution in the form of \(\omega(n)\).

For MAA and MZA we can plot \(\omega\) as a function of n.

../../_images/DR-omega-k-direct-continuous-maa-no-crossing.png — Fig. 1.77 \(\omega(n)\) for MAA solution. The spectrum is feeded in as \(\{\{\cos\theta_1,\cos\theta_2\}, G \}=\{\{0.9, 0.3\}, 1\}\).¶

../../_images/DR-omega-k-direct-continuous-mza-no-crossing.png — Fig. 1.78 \(\omega(n)\) for MZA solution. The spectrum is feeded in as \(\{\{\cos\theta_1,\cos\theta_2\}, G \}=\{\{0.9, 0.3\}, 1\}\).¶

On the other hand, we know \(k=n\omega\), so we have parametrized the dispersion relation using a parameter \(n\).

../../_images/DR-omega-k-direct-continuous-dr-no-crossing.png — Fig. 1.79 Dispersion relation out of Eq. (1.36). For large k, the relation becomes proportional. The discontinuties are at location of \(\lim_{n\to \infty}\{n\omega(n),\omega(n)\}\). In this example, \(\{ -0.184653, 0 \}\) (MAA solution) and \(\{ 0.729306, 0 \}\) (MZA solution).¶

The Limits

There are several limits in the dispersion relation.

From the figure of \(\omega\), we notice the singularities at the two ends of the distribution, \(1/c_1\) and \(1/c_2\). At these points, a tiny change of n will cause a significant change in \(\omega\) and \(k\). In fact the relation between them becomes a proportional relation since \(n\) remains almost constant.

Another limit is \(n\to \infty\). Since \(\ln(n)\) increases slower than \(n\), we have

\[\lim_{n\to\infty}\omega(n) =\lim_{n\to\infty} \frac{ I_0-I_2 }{4} = 0.\]

We can calculate \(\lim_{n\to \infty}\{n\omega(n),\omega(n)\}\), for \(c_1=0.9\) and \(c_1=0.3\),

\[\begin{split}\lim_{n\to\infty}\{n\omega(n),\omega(n)\} &= \{-0.184653, 0. \}, &\qquad \text{for MAA solution}\\ \lim_{n\to\infty}\{n\omega(n),\omega(n)\} &= \{0.729306, 0. \}, &\qquad \text{for MZA solution}.\end{split}\]

In general, we have the limits

\[\begin{split}\lim_{n\to\infty}\{n\omega(n),\omega(n)\} &= \{(c_1^2-c_2^2 + 2\ln \left\vert\frac{c_2}{c_1}\right\vert )/8, 0 \}, &\qquad \text{for MAA solution} \\ \lim_{n\to\infty}\{n\omega(n),\omega(n)\} &= \{ ( c_1^2-c_2^2 - 2\ln \left\vert \frac{c_2}{c_1} \right\vert )/4, 0 \}, &\qquad \text{for MZA solution}.\end{split}\]

1.2.4.5.2. Emission with Crossing¶

I have to break each of the integral into two parts. I calculate \(I_0-I_2\) for the first region then add to it the second region. Within a region \([\theta_1,\theta_2]\) and

\[\begin{split}G=\begin{cases} g_1, \qquad & \theta_1<\theta<\theta_0 \\ g_2, \qquad & \theta_0<\theta<\theta_2 \end{cases}.\end{split}\]

In other words, we have a box-like spetrum.

For MAA solution we define a function,

\[RHS(c_1,c_2,g)=\frac{I_0-I_2}{4} = \frac{g}{4} \left[ \left( \frac{1}{n} - \frac{1}{n^3} \right)\ln\left( \frac{1-n c_2}{1-n c_1} \right) - \frac{c_2-c_1}{n} \left( \frac{c_1+c_2}{2} + \frac{1}{n} \right) \right].\]

The dispersion relation is given by

\[\omega = RHS(c_1,c_0,g_1) + RHS(c_0,c_2,g_2).\]

Then we parametrically plot \(\{ n \omega(n), \omega(n)\}\) to get the dispersion relation, for MAA solution. Similarly I can find that of MZA solution.

Limits

Before we do any numerical calculations, we can calculate the limits first.

\[\begin{split}\lim_{n\to\pm\infty}\{n\omega(n),\omega(n)\} &= \{( g_1(c_1^2-c_0^2) + g_2(c_0^2-c_2^2) + 2 g_1 \ln \frac{c_0}{c_1} + 2 g_2\ln \frac{c_2}{c_0} )/8, 0 \}, &\qquad \text{for MAA solution} \\ \lim_{n\to\infty}\{n\omega(n),\omega(n)\} &= \{ \{two very complicated expression, try Mathematica \}, 0 \}, &\qquad \text{for MZA solution}.\end{split}\]

For \(g_1=g_2=1\) these limit match the homogeous result, which they should.

We also have the large k limit which are \(\omega = \frac{1}{c_i}k\).

Mathematically, we also have

\[\begin{split}\lim_{n\to 1/c_1} \{ n\omega(n),\omega(n) \} &= \{ \mathrm{DirectedInfinity}(g_1 (1-c_1^2) ), c_1 \mathrm{DirectedInfinity}(-\mathrm{Sign}(g_1 (-1+c_1^2) ) ) \}, &\qquad \text{for MAA solution} \\ \lim_{n\to 1/c_0} \{ n\omega(n),\omega(n) \} &= \{ \mathrm{DirectedInfinity}( (g_1 - g_2 ) (1-c_0^2) ), c_0 \mathrm{DirectedInfinity}( (g_1-g_2) (-1+c_0^2) ) \}, &\qquad \text{for MAA solution} \\ \lim_{n\to 1/c_2} \{ n\omega(n),\omega(n) \} &= \{ \mathrm{DirectedInfinity}(g_2 (1-c_2^2) ), c_2 \mathrm{DirectedInfinity}( g_2 (-1+c_2^2) ) \}, &\qquad \text{for MAA solution}.\end{split}\]

For simplicity, we choose \(g_2=-g_1=1\).

../../_images/DR-omega-k-direct-continuous-two-regions-dr-crossing.png — Fig. 1.80 Dispersion relation for spectral crossing. The discontinuties are at \(\{ 0.0944205, 0 \}\) (MAA solution) and \(\{ -0.098841, 0 \}\) (MZA solution). The spectrum I use is \(\{ \{\{\cos\theta_1,\cos\theta_0\}, g1\}, \{\{ \cos\theta_0,\cos\theta_2\},g2\}\}=\{\{\{0.9, 0.6\}, 3\}, \{\{0.6, 0.3\}, -3\}\}\).¶

I can also plot the MAA and MZA soltions for \(\omega(n)\).

../../_images/DR-omega-k-direct-continuous-two-regions-maa-crossing.png — Fig. 1.81 \(\omega(n)\) for MAA solution. The vertical grid lines are \(n=1/c_1, 1/c_2\). The spectrum I use is \(\{\{\{0.9, 0.6\}, 3\}, \{\{0.6, 0.3\}, -3\}\}\). At \(n\to \pm\infty\), both of the lines approaches 0.¶

../../_images/DR-omega-k-direct-continuous-two-regions-mza-crossing.png — Fig. 1.82 \(\omega(n)\) for MZA solution. The vertical grid lines are \(n=1/c_1, 1/c_2\). The spectrum I use is \(\{\{\{0.9, 0.6\}, 3\}, \{\{0.6, 0.3\}, -3\}\}\). At \(n\to \pm\infty\), all lines approach 0. At \(n\to 1/0.9\), + solution becomes \(-\infty\) and th - solution becomes a complex number \(-0.639474 - 36.1614 i\). (This should be irrelevant.) At \(n\to 1/0.3\), + solution approaches \(-\infty\) while - solution approaches 0.0445.¶

The reason we have no real values between \(1/0.9\) and \(1/0.3\) is because the argument of the ln function is negative within this regime.

../../_images/DR-omega-k-direct-continuous-two-regions-arg-ln.png — Fig. 1.83 The argument of the ln function. The vertical grid lines are \(n=1/0.3,1/0.6,1/0.9\). Even \(n=1/0.6\) gives us the zero argument, which means the ln function is infinite, I do not think this is some point that we can have a \(\omega,k\) relation.¶

Some Discussions about \(\omega(n)\)

It seems that if we plot \(-\omega(n)\) for MZA we will have some kind of similar plot compared to MAA. Well the singularities are all the same location. If we do that we will have the same \(\lvert\omega\rvert\) for \(n\) very close to the sigularities thus same \(\lvert k\rvert\). Of couse for MAA and MZA, \(\omega\) have different signs.

In other words, the slopes of \(\omega(k)\) will have the same value for MZA and MAA but in different quadrant.

1.2.4.6. General Discussions of Significance of Spectra¶

../../_images/DR-omega-of-n-for-different-spectra-table-0.9-0.6-0.3.png — Fig. 1.84 Function \(\omega(n)\) for different \(g_1,g_2\). The emission was set to \(G=\begin{cases} g_1, \qquad & \theta_1<\theta<\theta_0 \\ g_2, \qquad & \theta_0<\theta<\theta_2 \end{cases}\) where \(\cos\theta_1=0.9\) and \(\cos\theta_2=0.3\) \(\cos\theta_0=0.6\).¶

../../_images/DR-k-of-n-for-different-spectra-table-0.9-0.6-0.3.png — Fig. 1.85 \(k(n)\) for the same parameters.¶

The corresponding dispersion relations are shown in Fig. 1.86

../../_images/DR-omega-k-continuous-table-0.9-0.6-0.3.png — Fig. 1.86 DR for spectra \(\{ \{ \{ 0.9,0.6 \}, g1 \}, \{ \{ 0.6,0.3 \}, g2 \} \}\). Everytime we have \(g_1\) or \(g_2\) becomes 0, the results should be ignored.¶

../../_images/DR-omega-of-n-for-different-spectra-table-0.9-0.4-0.3.png — Fig. 1.87 Function \(\omega(n)\) for different \(g_1,g_2\). The emission was set to \(G=\begin{cases} g_1, \qquad & \theta_1<\theta<\theta_0 \\ g_2, \qquad & \theta_0<\theta<\theta_2 \end{cases}\) where \(\cos\theta_1=0.9\) and \(\cos\theta_2=0.3\) \(\cos\theta_0=0.4\)¶

../../_images/DR-k-of-n-for-different-spectra-table-0.9-0.4-0.3.png — Fig. 1.88 \(k(n)\) for the above case.¶

../../_images/DR-omega-k-continuous-table-0.9-0.4-0.3.png — Fig. 1.89 Dispersion relations for the above plot.¶

Why

It seems that crossing is important to a change in the number of solution to \(n(\omega)\). If we have crossing, then the number of solutions to \(n(\omega)\) will change for different values of \(\omega\). Otherwise, the number of solutions won’t change.

And a change of the number of solutions indicates a possible gap. I need some verification about the relation between such non-explicit gap and instabilities.

However, I also notice that the combination \({-2,-1}\) also has such change of number of solutions. In this case, I can spot that for some branches we have complex \(k\). What I don’t understand is that we always find real solutions to \(k\) for some other branches. Does it mean that only one such complex solutions is sufficient? And what exactly is the requirement?

Can I derive some expression for the \(\ln\) function for a continuous distribution \(G(\theta)\)?

1.2.4.6.1. Do we really need a gap?¶

As long as we have a point on \(\omega(k)\) plane that \(\omega'(k_b)=\infty\) and \(\omega'(k_{b,-})\omega'(k_{b,+})<0\), we might obtain a complex k for some range of it.

Or I can simply consider \(k(n)\) and find whether we have a region of \(k'(n)=0\) and \(k''(n)\neq 0\) which generates different numbers of solutions in different region of \(k\).

Imagary Part in k

For MAA solution, I can try to solve \(n\) given value of \(\omega\). If I obtain complex value of \(n\), then I have complex value of k thus instability.

Didn’t find any weird numbers here.

Fig. 1.86 also indicates that crossing probably change the number of solutions to \(omega\) given different \(n\) thus change the number of solutions to \(k\) given different \(n\). For example, the lower left panel shows that \(\omega\) has 1 solutions within \(n\in [-\infty,1/0.9]\) but has 3 solutions within \(n\in [1/0.3,\infty]\). The lower right panel have 3 solutions whithin \(n\in [-\infty,1/0.9]\) and \(n\in [1/0.3,\infty]\). This might be important when identifying the possible gaps.

1.2.4.7. How to Use the Dispersion Relation¶

I can analyze some points on the \(\omega\sim k\) plane.

../../_images/lsa-DR-963-MZA.png — Fig. 1.90 MZA solution for spectrum \(\{ \{ \{0.9,0.6\}, 3 \}, \{ \{0.6,0.3\}, -3 \} \}\). In principle I should have some real solutions fall onto the dispersion relation lines. Other points solved directely from linear stability analysis should be the unstable regions.¶

../../_images/lsa-963-kimag-omega-MZA-p.png — Fig. 1.91 for MZA+ solution for spectrum \(\{ \{ \{0.9,0.6\}, 3 \}, \{ \{0.6,0.3\}, -3 \} \}\). The points for \(\omega>0\) are not real because these are values returned without convergence. The values on the left are returned with convergence.¶

../../_images/lsa-963-kimag-omega-MZA-m.png — Fig. 1.92 for MZA- solution. The points for \(\omega<0\) are not real because these are values returned without convergence.¶

../../_images/lsa-DR-93-MZA-Omega-Real.png — Fig. 1.93 MZA solution for spectrum \(\{ \{ \{0.9,0.3\}, 1 \} \}\). Solve complex \(k\) for real \(\omega\).¶

../../_images/lsa-DR-93-MZA-K-Real.png — Fig. 1.94 MZA for spectrum \(\{ \{ \{0.9,0.3\}, 1 \} \}\). Solve complex \(\omega\) for real \(k\).¶

1.2.4.8. Discrete Case and Continuous Case¶

Will we have a continuous case if the number of beams is infinite.

For discrete case

(1.37)¶\[I_0 - I_2 = \sum_{i=1}^N G_i \frac{ 1-\cos^2\theta_i }{\omega - k \cos \theta_i} = \sum_i G_i \frac{1-u_i^2}{\omega-k u_i}.\]

The continuous case is

(1.38)¶\[I_0 - I_2 = \int d\cos\theta G(\cos\theta) \frac{1 - \cos^2\theta}{\omega- k \cos\theta} = \int du G(u) \frac{1-u^2}{\omega-k u} .\]

We notice that Eq. (1.37) and Eq. (1.38) are the same when number of beams becomes large.

\(G_i\) is in fact \(G_i = g_i \Delta u_i\), where \(\Delta u_i\) is the range of \(\cos \theta_i\) around \(\theta_i\).

1.2.4.9. Arbitrary Spectrum¶

../../_images/DR-continuous-spectrum-near-realistic.png — Fig. 1.95 Reproducing Fig. 3 of Raffelt’s paper.¶

../../_images/DR-continuous-spectrum-near-realistic-spectrum-created.png — Fig. 1.96 The spectrum I created and used. \(-11.2738 u^5-11.5219 u^4+9.1079 u^3+12.6683 u^2+4.51014 u+2.50958\).¶

The unstable regions can be calculated exactly by setting \(\omega\) to be real.

../../_images/DR-continuous-spectrum-near-realistic-no-crossing-maa.png — Fig. 1.97 MAA¶

../../_images/DR-continuous-spectrum-near-realistic-no-crossing-mzap.png — Fig. 1.98 MZA+¶

../../_images/DR-continuous-spectrum-near-realistic-no-crossing-mzam.png — Fig. 1.99 MZA-¶