# Energy flow and momentum under paraxial regime

R. Mahillo-Isla  Área de Ing. Proc. de Fab., University of Valladolid.
Paseo del Cauce 59, 47011, Valladolid, Spain<EMAIL_ADDRESS>M. J.
González-Morales  Teoría de la Señal y Comunicaciones e Ingeniería Telemática.
University of Valladolid.
Paseo de Belén 15, 47011, Valladolid, Spain

###### Abstract

The paraxial model of propagation is an approximation to the model described
by the d’Alembert equation. It is widely used to describe beam propagation and
near-field diffraction patterns. Therefore, its use in optics and acoustics
engineering is rather general. On the other hand, energetic balance and
momentum in the electromagnetic or acoustic frameworks are well-known and lay
in their own physical context. When dealing with paraxial solutions, these
analyses are not so clear since paraxial propagation is not supported by the
electromagnetic or mechanic theory. The present document establishes the
fundamental energy and momentum analysis for paraxial solutions based on the
classical approach by studying a Lagrangian density associated to the paraxial
equation. Solutions of the paraxial wave equation, such as plane waves,
Green’s function and Gaussian beams, are studied under this scheme.

###### pacs:

02.30.Jr, 42.60.Jf, 42.25.Bs

##  I Introduction

The paraxial model of propagation is widely used not only to describe beam
propagation Siegman (1986), but also diffraction problems under the Fresnel
approximation Grella (1982); Born and Wolf (1980). Gaussian beams and their
modes are its best known solutions and the Fresnel diffraction formula has
been applied to obtain near field diffraction patterns. The success of
paraxial propagation shows that it is a fine approximation since its is used
use several disciplines, more or less near physics, such as material
processing Ion (2005), spectroscopy Radziemski  _et al._ (1987) and medicine
Welch and van Gemert (2011) to cite a few. Nevertheless, up to our knowledge,
the energy flow and momentum have never been treated within this
approximation. As far as the paraxial approximation holds, the use of the
formulae obtained from the full propagation framework could be a fair
procedure for energy and momentum description. Nonetheless, since the
solutions do not fulfill this scheme, the conservation of these quantities is
not obtained. The classical procedure to study the conservation of these
quantities under the d’Alembert propagation scheme is to perform a
reconstruction method on the paraxial solution, that is, to obtain a solution
of the d’Alembert equation from a paraxial solution. But there is not a unique
method to obtain a d’Alembert solution from a paraxial one Mahillo-Isla and
González-Morales (2015, 2016). Therefore, it is worth developing, if possible,
the energy and momentum properties associated to the paraxial operator. The
derived formulae should be used when dealing with solutions of the paraxial
wave equation. Previous related work has focused on the modification of
Maxwell’s equations to obtain an electromagnetic description of paraxial
solutions Mukunda  _et al._ (1985); Simon  _et al._ (1986) where the paraxial
version of the Poynting vector is obtained. Although there is published
material regarding energy invariants in paraxial solutions Vaveliuk  _et al._
(2007), we have not found the treatment of momentum under paraxial regime
neither in the scalar nor in the vectorial frameworks.

The organization of the article is as follows. First, we present the partial
differential equation to be studied. Then we study a Lagrangian associated to
this equation where the expressions of energy flow and momentum are obtained
as well as their respective conservation equations. Based on these
expressions, we also propose a method for evaluating the quality of a paraxial
solution at every point in the space. Some fundamental solutions of the
paraxial wave equation are considered in this article as examples: paraxial
plane waves and paraxial Green’s function are the classical constituents of
any other solution in terms of plane wave spectrum and convolution integrals
respectively. We also deal with the most popular paraxial solution, the
Gaussian beam. The developed framework allows us to gain a better
understanding of the paraxial model of propagation.

##  II Paraxial propagation

To account not only for monochromatic waves, we use as generalized motion
equation a time dependent scheme that under time harmonic regime turns into
the classical paraxial wave equation. Thus, we deal with a wave and not with
the, by far, more usual complex envelope. The paraxial differential equation
is usually obtained by means of the application of the _slowly varying
envelope approximation_ to a time harmonic solution. This leads to the
parabolic equation for the complex envelope $\Psi$

$\nabla_{\perp}^{2}\Psi+2ik\partial_{z}\Psi=0.$ (1)

In this equation the axial direction is $+z$, the time harmonic dependence is
$e^{-i\omega t}$, the symbol $\perp$ denotes operation in the transverse plane
to $z$ and $k=\frac{\omega}{c}$ is the wavenumber. Note that, the complex
envelope is

$\Psi=\tilde{\psi}\exp(-ikz),$ (2)

and thus, the equation that must be fulfilled by the phasor $\tilde{\psi}$ is

$\nabla_{\perp}^{2}\tilde{\psi}+2ik\partial_{z}\tilde{\psi}+2k^{2}\tilde{\psi}=0.$
(3)

Since, in free space, the wavenumber is proportional to the angular frequency
and the angular frequency can be identified as a partial derivative of the
solution $\psi$ with respect to $t$; Eq. (3) comes from the generalized motion
equation,

$\nabla^{2}_{\perp}\psi-\frac{2}{c}\partial_{tz}\psi-\frac{2}{c^{2}}\partial_{tt}\psi=0,$
(4)

where the axial direction is $+z$ . The solutions with time harmonic
dependence will be analyzed later. Although this equation has been obtained
from a particular choice of time harmonic dependence, it does not depend on
this choice. Conversely, the partial differential equations for the phasor and
the complex envelope does depend on the particular choice of the time harmonic
dependence. Note that there are solutions of the full wave equation that are
also solutions of this equation, for instance, the forward propagating
d’Alembert solution:

$\psi=F(\mathbf{r}_{\perp})g(z-ct),\;\textrm{if}\left\\{\begin{array}[]{l}\nabla^{2}_{\perp}F=0\\\
g(u)\in\mathcal{C}^{2}\end{array}\right.$ (5)

This solution is not of our interest since they have been widely studied
within the framework of transmission lines. The most important point is that
such a solution represents a wave that propagates along the positive $z$ axis
and the energy density flow at any point is well defined at any time as well
as the momentum density. Furthermore, since Eq. (5) is solution of the
d’Alembert equation, which is a hyperbolic partial differential equation, it
is suspected that Eq. (4) should be hyperbolic as well. The matrix associated
to Eq. (4) is

$\bm{A}=\left[\begin{array}[]{cccc}1&0&0&0\\\ 0&1&0&0\\\ 0&0&0&-c^{-1}\\\
0&0&-c^{-1}&-2c^{-2}\end{array}\right].$ (6)

The eigenvalues $\lambda_{n}$ of $\bm{A}$ are

$\displaystyle\lambda_{1}$ $\displaystyle=$ $\displaystyle 1;$ (7a)
$\displaystyle\lambda_{2}$ $\displaystyle=$ $\displaystyle 1;$ (7b)
$\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{1+c^{-2}}-1}{c^{2}},\,>0;$ (7c)
$\displaystyle\lambda_{4}$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{1+c^{-2}}+1}{c^{2}}\,<0;$ (7d)

for any $c\in\mathbb{R}$. Therefore, Eq. (4) is a hyperbolic partial
differential equation as expected. The main goal of this paper is to establish
the general energy and momentum balance for any solution of Eq. (4) and,
consequently for Eq. (3) also. For that, we study the Lagrangian associated to
Eq. (4).

##  III Analysis of the Lagrangian density of the paraxial equation

The Lagrangian density associated to Eq. (4) can be written by simple
inspection of the differential equation,

$\mathscr{L}=\frac{1}{c}\partial_{t}\psi\left(\frac{1}{c}\partial_{t}\psi+\partial_{z}\psi\right)-\frac{1}{2}\left[(\partial_{x}\psi)^{2}+(\partial_{y}\psi)^{2}\right]$
(8)

In fact, we can recover Eq. (4) by applying Euler-Lagrange equation

$\nabla_{4}\cdot\left(\sum_{\mu=x,y,z,t}\frac{\partial\mathscr{L}}{\partial(\partial_{\mu}\psi)}\mathbf{e}_{\mu}\right)=\frac{\partial\mathscr{L}}{\partial\psi},$
(9)

where $\nabla_{4}=\mathbf{e}_{t}\partial_{t}+\nabla$ in terms of the
generalized coordinates. Let us calculate the terms:

$\sum_{\mu=x,y,z,t}\frac{\partial\mathscr{L}}{\partial(\partial_{\mu}\psi)}\mathbf{e}_{\mu}=\frac{1}{c}\left(2\frac{\partial_{t}\psi}{c}+\partial_{z}\psi\right)\mathbf{e}_{t}-\nabla_{\perp}\psi+\frac{\partial_{t}\psi}{c}\mathbf{e}_{z},$
(10)

and the left hand side of Eq. (9) produces

$\frac{2}{c^{2}}\partial_{tt}\psi+\frac{1}{c}\partial_{zt}\psi-\nabla_{\perp}^{2}\psi+\frac{1}{c}\partial_{tz}\psi=\frac{\partial\mathscr{L}}{\partial\psi}.$
(11)

By requiring the necessary conditions of regularity on $\psi$,
$\partial_{zt}\psi=\partial_{tz}\psi$,

$\nabla_{\perp}^{2}\psi-\frac{2}{c}\partial_{tz}\psi-\frac{2}{c^{2}}\partial_{tt}\psi=-\frac{\partial\mathscr{L}}{\partial\psi}.$
(12)

From Eq. (8), it follows that

$\frac{\partial\mathscr{L}}{\partial{\psi}}=0,$ (13)

and the expression in Eq. (4) is recovered. Since the Lagrangian does not
depend on the coordinates explicitly, the hypothesis of Noether’s theorem are
fulfilled. Thus, there are four primary conserved quantities. One of them is
the energy and the other three ones are momentum. The conservation of the
energy is found in the “time” row (usually the first row) of the stress-energy
tensor. This row is the four-component vector

$\mathbf{T}_{t}=\partial_{t}\psi\left(\sum_{\mu=x,y,z,t}\frac{\partial\mathscr{L}}{\partial(\partial_{\mu}\psi)}\mathbf{e}_{\mu}\right)-\mathscr{L}\mathbf{e}_{t}.$
(14)

The energy density $U$ is identified as the t-component of $\mathbf{T}_{t}$,

$U=\frac{(\partial_{t}\psi)^{2}}{c^{2}}+\frac{1}{2}\left[(\partial_{x}\psi)^{2}+(\partial_{y}\psi)^{2}\right].$
(15)

The rest of the components of $\mathbf{T}_{t}$ forms the energy density flow
$\mathbf{S}$,

$\mathbf{T}_{t}=\mathbf{S}+U\mathbf{e}_{t}.$ (16)

Thus, the energy density flow is calculated as

$\mathbf{S}=-\partial_{t}\psi\left(\nabla_{\perp}\psi-\frac{\partial_{t}\psi}{c}\mathbf{e}_{z}\right),$
(17)

Upon Noether’s theorem

$\nabla\cdot\mathbf{S}+\partial_{t}U=0.$ (18)

The rest of the stress-energy tensor allows us to define the conservation of
momentum. The momentum density is

$\mathbf{p}=\left[\frac{2}{c^{2}}(\partial_{t}\psi)+\frac{1}{c}(\partial_{z}\psi)\right]\nabla\psi.$
(19)

The momentum density flow is a three tensor

$\bm{M}=[\mathbf{M}_{x},\mathbf{M}_{y},\mathbf{M}_{z}]$ (20)

with columns

$\displaystyle\mathbf{M}_{x}$ $\displaystyle=$
$\displaystyle-(\partial_{x}\psi)\left[\nabla_{\perp}\psi-\frac{1}{c}(\partial_{t}\psi)\mathbf{e}_{z}\right]-\mathscr{L}\mathbf{e}_{x},$
(21) $\displaystyle\mathbf{M}_{y}$ $\displaystyle=$
$\displaystyle-(\partial_{y}\psi)\left[\nabla_{\perp}\psi-\frac{1}{c}(\partial_{t}\psi)\mathbf{e}_{z}\right]-\mathscr{L}\mathbf{e}_{y},$
(22) $\displaystyle\mathbf{M}_{z}$ $\displaystyle=$
$\displaystyle-(\partial_{z}\psi)\nabla_{\perp}\psi-\frac{1}{c^{2}}(\partial_{t}\psi)\mathbf{e}_{z}-\mathscr{L}\mathbf{e}_{z}.$
(23)

Noether’s theorem states that

$\partial_{t}\mathbf{p}=-\nabla\bm{M},$ (24)

where $\nabla$ is a row vector in this last equation.

Momentum analysis can be taken further than the coverage given in this
document. Other relevant analysis, such as angular momentum, derive directly
from ours since angular momentum density is derived from the momentum density
by its definition.

$\mathbf{l}=\mathbf{r}\times\mathbf{p}$ (25)

One of the most useful calculation could be the $z$ component of the angular
momentum Allen  _et al._ (2016). In this case, the angular momentum density
with regard to the $z$ axis in cylindrical coordinates is,

$l_{z}=(\partial_{\phi}\psi)\left[\frac{2}{c^{2}}(\partial_{t}\psi)+\frac{1}{c}(\partial_{z}\psi)\right].$
(26)

This quantity is conserved by means of the angular momentum density flow
$\mathbf{M}_{\phi}$, whose expression is

$\mathbf{M}_{\phi}=-(\partial_{\rho}\psi)(\partial_{\phi}\psi)\mathbf{e}_{\rho}+\left[\frac{1}{2}\left(\partial_{\rho}\psi-\frac{\partial_{\phi}\psi}{\rho^{2}}\right)-\frac{\partial_{t}\psi}{c}\left(\frac{\partial_{t}\psi}{c}+\partial_{z}\psi\right)\right]\mathbf{e}_{\phi}+\frac{\partial_{t}\psi}{c}\partial_{\phi}\psi\mathbf{e}_{z}.$
(27)

The conservation equation of the angular momentum around the $z$ axis is

$\frac{1}{\rho}\partial_{\rho}(\rho\mathbf{M}_{\phi}\cdot\mathbf{e}_{\rho})+\partial_{\phi}(\mathbf{M}_{\phi}\cdot\mathbf{e}_{\phi})+\partial_{z}(\mathbf{M}_{\phi}\cdot\mathbf{e}_{z})+\partial_{t}l_{z}=0$
(28)

We give the previous equation to help further studies with regard to orbital
angular momentum of light beams.

By comparing Eqs. (19) and (17) it can be seen that $\mathbf{S}$ and
$\mathbf{p}$ are not antiparallel, in general, within the paraxial framework.
Note that this is a key difference with regard to d’Alembert equation, as
shown in Eqs. (89) and (94) where energy density flow and momentum have a
simple relation of proportionality between them. A solution of the paraxial
wave equation can be considered as a tool for describing propagation if it
behaves approximately as a solution of the d’Alembert wave equation.
Therefore, any parameter related to $\mathbf{S}$ and/or $\mathbf{p}$ can be
seen as a measure of the quality of a paraxial solution. For instance,

$\mathbf{D}=c^{2}\mathbf{p}+\mathbf{S}$ (29)

will measure where a solution of the paraxial wave equation can be a good
approximation to wave propagation. Note that the last quantity is identically
zero for any solution of the d’Alembert equation. Thus, the paraxial quality
of a solution can be obtained by comparing $\mathbf{D}$ with $\mathbf{S}$; for
instance

$Q=\frac{\|\mathbf{D}\|}{\|\mathbf{S}\|}.$ (30)

##  IV Time harmonic solutions

We let $\psi$ have a time harmonic dependence. Hence,

$\psi(\mathbf{r},t)=\tilde{\psi}(\mathbf{r})e^{-i\omega t},$ (31)

where $\tilde{\psi}(\mathbf{r})$ is a solution of (3). To recover the
practical solution of (4) one has to extract, for example, the real part of
(31),

$\psi=\Re\left\\{\tilde{\psi}e^{-i\omega
t}\right\\}=\frac{1}{2}\left(\tilde{\psi}\,e^{-i\omega
t}+\tilde{\psi}^{*}e^{i\omega t}\right),$ (32)

where the asterisk superscript denotes the complex conjugate. This form of
expressing the paraxial solution leads to a straightforward calculation of
time average energy and flow densities. As usual, time average value of a
given time harmonic function $A(t)$ is obtained by doing the integral

$\langle A\rangle=\frac{1}{T}\int_{<T>}\mathrm{d}t\,A(t),$ (33)

where $T=\frac{2\pi}{\omega}$. Let us start from energy density flow
calculation. First we find that

$\partial_{t}\psi=\frac{\omega}{2i}\left(\tilde{\psi}\,e^{-i\omega
t}-\tilde{\psi}^{*}e^{i\omega t}\right).$ (34)

The expressions of the partial derivatives with regard to the spatial
coordinates are also needed:

$\partial_{\mu}\psi=\frac{1}{2}\left(\partial_{\mu}\tilde{\psi}\,e^{-i\omega
t}+\partial_{\mu}\tilde{\psi}^{*}e^{i\omega t}\right),\quad\mu=x,y,z.$ (35)

Thus, energy density flow for time harmonic solutions is

$\mathbf{S}=\frac{\omega}{2}\left[\Im\left\\{\tilde{\psi}^{*}\left(\nabla_{\perp}\tilde{\psi}+ik\tilde{\psi}\mathbf{e}_{z}\right)\right\\}-\Im\left\\{\tilde{\psi}\left(\nabla_{\perp}\tilde{\psi}+ik\tilde{\psi}\mathbf{e}_{z}\right)e^{-i2\omega
t}\right\\}\right]$ (36)

The last expression is adequate for calculating time average density flow

$\langle\mathbf{S}\rangle=\frac{\omega}{2}\Im\left\\{\tilde{\psi}^{*}\left(\nabla_{\perp}\tilde{\psi}+ik\tilde{\psi}\mathbf{e}_{z}\right)\right\\}.$
(37)

The same procedure is performed with respect to the energy density $U$. The
quantities involved in the calculation have been already obtained, and
therefore

$U=\frac{k^{2}}{2}|\tilde{\psi}|^{2}+\frac{1}{4}\left\|\nabla_{\perp}\tilde{\psi}\right\|^{2}+\frac{1}{2}\Re\left\\{\left\\{\frac{1}{2}\left[(\partial_{x}\tilde{\psi})^{2}+(\partial_{y}\tilde{\psi})^{2}\right]-k^{2}\tilde{\psi}^{2}\right\\}e^{-i2\omega
t}\right\\}.$ (38)

Time average energy density for time harmonic solutions is obtained from this
last expression as,

$\langle
U\rangle=\frac{1}{2}\left(k^{2}|\tilde{\psi}|^{2}+\frac{1}{2}\left\|\nabla_{\perp}\tilde{\psi}\right\|^{2}\right).$
(39)

The formalism followed gives rise to the definition of a complex flow density
for $\tilde{\psi}$

$\mathbf{S}_{c}=-\frac{i\omega}{2}\tilde{\psi}^{*}(\nabla_{\perp}\tilde{\psi}+ik\tilde{\psi}\,\mathbf{e}_{z}.).$
(40)

It follows from its definition that

$\langle\mathbf{S}\rangle=\Re\\{\mathbf{S}_{c}\\}.$ (41)

The imaginary part of $\mathbf{S}_{c}$ will lead to the reactive energy
balance associated to the solution $\psi$.

Momentum density can also be calculated for time harmonic solutions
explicitly. The procedure is quite similar to the one carried out to obtain
$\mathbf{S}$ and yields

$\mathbf{p}=\frac{\omega}{c^{2}}\left[\Im\left\\{\tilde{\psi}\nabla\tilde{\psi}^{*}+\tilde{\psi}\nabla\tilde{\psi}e^{-i2\omega
t}\right\\}\right]+\frac{1}{2c}\left[\Re\left\\{\partial_{z}\tilde{\psi}\nabla\tilde{\psi}^{*}+\partial_{z}\tilde{\psi}\nabla\tilde{\psi}e^{-i2\omega
t}\right\\}\right].$ (42)

From this last equation, the average momentum density is

$\langle\mathbf{p}\rangle=\frac{\omega}{c^{2}}\Im\left\\{\tilde{\psi}\nabla\tilde{\psi}^{*}\right\\}+\frac{1}{2c}\Re\left\\{\partial_{z}\tilde{\psi}\nabla\tilde{\psi}^{*}\right\\}.$
(43)

Since energy density flow and momentum density do not follow a simple
proportionality relation in the paraxial approximation scheme, one issue about
using raytracing techniques may arise within this framework. In the full wave
propagation scheme, it can be shown how, for time harmonic regime, the phase
of a solution of the Helmholtz equation contains some information about the
vector behavior of average momentum density and energy density flow (see
Appendix A.1). From Eqs. (37) and (43), it can be seen how average energy flow
and momentum densities depend on the phase of $\tilde{\psi}$. The
decomposition of the phasor in terms of real amplitude $A(\mathbf{r})$ and
real phase $\phi(\mathbf{r})$ (the extraction of the constant $\phi_{0}$
allows us to define $A$ and $\phi$ as real functions),

$\tilde{\psi}=A(\mathbf{r})\exp\\{i[k\phi(\mathbf{r})+\phi_{0}]\\},$ (44)

in the aforementioned equations, gives the dependence of the average energy
density flow and momentum density in terms of amplitudes and phase of
$\tilde{\psi}$. For the average energy density flow, we obtain

$\langle\mathbf{S}\rangle=\frac{\omega^{2}}{2c}\left[A(\mathbf{r})\nabla_{\perp}\phi(\mathbf{r})+A^{2}(\mathbf{r})\mathbf{e}_{z}\right].$
(45)

The average momentum density in terms of amplitude and phase is

$\langle\mathbf{p}\rangle=\frac{1}{2c}\left\\{k^{2}[\partial_{z}\phi(\mathbf{r})-2A(\mathbf{r})]\nabla\phi(\mathbf{r})+\partial_{z}A(\mathbf{r})\nabla
A(\mathbf{r})\right\\}$ (46)

These last equations show that neither $\langle\mathbf{S}\rangle$ nor
$\langle\mathbf{p}\rangle$ are proportional, in general, to $\nabla\phi$.
Since we usually have the solution $\tilde{\psi}$ in this framework,
characteristics of $\langle\mathbf{S}\rangle$ and $\langle\mathbf{p}\rangle$
cannot be inferred by studying $\nabla\phi$. The energy density flow depends
on $\nabla_{\perp}\phi$ instead of $\nabla\phi$. The momentum density does
depend on $\nabla\phi$ but there is another term that contributes with $\nabla
A$. Thus, if the following condition is met,

$k^{2}[\partial_{z}\phi(\mathbf{r})-2A(\mathbf{r})]\gg\partial_{z}A(\mathbf{r}),$
(47)

the study of $\nabla\phi$ will help us to see features of the momentum density
approximately. The condition expressed in this last equation requires the
amplitude $A(\mathbf{r})$ to be a weak function of the coordinates and the
high frequency condition $k\gg 1$.

##  V Examples

###  V.1 Generic traveling wave towards $+z$

Let us compute energy density flow with Eq. (17) of the solution in Eq. (5);

$\mathbf{S}=c\,F(\mathbf{r}_{\perp})g^{\prime}(z-ct)\left[g(z-ct)\nabla_{\perp}F(\mathbf{r}_{\perp})+F(\mathbf{r}_{\perp})g^{\prime}(z-ct)\mathbf{e}_{z}\right].$
(48)

Also, the momentum density can be obtained

$\mathbf{p}=\frac{-1}{c}F(\mathbf{r}_{\perp})g^{\prime}(z-ct)\left[g(z-ct)\nabla_{\perp}F(\mathbf{r}_{\perp})+F(\mathbf{r}_{\perp})g^{\prime}(z-ct)\mathbf{e}_{z}\right].$
(49)

These expressions are also obtained if Eqs. (89) and (94) are used instead.
Note that, since the solution in (5) is also a solution to Eq. (87) the
procedure explained can be done as well. It follows that the energy and
momentum balance obtained in both cases is exactly the same, as it should be.

###  V.2 Paraxial plane waves

Under the time dependence $e^{-i\omega t}$, the phasor associated to a plane
wave is

$\tilde{\psi}=\exp(i\mathbf{k}\cdot\mathbf{r}),$ (50)

where

$\mathbf{k}=k\left[q\left(\cos\beta\,\mathbf{e}_{x}+\sin\beta\,\mathbf{e}_{y}\right)+\left(1-\frac{q^{2}}{2}\right)\mathbf{e}_{z}\right].$
(51)

Phase fronts are planes normal to the direction of $\mathbf{k}$. The zenith
angle of $\mathbf{k}$ (angle with the $+z$ axis) is

$\alpha=\arccos\left(\frac{1-\frac{q^{2}}{2}}{\sqrt{1+\frac{q^{4}}{4}}}\right),$
(52)

and the azimuthal angle is $\beta$. The complex energy density flow is a
constant vector for such solution

$\mathbf{S}_{c}=\frac{k^{2}c}{2}\left[\mathbf{n}+\frac{q^{2}}{2}\mathbf{e}_{z}\right],$
(53)

where

$\mathbf{n}=\frac{\mathbf{k}}{\|\mathbf{k}\|},$ (54)

as happens with propagative usual plane waves. For real parameters $q$ and
$\beta$, the vector is a real one and there is not reactive energy density.
The average momentum density can also be computed:

$\langle\mathbf{p}\rangle=\frac{-k^{2}}{2c}\left(1+\frac{q^{2}}{2}\right)\mathbf{n}.$
(55)

For this kind of solutions with constant amplitude, the average momentum
density is parallel to $\mathbf{k}$ as predicted by Eq. (46). Conversely, the
zenith angle of $\mathbf{S}_{c}$ is

$\alpha_{\mathbf{S}}=\arccos\left(\frac{1}{\sqrt{1+q^{2}}}\right)$ (56)

which, unlike usual homogeneous Helmholtz plane waves, is different from
$\alpha$. The angle $\alpha_{\mathbf{S}}$ is in $[0,\frac{\pi}{2})$ as $q$ is
in $[0,+\infty)$ (this set of $q$ leads to a basis in terms of the Fourier
transform for the paraxial equation). In other words, $\mathbf{S}_{c}$ cannot
have components towards $-z$. The norms of both $\langle\mathbf{p}\rangle$ and
$\mathbf{S}_{c}$ are functions of the direction which is also a difference
with respect to Helmholtz plane waves. This fact can be explained because the
medium gets denser as the zenith angle of $\mathbf{k}$ increases,

$\|\mathbf{k}\|=k\sqrt{1+\frac{q^{4}}{4}}.$ (57)

As a consequence, the phase velocity

$v_{p}\leq c,$ (58)

for paraxial plane waves. Besides, except for $q=0$ constant phase curves are
not orthogonal to energy density flow as was previously predicted. Momentum
density can be written in terms of energy density flow as,

$\langle\mathbf{p}\rangle=\frac{1}{c^{2}}\left(\mathbf{D}-\mathbf{S}_{c}\right)$
(59)

where $\mathbf{D}$ is obtained as

$\mathbf{D}=\frac{k^{2}c}{2}\frac{q^{2}}{2}\left(\mathbf{n}-\mathbf{e}_{z}\right).$
(60)

Note that as $q$ increases, the vector $\mathbf{D}$ gets larger, showing that
the idea of measuring the quality of the paraxial solution based on this
vector is a fair procedure. The case $q=0$ is a special case of the example in
Subsection V.1.

###  V.3 Green’s function of the time harmonic paraxial wave equation

Green’s function associated to the operator in Eq. (3) is

$\tilde{\psi}=\frac{1}{4\pi
z}\exp\left[ik\left(z+\frac{\rho^{2}}{2z}\right)\right].$ (61)

The complex energy density flow can be obtained by using Eq. (40),

$\displaystyle\mathbf{S}_{c}$ $\displaystyle=$
$\displaystyle\frac{k^{2}c}{32\pi^{2}}\frac{1}{z^{2}}\left(\frac{\rho}{z}\mathbf{e}_{\rho}+\mathbf{e}_{z}\right),\quad(\textrm{cyl.});$
(62a) $\displaystyle\mathbf{S}_{c}$ $\displaystyle=$
$\displaystyle\frac{k^{2}c}{32\pi^{2}}\frac{1}{r^{2}\cos^{3}\theta}\mathbf{e}_{r},\quad(\textrm{sph.}).$
(62b)

There is not reactive energy in this solution since the complex energy density
flow is a real quantity. As shown in Eq. (62b), field lines of
$\mathbf{S}_{c}$ are radial. In the semispaces of $z<0$, the energy goes
towards the origin, where it concentrates and then, the energy is radiated
towards the $z>0$ semispace. Formally, this statement can be shown by means of
the divergence. It is easy to see that, in general,

$\nabla\cdot\mathbf{S}_{c}=0.$ (63)

From the viewpoint of the energy, this last equation shows that there is not
energy storage, only energy flow. The origin in such solutions must be treated
by means of the definition of the divergence. Let $B_{R}(0)$ be a ball with
center at the origin and radius $R$. Then

$\left.\nabla\cdot\mathbf{S}_{c}\right|_{\mathrm{origin}}=\lim_{R\to
0}\frac{3}{4\pi}\frac{1}{R^{3}}\iint_{\partial
B_{R}(0)}\mathrm{d}a\,\mathbf{S}_{c}\cdot\mathbf{e}_{r}.$ (64)

Using spherical coordinates,

$\left.\nabla\cdot\mathbf{S}_{c}\right|_{\mathrm{origin}}=\frac{3k^{2}c}{128\pi^{3}}\lim_{R\to
0}\frac{1}{R^{3}}\int_{0}^{\pi}R\mathrm{d}\theta\int_{0}^{2\pi}R\sin\theta\mathrm{d}\varphi\,\frac{1}{R^{2}\cos^{3}\theta}.$
(65)

By operating in this last equation, it is found that,

$\left.\nabla\cdot\mathbf{S}_{c}\right|_{\mathrm{origin}}=\frac{3k^{2}c}{64\pi^{2}}\lim_{R\to
0}\frac{1}{R^{3}}\int_{0}^{\pi}\mathrm{d}\theta\,\frac{\sin\theta}{\cos^{3}\theta}.$
(66)

The integral in this last equation is equal to zero and thus,

$\left.\nabla\cdot\mathbf{S}_{c}\right|_{\mathrm{origin}}=0,$ (67)

showing that Eq. (63) holds for every point in the space. This solution is a
very special one because, from the viewpoint of the energy, it lacks any
source even though it was generated by an impulse function. Therefore, by
means of the convolution integral, any radiation problem in time harmonic
regime is automatically transformed into a propagation problem without
sources. It follows, that this feature is also naturally inherited by
solutions of general time dependence.

The average momentum density of Green’s function can be obtained as

$\langle\mathbf{p}\rangle=\frac{-k^{2}}{32c\pi^{2}}\frac{1}{z^{2}}\left(1+\frac{\rho^{2}}{2z^{2}}\right)\left[\frac{\rho}{z}\mathbf{e}_{\rho}+\left(1-\frac{\rho^{2}}{2z^{2}}-\frac{k^{-2}}{z^{2}+\frac{\rho^{2}}{2}}\right)\mathbf{e}_{z}\right].$
(68)

Note that, it is easy to write $\langle\mathbf{p}\rangle$ in terms of
$\mathbf{S}_{c}$ plus another term that accounts for the difference between
them,

$\langle\mathbf{p}\rangle=-\frac{1}{c^{2}}\mathbf{S}_{c}+\frac{k^{2}}{32c\pi^{2}}\frac{1}{z^{2}}\left[\frac{-\rho^{3}}{2z^{3}}\mathbf{e}_{\rho}+\left(\frac{\rho^{4}}{4z^{4}}+\frac{k^{-2}}{z^{2}+\frac{\rho^{2}}{2}}+\frac{k^{-2}\rho^{2}}{2z^{4}+z^{2}\rho^{2}}\right)\mathbf{e}_{z}\right].$
(69)

Therefore, the proposed vectorial quantity defined in Eq. (29) is

$\mathbf{D}=\frac{k^{2}c}{32\pi^{2}}\frac{1}{z^{2}}\left[\frac{-\rho^{3}}{2z^{3}}\mathbf{e}_{\rho}+\left(\frac{\rho^{4}}{4z^{4}}+\frac{k^{-2}}{z^{2}+\frac{\rho^{2}}{2}}+\frac{k^{-2}\rho^{2}}{2z^{4}+z^{2}\rho^{2}}\right)\mathbf{e}_{z}\right].$
(70)

This last equation shows that momentum density and energy density flow are not
antiparallel in general. In high frequency conditions ($k\gg 1$), $\mathbf{D}$
can be neglected with regard to $\mathbf{S}_{c}$ in the regions where
$|z|\gg\rho$, that is, points near the axial direction. This is a well known
result since the paraxial Green’s function can be obtained by the appropriate
approximations of the distance in the Helmholtz Green’s function. Again, this
result shows that the comparison between momentum density and energy density
flow gives a fair procedure to obtain a measure of the quality of a paraxial
solution. Last equations are rather cumbersome to establish further analytical
differences between momentum density and energy density flow. Nonetheless,
since we know that energy density flow of Green’s function is radial, these
differences can be seen graphically.

Figure 1: Momentum density of paraxial Green’s function: Field lines and
intensity. The level of gray indicates the norm of the momentum density at
that point, being black the minimum values (at $z=\pm k^{-1}$ in the z axis)
and white the maximum ones (concentrated in the $z=0$ plane). The intensity is
in logarithmic scale.

In Fig. 1 a sketch of momentum density of paraxial Green’s function is shown.
All the field lines start at the point $(\rho=0,z=-k^{-1})$ and end at
$(\rho=0,z=k^{-1})$ and cross the $z=0$ plane in its normal direction. At this
plane, momentum density and energy density flow are ill-defined just like the
phasor. No other points, except those fulfilling $|z|\gg\rho$ and $|z|>k^{-1}$
can be seen as a region where the paraxial Green’s function is a fair tool for
representing wave propagation. Note that the origin is excluded from this
region since near this region $\mathbf{S}_{c}$ and $\langle\mathbf{p}\rangle$
are parallel and not antiparallel.

To finish this example, let us study the field of phase velocity of the
paraxial Green’s function. The objective is to see if the phase of the
solution has relevant information. From the previous example, we found that
phase velocities below the propagation speed in the medium were natural in
this scheme of propagation. The phase velocity is a field that can be defined
as

$\mathbf{v}_{p}=\omega\frac{\nabla\angle\\{\tilde{\psi}\\}}{\left\|\nabla\angle\\{\tilde{\psi}\\}\right\|^{2}}.$
(71)

In this case,

$\angle\\{\tilde{\psi}\\}=k\left(z+\frac{\rho^{2}}{2z}\right);$ (72)

and therefore, we obtain

$\mathbf{v}_{p}=c\left(1+\frac{\rho^{4}}{4z^{4}}\right)^{-1}\left[\frac{\rho}{z}\mathbf{e}_{\rho}+\left(1-\frac{\rho^{2}}{2z^{2}}\right)\mathbf{e}_{z}\right].$
(73)

By comparing this last expression to Eq. (68), it can be seen that under high
frequency conditions $\mathbf{v}_{p}$ is antiparallel to
$\langle\mathbf{p}\rangle$ as was predicted by Eq. (47). This feature can also
be seen in Fig.2.

Figure 2: Phase velocity field of the paraxial Green’s function. Line fields
are ray trajectories and the level of gray indicates the norm of the phase
velocity, being black null phase velocity and white wave propagation speed in
the medium. The general form of phase velocity confirms that, apparently there
is not a source at the origin. The wave travels from the $z=0^{-}$ plane where
phase velocity vanishes, concentrates in the origin as if there were a source-
sink combination and, then travels to the $z=0^{+}$ plane

The comparison of the field lines in Figs. 1 and 2 not only shows that, under
high frequency condition, momentum density and phase velocity are
antiparallel, but also where weak dependence of the amplitude function
$A(\mathbf{r})$ requirement holds in this example. Note that phase velocity in
the paraxial Green’s function also fulfills

$\|\mathbf{v}_{p}\|\leq c,$ (74)

as happened with paraxial plane waves.

###  V.4 Fundamental Gaussian beam

The phasor associated to the fundamental paraxial Gaussian beam can be written
in terms of the Rayleigh distance $b$ as follows

$\tilde{\psi}=\frac{1}{z-ib}\exp\left[ik\left(z+\frac{1}{2}\frac{\rho^{2}}{z-ib}\right)\right].$
(75)

Thus, the complex energy density flow is found by applying Eq. (40) as

$\mathbf{S}_{c}=\frac{\omega}{2}k|\tilde{\psi}|^{2}\left(\frac{x}{z-ib}\mathbf{e}_{x}+\frac{y}{z-ib}\mathbf{e}_{y}+\mathbf{e}_{z}\right).$
(76)

This solution, unlike the previous examples, has reactive energy density flow,

$\mathbf{Q}=\Im\\{\mathbf{S}_{c}\\}=\frac{k^{2}c}{2}|\tilde{\psi}|^{2}\frac{b\rho}{z^{2}+b^{2}}\mathbf{e}_{\rho}.$
(77)

From Eq. (76), the average energy density flow is

$\langle\mathbf{S}\rangle=\Re\\{\mathbf{S}_{c}\\}=\frac{k^{2}c}{2}|\tilde{\psi}|^{2}\left(\frac{\rho
z}{z^{2}+b^{2}}\mathbf{e}_{\rho}+\mathbf{e}_{z}\right).$ (78)

A sketch of the average energy density flow can be seen in Fig. 3

Figure 3: Example of energy density flow of a Gaussian beam. The level of gray
indicates the value of $\langle\mathbf{S}\rangle$: the lighter gray, the
higher values. It is straightforward to see from Eq. (78) that field lines in
the $\rho-z$ plane are one branch hyperbolas.

The average momentum density can be calculated as well, yielding

$\displaystyle\langle\mathbf{p}\rangle=-\frac{k^{2}}{2c}$
$\displaystyle|\tilde{\psi}|^{2}$ $\displaystyle\left\\{\frac{\rho
z}{z^{2}+b^{2}}\left(1+\frac{1}{2}\frac{\rho^{2}}{z^{2}+b^{2}}\right)\mathbf{e}_{\rho}+\ldots\phantom{\left(\frac{1}{1}\right)^{2}}\right.$
$\displaystyle\ldots$ $\displaystyle+$
$\displaystyle\left.\left[1+\frac{2kb-1}{k^{2}}\frac{1}{z^{2}+b^{2}}-\left(\frac{1}{2}\frac{\rho^{2}}{z^{2}+b^{2}}\right)^{2}\right]\mathbf{e}_{z}\right\\},$
(79)

which is a rather cumbersome expression. An example of momentum density of
this solution can be seen in Fig. 4. The resemblance with energy density flow
is rather high in contrast to what happened in the previous example.

Figure 4: Example of momentum density of a Gaussian beam. The level of gray
indicates the value of $\langle\mathbf{p}\rangle$: the lighter gray, the
higher values.

In any case, as in the previous example, we can extract the average energy
density flow dependence in the momentum density, yielding the vector
$\mathbf{D}$ that measures the quality of the solution at each point,

$\mathbf{D}=\frac{k^{2}c}{2}|\tilde{\psi}|^{2}\left\\{\frac{-\rho^{3}z}{2(z^{2}+b^{2})^{2}}\mathbf{e}_{\rho}+\left[\frac{1-2kb}{k^{2}}\frac{1}{z^{2}+b^{2}}+\left(\frac{1}{2}\frac{\rho^{2}}{z^{2}+b^{2}}\right)^{2}\right]\mathbf{e}_{z}\right\\}.$
(80)

Inasmuch as $\mathbf{D}$ can be neglected with respect to
$\langle\mathbf{S}\rangle$, the Gaussian beam can be considered, in principle,
as a fair approximation to the full wave propagation model. At first glance,
$\mathbf{D}$ vanishes if $z\gg\rho$ and $k\gg 1$, which are the classical
“near the axis” and “high frequency” approximations, respectively.
Nonetheless, this solution is required to represent a beam at its waist: in
the $z=0$ plane, under high frequency conditions, $\mathbf{D}$ can be
neglected with regard to $\langle\mathbf{S}\rangle$ for some “low” values of
$\rho$ but not the whole plane. All these issues are shown in Fig. 5 where the
parameter in Eq. (30) is taken as a quality measure for each point.

Figure 5: Quality in an example of Gaussian beam. The level of gray indicates
the quality being the lighter gray, the better quality.

It is noteworthy that points near the origin in Fig. 5 are not within the
lighter ones in this example where the high frequency condition is not
fulfilled. If the high frequency condition is met, all the points “near”
$\rho=0$ in the $z=0$ plane will raise their quality. Clearly, only a centered
disk at the origin in the $z=0$ plane can be considered as a good
approximation for propagation purposes but not the entire plane.

The study of the phase velocity is worth it in this example also. By taking
the definition in Eq. (71), the gradient of the phase of the Gaussian beam in
Eq. (75) is needed to obtain the field of phase velocity. Thus, in this
example,

$\nabla\angle\\{\tilde{\psi}\\}=k\left\\{\frac{z\rho}{z^{2}+b^{2}}\mathbf{e}_{\rho}+\left[1+\left(\frac{\rho^{2}}{2}-\frac{b}{k}\right)\frac{1}{z^{2}+b^{2}}-\frac{z^{2}\rho^{2}}{(z^{2}+b^{2})^{2}}\right]\mathbf{e}_{z}\right\\},$
(81)

which yields

$\mathbf{v}_{p}=c\frac{\frac{z\rho}{z^{2}+b^{2}}\mathbf{e}_{\rho}+\left[1+\left(\frac{\rho^{2}}{2}-\frac{b}{k}\right)\frac{1}{z^{2}+b^{2}}-\frac{z^{2}\rho^{2}}{(z^{2}+b^{2})^{2}}\right]\mathbf{e}_{z}}{\frac{z^{2}\rho^{2}}{(z^{2}+\rho^{2})^{2}}+\left[1+\left(\frac{\rho^{2}}{2}-\frac{b}{k}\right)\frac{1}{z^{2}+b^{2}}-\frac{z^{2}\rho^{2}}{(z^{2}+b^{2})^{2}}\right]^{2}}.$
(82)

The evaluation of the phase velocity at the axis is more straightforward,

$\left.\mathbf{v}_{p}\right|_{\rho=0}=c\frac{z^{2}+b^{2}}{z^{2}+b^{2}-\frac{b}{k}}\mathbf{e}_{z}.$
(83)

From this last equation, it is easy to see that

$\left.\|\mathbf{v}_{p}\|\right|_{\rho=0}>c.$ (84)

In other words, near the axis, the phase velocity only tends to the
propagation speed in the medium asymptotically. Let us evaluate phase velocity
at the plane of the beam waist ($z=0$), obtaining

$\left.\mathbf{v}_{p}\right|_{z=0}=\frac{c}{1+\frac{\rho^{2}}{2b^{2}}-\frac{1}{kb}}\mathbf{e}_{z}.$
(85)

From this last equation, one may see that

$\left.\|\mathbf{v}_{p}\|\right|_{\rho=w_{0},z=0}=c,$ (86)

where $w_{0}=\sqrt{\frac{2b}{k}}$ is the parameter that defines the beam width
at its waist. Phase velocity at points of the plane out of this circumference
is lower than $c$ whilst, at points contained by the circumference, the phase
velocity will be higher than the propagation speed in the medium. An example
of the field of phase velocity is shown in Fig. 6.

Figure 6: (Color) Example of phase velocity of a Gaussian beam. Phase velocity
below $c$ is shown in gray as in the previous example, while values over $c$
are drawn in red. The locus where $\|\mathbf{v}_{p}\|=c$ is shown in a dashed
line.

Note that the locus $\|\mathbf{v}_{p}\|=c$ does not collapse into the $z$
axis. Phase velocity above the propagation speed in the medium is a feature
that did not arise in the previous examples. Nonetheless, as it happens in the
full wave propagation framework, the presence of reactive energy is associated
to phase velocity greater than the propagation speed in the medium. Note that,
on account of the properties of the paraxial Green’s function with respect to
the sources, unlike Helmholtz waves, reactive energy is not associated to a
given source distribution, either real or virtual. It is possible to reproduce
the behavior of a Gaussian beam around its axis in a solution of the Helmholtz
equation, giving rise to a complex source point solution Couture and Bélanger
(1981); González-Morales  _et al._ (2010). This reconstruction method is a
valid one with the criteria exposed in Mahillo-Isla and González-Morales
(2016). But, the complex source point solution has to be regularized near the
ring determined by the branch point of the complex distance function Kaiser
(2003); Tagirdzhanov and Kiselev (2015). This fact shows that the reactive
power in the paraxial framework does not present a problem in terms of
regularity of the solutions but its counterpart in the Helmholtz scheme
presents some difficulties.

##  VI Conclusions

The application of the _slowly varying envelope approximation_ to the
d’Alembert equation results in another hyperbolic problem that still holds
invariant quantities in terms of energy and momentum. It has been shown that
these quantities do not follow the expressions of their counterparts in the
full wave propagation scheme. The comparison of these invariants with the
expected behavior in terms of the full propagation framework allows us to
obtain quality parameters of paraxial solutions. We have built a vector one
that takes into account the momentum density and energy density flow. Although
the quality is measured locally, a global parameter could be obtained by
integration. Nonetheless, our local measure of quality shows that the high-
frequency condition turns out to be relevant. We have gained new insights with
regard to the behavior of paraxial solutions: rays (defined by phase fronts)
are parallel neither to energy flow nor to momentum. Phase velocity can be
either higher or lower than the propagation speed in the medium. Phase
velocity above the propagation speed in the medium is associated with the
presence of reactive energy, as happens in the d’Alembert equation, while
phase velocity below the propagation speed seems to be “natural” in paraxial
propagation. Therefore, general phase matching when reconstructing a solution
form a paraxial one is a delicate procedure. The study of paraxial Green’s
function in time harmonic regime showed that any problem with sources turns
into a simple propagation problem. Altogether, the paraxial propagation
framework proposes a sourceless framework where the medium seems to be more
dispersive than in the full propagation scheme. The study of the Gaussian beam
shows that, as a solution it may have a remarkable high quality region around
its axis. The analysis of the phase velocity shows that it exceeds the
propagation speed in the medium in the high quality region. Although the
particular feature with regard to the sources in paraxial solutions does not
present a problem, it follows that a solution of the d’Alembert equation may
need a non-physical distribution of sources in order to maintain that phase
speed far from them.

###### Acknowledgements.

This work is supported by the Ministerio de Econom a y Competitividad (MINECO)
(MINECO/FEDER, EU) (TEC2015-69665-R). The authors want to thank Prof. Carlos
Dehesa-Martínez not only for his help in this article but also for his
knowledge, friendship and patience during the last years. Without his
guidance, the good part of this paper would have never been possible.

##  Appendix A Summary of results for the d’Alembert equation

A Lagrangian density associated to the d’Alembert equation

$\nabla^{2}\psi-\frac{1}{c^{2}}\partial_{tt}\psi=0$ (87)

can be found,

$\mathscr{L}=\frac{1}{2}\frac{(\partial_{t}\psi)^{2}}{c^{2}}-\frac{1}{2}\left[(\partial_{x}\psi)^{2}+(\partial_{y}\psi)^{2}+(\partial_{z}\psi)^{2}\right].$
(88)

Thus, the energy density flow is

$\mathbf{S}=-\partial_{t}\psi\,\nabla\psi.$ (89)

The phasor associated to a time harmonic solution of the d’Alembert equation
must fulfill the well-known Helmholtz equation

$\nabla^{2}\tilde{\psi}+k^{2}\tilde{\psi}=0.$ (90)

In this case,the energy density flow depends only on the phasors as

$\mathbf{S}=\frac{\omega}{2}\left(\Im\left\\{\tilde{\psi}^{*}\nabla\tilde{\psi}\right\\}-\Im\left\\{\tilde{\psi}\nabla\tilde{\psi}\,e^{-i2\omega
t}\right\\}\right).$ (91)

The average energy density flow is

$\langle\mathbf{S}\rangle=\frac{\omega}{2}\Im\left\\{\tilde{\psi}^{*}\nabla\tilde{\psi}\right\\}.$
(92)

This leads to the definition of the complex energy density flow as

$\mathbf{S}_{c}=-i\frac{\omega}{2}\tilde{\psi}^{*}\nabla\tilde{\psi},$ (93)

It follows from its construction that the real part of $\mathbf{S}_{c}$ is the
average energy density flow. The imaginary part gives the reactive energy
density associated to the solution whose phasor is $\tilde{\psi}$. The
momentum density associated to a solution of the d’Alembert equation is
obtained as follows from the Noether’s theorem

$\mathbf{p}=\frac{1}{c^{2}}\partial_{t}\psi\,\nabla\psi,$ (94)

which is almost equal to the energy density flow in terms of the dependence
with $\psi$, see Eq. (89). From Eqs. (89) and (94) it is clear that
$\mathbf{S}$ and $\mathbf{p}$ are antiparallel vectors. Momentum can be also
calculated explicitly for time harmonic solutions. The procedure is quite
similar to $\mathbf{S}$ and produces

$\mathbf{p}=\frac{1}{c^{2}}\frac{\omega}{2}\left(\Im\left\\{\tilde{\psi}\nabla\tilde{\psi}^{*}\right\\}+\Im\left\\{\tilde{\psi}\nabla\psi\,e^{-i2\omega
t}\right\\}\right).$ (95)

From this last equation, the average momentum is

$\langle\mathbf{p}\rangle=\frac{1}{c^{2}}\frac{\omega}{2}\Im\left\\{\tilde{\psi}\nabla\tilde{\psi}^{*}\right\\}.$
(96)

###  A.1 The connection of energy density flow and momentum with raytracing
techniques

For time harmonic waves, the study of phasors gives an idea of momentum and
energy flow. A phasor has the general form of

$\tilde{\psi}(\mathbf{r})=A(\mathbf{r})\exp[i\Phi(\mathbf{r})],$ (97)

where $A(\mathbf{r})$ and $\Phi(\mathbf{r})$ are the amplitude and phase
functions respectively. By substituting this last expression into Eqs. (92)
and (96), it follows that

$\left.\begin{array}[]{c}\langle\mathbf{S}\rangle\\\
\langle\mathbf{p}\rangle\end{array}\right\\}\propto|A(\mathbf{r})|^{2}\,\nabla\Phi(\mathbf{r}).$
(98)

Thus,

$\left.\begin{array}[]{c}\langle\mathbf{S}\rangle\\\
\langle\mathbf{p}\rangle\end{array}\right\\}\parallel\nabla\Phi(\mathbf{r}).$
(99)

This last equation explains why the eikonal equation and raytracing techniques
are so important. Relevant information of the vector structure of the energy
flow and momentum is contained in the phase of the phasor without even
assuming a high frequency condition. If the application of these techniques
gives an accurate description of $\Phi(\mathbf{r})$, by applying conservation
procedures both momentum and energy flow can be obtained and, by extension,
the solution $\psi$.

## References

  * Siegman (1986) A. E. Siegman,  _Lasers_ , 1st ed. (University Science Books, California, 1986).
  * Grella (1982) R. Grella, J. Opt. 13, 367 (1982).
  * Born and Wolf (1980) M. Born and E. Wolf,  _Principles of Optics_ , 6th ed. (Pergamon Press Ltd., Oxford, 1980).
  * Ion (2005) J. C. Ion,  _Laser Processing of Engineering Materials_ (Butterworth-Heinemann, Oxford, 2005).
  * Radziemski  _et al._ (1987) L. J. Radziemski, R. W. Solarz, and J. A. Paisner, eds.,  _Laser Spectroscopy and its Applications_ (Marcel Dekker Inc, New York, 1987).
  * Welch and van Gemert (2011) A. J. Welch and M. J. C. van Gemert, eds.,  _Optical-Thermal Response of Laser-Irradiated Tissue_ , 2nd ed. (Springer, Netherlands, 2011).
  * Mahillo-Isla and González-Morales (2015) R. Mahillo-Isla and M. J. González-Morales, J. Opt. Soc. Am. A 32, 1236 (2015).
  * Mahillo-Isla and González-Morales (2016) R. Mahillo-Isla and M. J. González-Morales, J. Opt. Soc. Am. A 33, 1735 (2016).
  * Mukunda  _et al._ (1985) N. Mukunda, R. Simon, and E. C. G. Sudarshan, J. Opt. Soc. Am. A 2, 416 (1985).
  * Simon  _et al._ (1986) R. Simon, E. C. G. Sudarshan, and N. Mukunda, J. Opt. Soc. Am. A 3, 536 (1986).
  * Vaveliuk  _et al._ (2007) P. Vaveliuk, B. Ruiz, and A. Lencina, Opt. Lett. 32, 927 (2007).
  * Allen  _et al._ (2016) L. Allen, S. M. Barnett, and M. J. Padgett,  _Optical Angular Momentum_ (CRC Press, Boca Raton, 2016).
  * Couture and Bélanger (1981) M. Couture and Pierre-A. Bélanger, Phys. Rev. A 24, 355 (1981).
  * González-Morales  _et al._ (2010) M. J. González-Morales, R. Mahillo-Isla, E. Gago-Ribas, and C. Dehesa-Martínez, J. of Electromagn. Waves and Appl. 24, 1103 (2010).
  * Kaiser (2003) G. Kaiser, J. Phys. A 36, R291 (2003).
  * Tagirdzhanov and Kiselev (2015) A. M. Tagirdzhanov and A. P. Kiselev, Opt. Spectrosc. 119, 257 (2015).


#  The Scattering of Dirac and Majorana Fermions in Spherically Symmetric
Gravitational Field and Torsion Field

Junhui Laia, Xun Xuea,b,111Corresponding<EMAIL_ADDRESS>
###### Abstract

The possibility to distinguish Dirac from Majorana fermion by gravitational
interaction and torsion coupling is discussed. The scattering of both Dirac
and Majorana particles are the same in the spherically symmetric gravitational
field up to the first order perturbation of the fermion-gravity interaction.
The scattering behaviors of Dirac and Majorana fermion are also the same by
axial vector torsion but different by vector torsion. We get the conclusion
that the scattering by vector torsion field can distinguish the Dirac from
Majorana neutrino.

aDepartment of Physics, East China Normal University, Shanghai 200241, China

bCenter for Theoretical Physics, Xinjiang University, Urumqi 830046, China

##  1 Introduction

Majorana Fermion is first proposed by Majorana in 1937 that the neutral spin
one half fermion can be described by the real valued wave function and
therefore the fermion is its own antiparticle[1]. The idea of Majorana nature
of neutrino is proposed by Furry in 1939[2], the neutrinoless double beta
decay, which only occurs for the Majorana type of neutrino particle. The
Standard Model fermions are all known to behave as Dirac particles at the
energy scale lower than the electroweak symmetry breaking scale except
neutrino, the nature of which is not known. The behavior of neutrino at low
energy scale depends on the fermion type of right-handed neutrino introduced
to explain neutrino oscillation. However the Majorana nature of neutrino
violates the lepton number conservation and even the $B-L$ conservation.
Though the neutrinoless double beta decay play an important role in solving
the absolute neutrino masses and the neutrino mass hierarchy problems, it has
not been discovered yet[3]. In principle, magnetic field could also
distinguish the Dirac or the Majorana type of neutrinos for their unique
magnetic moment structures, though the required fields are extremely strong.

The idea of distinguishing the types of neutrino by the their response to the
gravitational effect is proposed since 1990s[4, 5, 6, 7, 8]. The spin-gravity
effect is proposed to distinguish the fermion type of the neutrino by the
spin-flip matrix elements induced by the spin-gravity coupling in a Lense-
Thirring background which is different for Dirac type of neutrino from
Majorana one[6, 9]. However, the requirement of self charge conjugation of
spinor rather than the neutrino field can only apply to the quantum mechanical
Majorana type of wave packet rather than a Majorana particle in propagation,
which is a plane-wave expansion in quantum field theory approach just as
pointed out in reference[10]. The conclusion about the spin-gravity coupling
induced spin-flip matrix elements does not follow for the quantum field theory
approach and the difference between Majorana and Dirac particle is suppressed
by the factor $m/E$ and undetectable in the case[10].

As a matter of fact, a momentum eigenstate can not be charge conjugation
invariant at the same time in a quantum mechanical sense. There is the similar
problem in the wave function treatment of neutral scalar particle that there
is no momentum eigenfunction which is charge conjugation invariant. The proper
treatment framework is the second quantization scheme. The response of
neutrino to the gravity is very weak in whatever the direct gravitational
vertex coupling or the spin-gravity coupling sense so the effect can only be
significant in the astronomical scale propagation of the neutrino. The wave
packet approach with different momentum components is not applicable at such
large scale propagation. We tend to consider the scattering of asymptotic free
neutrino by a weak spherically symmetric gravitational field for both
torsionless and torsionness cases in this paper.

The connection is an independent degree of freedom of geometry in general. The
Platini formalism of general relativity constrains the connection of spacetime
to Levi-Civita one in the absence of spinor matter source. The Einstein-Cartan
spacetime is a natural extension of gravitation theory beyond general
relativity. The general relativity is usually considered as an effective
theory at low energy scale. The connection is constrained or condensed to
torsionless Levi-Civita one. At high energy, the connection may be an active
degree of freedom and deviates from Levi-Civita one by excitation[11].
Moreover, the non-trivial distribution of spin angular momentum tensor also
induces curved spacetime with torsion. Even in the classical gravitation
theory which is considered as the low energy effective theory of quantum
gravity, torsion can also play an important role. It is revealed that gravity
can be described equivalently by either the torsionless curved spacetime with
zero nonmetricity or the curvature-less spacetime with torsion and zero
nonmetricity or even both torsionless and curvature-less spacetime with
dynamical nonmetricity[12, 13]. In the teleparalelle gravity framework,
gravity is described by dynamical torsion of the flat spacetime with zero
nonmetricity[14].

Torsion can be decomposed into vector, axial vector and pure tensor parts
according to the irreducible representation of Lorentz group[14]. The pure
tensor part does not appear in the coupling of torsion with fermion and so
does the vector part in the minimal coupling scheme. However, the coupling
between torsion and fermion is an effective one at low energy scale in the
effective theory of high energy gravitational theory such as quantum gravity.
The non-minimal coupling between vector torsion and fermion is possible in
general[15]. Moreover, in the equivalent description of gravitation by torsion
rather than curvature, the vector torsion coupling with fermion does appear.
We take both vector torsion and axial vector torsion into consideration in
this paper.

The torsion free Levi-Civita connection in general relativity is the
consequence of both the absence of spinor source matter and the local Lorentz
invariance. Torsion naturally appears in a gravitation theory with Lorentz
violation even in the absence of spinor matter source[16, 17, 18, 19]. It is
proposed that the Lorentz violation of quantum gravity may be frozen into
large scale one by inflation and the gravitation theory governing the
evolution of the universe should take the large scale Lorentz violation into
account at the cosmic scale. The induced torsion can contribute an additional
energy-momentum distribution accompanying the cosmic media source and gives
rise to the accelerated expansion of the universe[20, 21, 22, 23, 24]. The
cosmic torsion can influence the propagation path of the neutrino traveling at
the cosmic distance. The possible different responce manners of Majorana and
Dirac particle to torsion may also distinguish the fermion type of neutrino.

##  2  Tetrad fields and spin connection

The tetrad fields ${e^{\mu}}_{a}(x),a=0,1,2,3$ are defined as

${e^{\mu}}_{a}{e^{\nu}}_{b}{g_{\mu\nu}}={\eta_{ab}}\;,$ (1)

where we use Greek and Latin letters to represent space-time and tangent space
index respectively. Then we can define Gamma matrices with space-time index

${\gamma^{\mu}}(x)={e^{\mu}}_{a}(x){\gamma^{a}}$ (2)

and spin connection

${\Gamma_{\mu
a}}^{b}={e_{\beta}}^{b}{\Gamma_{\mu\alpha}}^{\beta}{e^{\alpha}}_{a}+{e_{\alpha}}^{b}{\partial_{\mu}}{e^{\alpha}}_{a}={e_{\alpha}}^{b}{\nabla_{\mu}}{e^{\alpha}}_{a}\;.$
(3)

The connection can be decomposed into Levi-Civita one and contortion[25]

${\Gamma_{\mu\nu}}^{\lambda}={{\tilde{\Gamma}}_{\mu\nu}}^{\lambda}+{K_{\mu\nu}}^{\lambda}\;,$
(4)

where ${{\tilde{\Gamma}}_{\mu\nu}}^{\lambda}$ is the Christoffel symbol,
spacetime indexed Levi-Civata connection and ${K_{\mu\nu}}^{\lambda}$ is the
contortion tensor. It can be derived from the formula for the spin connection
(3) that

${\Gamma_{\mu ab}}={e_{\alpha
b}}{\nabla_{\mu}}{e^{\alpha}}_{a}={\nabla_{\mu}}\left({e_{\alpha
b}}{e^{\alpha}}_{a}\right)-{e^{\alpha}}_{a}{\nabla_{\mu}}{e_{\alpha
b}}=-{e^{\alpha}}_{a}{\nabla_{\mu}}{e_{\alpha b}}=-{\Gamma_{\mu ba}}\;,$ (5)

i.e. the last two indices of spin connection are anti-symmetric if metricity
condition ${\nabla_{\mu}}\left({e_{\alpha b}}{e^{\alpha}}_{a}\right)=0$ holds.
The metric preserving property of spin connection indicates that the spin
connection is Lorentzian, i.e. it takes value in the Lorentz Lie algebra. The
spin connection defines the Fock-Ivanenko covariant derivative

${D_{\mu}}={\partial_{\mu}}-\frac{i}{2}{\Gamma_{\mu ab}}{S^{ab}}\;,$ (6)

where the second part acts only on the tangent space indices and the Lorentz
group generators $S^{ab}$ act on field according to the specific Lorentz
representation that the field belongs. The Fock-Ivanenko covariant derivative
of a field transforming according to some specific representation of local
Lorentz group falls to the tensor product of vector and the representation of
the field in consideration so it is the covariant derivative of local Lorentz
gauge transformation in this sense. For example, the covariant derivative of
vector fields and spinor fields are

${D_{\mu}}{A^{d}}={\partial_{\mu}}{A^{d}}-\frac{i}{2}{\Gamma_{\mu
ab}}{({V^{ab}})_{c}}^{d}{A^{c}}={\partial_{\mu}}{A^{d}}+{\Gamma_{\mu
b}}^{d}{A^{b}}$ (7)

and

${D_{\mu}}\psi={\partial_{\mu}}\psi-\frac{i}{2}{\Gamma_{\mu ab}}{S^{ab}}\psi$
(8)

respectively, where

${({V^{ab}})_{c}}^{d}=i({\eta^{bd}}{\delta^{a}}_{c}-{\eta^{ad}}{\delta^{b}}_{c})$
(9)

and

${S^{ab}}=\frac{1}{2}{\sigma^{ab}}=\frac{i}{4}[{\gamma^{a}},{\gamma^{b}}]$
(10)

are the vector representation and spinor representation of generators
respectively. For scalar field $\phi$, the generators in the Fock-Ivanenko
covariant derivative (6) are

${S^{ab}}=0\;.$ (11)

The spacetime tensor field transforms in the way of the same type of tensor
under local Lorentz transformation. Hence the covariant derivative (6) acts on
spacetime indices and tangent space ones separately according to different
Lorentz representations that the index belongs, e.g.

${D_{\mu}}({\partial_{\nu}}\psi)={\partial_{\mu}}({\partial_{\nu}}\psi)-\frac{i}{2}{\Gamma_{\mu
ab}}{S^{ab}}{\partial_{\nu}}\psi+\frac{i}{2}{\Gamma_{\mu
ab}}{({V^{ab}})_{\nu}}^{\lambda}{\partial_{\lambda}}\psi={\nabla_{\mu}}({\partial_{\nu}}\psi)-\frac{i}{2}{\Gamma_{\mu
ab}}{S^{ab}}{\partial_{\nu}}\psi\;.$ (12)

Utilizing this property of the covariant derivative, one can prove

$[{D_{\mu}},{D_{\nu}}]=-{T_{\mu\nu}}^{\lambda}{D_{\lambda}}-\frac{i}{4}{R_{\mu\nu
ab}}{\sigma^{ab}}\;,$ (13)

where ${T_{\mu\nu}}^{\lambda}$ is the torsion tensor and $R_{\mu\nu ab}$ is
the curvature tensor respectively. The torsion and curvature hence can be
viewed as the gauge field strength of local Lorentz group and local
translation group respectively[25]. Detailed derivation reveals

$\displaystyle\left[{D_{\mu}},{\gamma^{\nu}}\right]$
$\displaystyle=({\nabla_{\mu}}{e^{\nu}}_{c}){\gamma^{c}}-\frac{i}{2}{\Gamma_{\mu
ab}}{e^{\nu}}_{c}[{S^{ab}},{\gamma^{c}}]$ (14)
$\displaystyle=({\nabla_{\mu}}{e^{\nu}}_{c}){\gamma^{c}}+{\Gamma_{\mu
ab}}{e^{\nu a}}{\gamma^{b}}$ $\displaystyle={e^{\nu a}}\left[{({e_{\lambda
a}}{\nabla_{\mu}}{e^{\lambda}}_{c}){\gamma^{c}}+{\Gamma_{\mu
ab}}{\gamma^{b}}}\right]$ $\displaystyle={e^{\nu a}}\left[{{\Gamma_{\mu
ba}}+{\Gamma_{\mu ab}}}\right]{\gamma^{b}}$ $\displaystyle=0\;.$

So for a spinor field $\psi$, we have

${D_{\mu}}{\gamma^{\nu}}\psi={\gamma^{\nu}}{D_{\mu}}\psi$ (15)

or more generally

${D_{\mu}}{\gamma^{\nu}}=0\;.$ (16)

##  3 Fermion scattering cross section by gravitational field

###  3.1 Fermion Lagrangian in curved spacetime background

To investigate the Fermion scattering cross section by gravitational field in
quantum field theory scheme, one needs to write down the Lagrangian of Fermion
in the curved spacetime background. Replacing the partial derivative with
covariant derivative in the Lagrangian of fermion in flat space-time,

$\mathcal{L}_{0}=\bar{\psi}({\partial}-m)\psi\;,$ (17)

one can get the Lagrangian in curved space-time,

$\mathcal{L}=\bar{\psi}(i{\gamma^{\mu}}{D_{\mu}}-m)\psi\;.$ (18)

Observing that the Lagrangian (18) is not self-adjoint,

$\displaystyle{\mathcal{L}^{*}}$ $\displaystyle={\mathcal{L}^{\dagger}}$ (19)
$\displaystyle={\left[{\bar{\psi}\left({i{\gamma^{\mu}}(x){\partial_{\mu}}+\frac{1}{4}{\gamma^{\mu}}(x){\Gamma_{\mu}}-m}\right)\psi}\right]^{\dagger}}$
$\displaystyle=-{\partial_{\mu}}\bar{\psi}i{\gamma^{\mu}}(x)\psi+\frac{1}{4}{\psi^{+}}{\gamma^{0}}{\Gamma_{\mu}}{\gamma^{0}}{\gamma^{0}}{\gamma^{\mu}}(x)\psi-m\bar{\psi}\psi$
$\displaystyle=-{\nabla_{\mu}}\left({\bar{\psi}i{\gamma^{\mu}}(x)\psi}\right)+\bar{\psi}i{\nabla_{\mu}}{\gamma^{\mu}}(x)\psi+\bar{\psi}i{\gamma^{\mu}}(x){\partial_{\mu}}\psi+\frac{1}{4}\bar{\psi}{\Gamma_{\mu}}{\gamma^{\mu}}(x)\psi-m\bar{\psi}\psi$
$\displaystyle=\mathcal{L}-{\nabla_{\mu}}\left({\bar{\psi}i{\gamma^{\mu}}(x)\psi}\right)+\bar{\psi}i[{D_{\mu}},{\gamma^{\mu}}(x)]\psi$
$\displaystyle=\mathcal{L}-{\nabla_{\mu}}\left({\bar{\psi}i{\gamma^{\mu}}(x)\psi}\right)\;,$

where ${\Gamma_{\mu}}={\Gamma_{\mu ab}}{\sigma^{ab}}$. It differs from its
real part only a total divergence term,
${\nabla_{\mu}}\left({\bar{\psi}i{\gamma^{\mu}}(x)\psi}\right)$. $\mathcal{L}$
can be re-written as

$\mathcal{L}=\frac{1}{2}(\mathcal{L}+\mathcal{L}^{*})+\frac{1}{2}{\nabla_{\mu}}\left({\bar{\psi}i{\gamma^{\mu}}(x)\psi}\right)$
(20)

and is indeed a well defined Lagrangian. The corresponding action is

$S=\int{{d^{4}}x\sqrt{-g}\mathcal{L}}\;.$ (21)

We will concentrate on the weak field limit of the static spherically
symmetric gravitational field background. To take account the weak field
limit, we will consider the metric deviation

${h_{\mu\nu}}={g_{\mu\nu}}-{\eta_{\mu\nu}}$ (22)

and spin connection ${\Gamma_{\mu ab}}$ as infinitesimal quantities. Up to the
first order approximation, we have

${{\sqrt{-g}}\cong 1-\frac{1}{2}{h^{\mu}}_{\mu}}$ (23)

and

${e^{\mu}}_{a}\cong{\delta^{\mu}}_{a}-\frac{1}{2}{h^{\mu}}_{a}+{\omega^{\mu}}_{a}\;,$
(24)

where ${\omega^{\mu}}_{a}$, satisfying ${\omega_{ab}}=-{\omega_{ba}}$, are the
pure inertial part of tetrad fields with six degrees of freedom generated by
local Lorentz transformation. One can always remove the inertial
${\omega^{\mu}}_{a}$ by choosing the appropriate local tetrad fields without
losing generality, i.e. choosing the preferred tetrad fields satisfying
${\omega^{\mu}}_{a}=0$. The spacetime indices and the tangent space one play
the same role for the covariant quantities in the first order approximation,
e.g.
${h^{\mu}}_{a}={e^{\mu}}_{b}{h^{b}}_{a}\cong{\delta^{\mu}}_{b}{h^{b}}_{a}$. Up
to the first order approximation, the action (21) can be expanded as

$S=\int{{d^{4}}x\bar{\psi}\left({i\not\partial-m-\frac{1}{2}{\gamma^{a}}{h^{\mu}}_{a}{P_{\mu}}+\frac{1}{4}{\gamma^{a}}{\Gamma_{a}}-\frac{1}{2}{h^{\mu}}_{\mu}(i\not\partial-m)}\right)\psi}\;.$
(25)

So the interaction part between fermion and background gravitational field can
be read from formula (25),

${S_{int}}=\int{{d^{4}}x\bar{\psi}\left({-\frac{1}{2}{\gamma^{a}}{h^{\mu}}_{a}{P_{\mu}}+\frac{1}{4}{\gamma^{a}}{\Gamma_{a}}-\frac{1}{2}{h^{\mu}}_{\mu}(i\not\partial-m)}\right)\psi}\;.$
(26)

Take the assumption that the spacetime is asymptotic Minkowskian and the
initial and final states are both free fermion state in the flat spacetime,
the last term $-\frac{1}{2}{h^{\mu}}_{\mu}\bar{\psi}(i\not\partial-m)\psi$ of
(26) does not contribute to the scattering amplitude actually.

###  3.2 The Scattering of Dirac Fermion

The plane wave expansion of field operator for Dirac fermion is

$\psi=\sum\limits_{s}{\int{{d^{3}}p}}\frac{1}{{\sqrt{EV}}}\left[{{a_{s}}(p){u_{s}}(p){e^{-ipx}}+{b^{\dagger}}_{s}(p){v_{s}}(p){e^{ipx}}}\right]\;,$
(27)

where the normalization condition is set to

${{\bar{u}}_{s}}(p){u_{r}}(p)=m{\delta_{sr}}$ (28)

and

${{\bar{v}}_{s}}(p){v_{r}}(p)=-m{\delta_{sr}}\;.$ (29)

To the first order approximation of spherically symmetric gravitational field,
the Schwarzschild metric can be approximated by

${h_{\mu\nu}}=2\phi(r){\delta^{\mu}}_{\nu}$ (30)

and the spin connection takes the form[26, 27]

${\Gamma_{\mu}}=-\frac{1}{2}{\phi_{,\nu}}{\sigma^{\mu\nu}}\;,$ (31)

where $\phi(r)=-GM/r$ and the imbalanced indices is just a matter of
convention. Denoting the momentum of the particle in the initial state
$\left|i\right\rangle$ and final state $\left|f\right\rangle$ with $\vec{k}$
and $\vec{k}^{\prime}$, the corresponding spin state with $r$ and $s$
respectively, the amplitude for the transition
$\left|i\right\rangle\to\left|f\right\rangle$ is given by

$\left\langle f\right|S\left|i\right\rangle=\left\langle f\right|Texp(-i\int
d^{4}x\mathcal{H}_{I}(x))\left|i\right\rangle\;,$ (32)

where $S$ and $\mathcal{H}_{I}(x)$ are the $S$-matrix and the interacting
Hamiltonian respectively. For simplicity, we take the lowest order non-trivial
expansion of the scattering amplitude, i.e. the first order perturbation, into
consideration,

$\displaystyle M_{Dirac}$ $\displaystyle=i\left\langle
f\right|{S_{{\mathop{\rm int}}}}\left|i\right\rangle$ (33)
$\displaystyle=i\left\langle
f\right|\int{{d^{4}}x\bar{\psi}\left({-\frac{1}{2}{\gamma^{a}}{h^{\mu}}_{a}{P_{\mu}}+\frac{1}{4}{\gamma^{a}}{\Gamma_{a}}}\right)\psi}\left|i\right\rangle$
$\displaystyle=\frac{{2\pi
i}}{{VE}}\delta(E-E^{\prime}){{\bar{u}}_{s}}(k^{\prime})\left({\frac{i}{8}{\gamma^{a}}{q_{\nu}\phi}(q){\sigma^{a\nu}}-\phi(q){\eta_{\mu
a}}{\gamma^{a}}{k_{\mu}}}\right){u_{r}}(k)$ $\displaystyle=-\frac{{2\pi
i}}{{VE}}\delta(E-E^{\prime})\phi(q){{\bar{u}}_{s}}(k^{\prime})\not\tilde{k}{u_{r}}(k)\;,$

where $\vec{q}=\vec{k}^{\prime}-\vec{k}$, the transferred momentum,
$|\vec{q}|=2|\vec{k}|\sin(\theta/2)$ if the scattering angle is $\theta$,
$\phi(q)={4\pi GM}/{|\vec{q}{|^{2}}}$, the three dimensional Fourier component
of $\phi(r)$, $q^{0}=0$ for the energy conservation and
${{\tilde{k}}^{\mu}}=(E,-\vec{k})={{k}_{\mu}}$. To derive the scattering cross
section for unpolarized particles, one needs to make summation over final spin
states and to take average on the initial spin states. The standard derivation
in quantum field theory(e.g. [28]) gives

$\frac{{d\sigma}}{{d\Omega}}=\frac{{{G^{2}}{M^{2}}}}{{4{v^{4}}{{\sin}^{4}}(\theta/2)}}\left[{{{\left({1+{v^{2}}}\right)}^{2}}-{v^{2}}\left({3+{v^{2}}}\right){{\sin}^{2}}\frac{\theta}{2}}\right]\;,$
(34)

where $v$ is velocity of particle($c=1$). Obviously, the cross section (34)
only depends on the velocity of the scattered particle. The non-relativistic
limit of the cross section (34) gives

$\mathop{\lim}\limits_{v\to
0}\frac{{d\sigma}}{{d\Omega}}=\frac{{{G^{2}}{M^{2}}{m^{2}}}}{{16{E_{k}}^{2}{{\sin}^{4}}(\theta/2)}}\;,$
(35)

where ${E_{k}}=m{v^{2}}/2$ is non-relativistic kinetic energy. Formula (35) is
just the familiar form of Rutherford scattering cross section. Considering the
small angle scattering($\sin\theta\cong\theta$), the ultra-relativistic limit
of the cross section (34) reduces to

$\underset{v\to
1}{\mathop{\lim}}\,\frac{d\sigma}{d\Omega}=\frac{16{{G}^{2}}{{M}^{2}}}{{{\theta}^{4}}}\;,$
(36)

which gives the familiar relationship between light deflection angle and
aiming distance $b$,

$\theta=\frac{{4GM}}{b}\;.$ (37)

The fact that the cross section (34) coincides with result in classical theory
in both the relativistic and the non-relativistic limit indicates that the
lowest order expansion of the scattering amplitude (33) can describe the most
important feathers of the fermion scattering by gravitational field.

###  3.3 The Scattering of Majorana Fermion

The plane wave expansion of the field operator for Majorana fermion is[29]

${\psi_{M}}=\sum\limits_{s}{\int{{d^{3}}k}}\frac{1}{{\sqrt{EV}}}\left[{{a_{s}}(k){u_{s}}(k){e^{-ikx}}+{a^{\dagger}}_{s}(k){v_{s}}(k){e^{ikx}}}\right]\;.$
(38)

The spinor wave function ${v_{s}}(k)$ and $u_{s}(k)$ are the charge
conjugation of each other,

${v_{s}}(k)={u^{c}}_{s}(k)=C{\bar{u}^{T}}_{s}(k)\;,$ (39)

where $C=i{\gamma^{0}}{\gamma^{2}}$ is the charge conjugation transformation.
The Majorana field operator is self charge conjugated,

${\psi_{M}}=C{\bar{\psi}_{M}}^{T}={\psi_{M}}^{c}\;,$ (40)

which means Majorana fermion is the antiparticle of itself. There is almost no
difference between the actions of Dirac and Majoranan fermion except for a
factor one half,

$\mathcal{L}_{M}=\frac{1}{2}{\bar{\psi}_{M}}(i{\gamma^{\mu}}{D_{\mu}}-m){\psi_{M}}\;,$
(41)

for Dirac field describes two kinds of particle while only one kind of
particle does Majorana field. Similar derivation of formula (26) gives the
interaction action between Majorana fermion and background gravitational
field,

${S_{int,M}}=\frac{1}{2}\int{{d^{4}}x{{\bar{\psi}}_{M}}\left({-\frac{1}{2}{\gamma^{a}}{h^{\mu}}_{a}{P_{\mu}}+\frac{1}{4}{\gamma^{a}}{\Gamma_{a}}-\frac{1}{2}{h^{\mu}}_{\mu}(i\not\partial-m)}\right){\psi_{M}}}\;.$
(42)

By the same consideration in deriving scattering amplitude for Dirac fermion
(33), the scattering amplitude of Majorana fermion by spherically symmetric
gravitational field in the lowest order expansion of the interaction can be
derived

$\displaystyle{{M}_{Majorana}}=\frac{\pi
i}{VE}\delta\left(E-{{E}^{\prime}}\right)\left[{{{\bar{u}}}_{s}}\left({{k}^{\prime}}\right)\left(\frac{1}{4}{{\gamma}^{a}}{{\Gamma}_{a}}(q)-\frac{1}{2}{{h}^{\mu}}_{a}(q){{\gamma}^{a}}{{k}_{\mu}}\right){{u}_{r}}(k)\right.$
(43)
$\displaystyle\left.-{{{\bar{v}}}_{r}}(k)\left(\frac{1}{4}{{\gamma}^{a}}{{\Gamma}_{a}}(q)+\frac{1}{2}h_{a}^{\mu}(q){{\gamma}^{a}}k_{\mu}^{\prime}\right){{v}_{s}}\left({{k}^{\prime}}\right)\right]\;.$

where ${{\Gamma_{a}}(q)}$ and ${{h^{\mu}}_{a}(q)}$ are the three dimensional
Fourier components of the corresponding ${{\Gamma_{a}}(x)}$ and
${{h^{\mu}}_{a}(x)}$. The first part of $M_{Majorana}$ contains $u$ spinor
only and is the same as in $M_{Dirac}$ while the second part contains $v$
spinor only. Utilizing the charge conjugation relation (40) of spinors, a
similar derivation of (33) gives

$\displaystyle{{\bar{v}}_{r}}(k)\left({\frac{1}{4}{\gamma^{a}}{\Gamma_{a}}(q)+\frac{1}{2}{h^{\mu}}_{a}(q){\gamma^{a}}{{k^{\prime}}_{\mu}}}\right){v_{s}}(k^{\prime})$
(44)
$\displaystyle=\phi(q){{{\bar{v}}}_{r}}(k)\left(\frac{1}{8}\not{q}+{\tilde{\not{k}}}^{\prime}\right){{v}_{s}}({k}^{\prime})$
$\displaystyle=\phi(q){{{\bar{v}}_{r}}(k)\tilde{\not}k^{\prime}{v_{s}}(k^{\prime})}$
$\displaystyle=-\phi(q){{u}_{r}}^{T}(k){{C}^{-1}}{\tilde{\not{k}}}^{\prime}C{{{\bar{u}}}^{T}}_{s}({k}^{\prime})$
$\displaystyle=\phi(q){{u}_{r}}^{T}(k){{({\tilde{\not{k}}}^{\prime})}^{T}}{{{\bar{u}}}^{T}}_{s}({k}^{\prime})$
$\displaystyle=\phi(q){{\bar{u}}_{s}}(k^{\prime})\left({+\tilde{\not}k^{\prime}}\right){u_{r}}(k)$
$\displaystyle=\phi(q){{\bar{u}}_{s}}(k^{\prime})\left({\tilde{\not}k-\not
q}\right){u_{r}}(k)$
$\displaystyle=\phi(q){{\bar{u}}_{s}}(k^{\prime})\tilde{\not}k{u_{r}}(k)\;,$

where ${\tilde{k}}^{\prime\mu}=k^{\prime}_{\mu}$. It leads to
$M_{Majorana}=M_{Dirac}$, which means two kinds of fermions will have the same
scattering behavior in spherically symmetric gravitational field in the lowest
order expansion of the interaction.

The treatment of Dirac and Majorana fermion scattering (33) and (43) can apply
to the neutrino scattering case with the insertion of chiral projection
operator ${P_{L}}=(1-{\gamma_{5}})/2$ for neutrino is only left handed with
positive left hand helicity. Neglecting the small mass of neutrino, the
scattering amplitude of Dirac neutrino is

${M_{Dirac}^{N}}=-\frac{{2\pi
i}}{{VE}}\delta(E-E^{\prime})\phi(q){\bar{u}_{s}}(k^{\prime})\tilde{\not}k{P_{L}}{u_{r}}(k)\;.$
(45)

For Majorana neutrino, the scattering amplitude is

$\displaystyle M_{Majorana}^{N}=\frac{\pi
i}{VE}\delta(E-{E}^{\prime})\left[{{{\bar{u}}}_{s}}({k}^{\prime})\left(\frac{1}{4}{{\gamma}^{a}}{{\Gamma}_{a}}(q)-\frac{1}{2}{{h}^{\mu}}_{a}(q){{\gamma}^{a}}{{k}_{\mu}}\right){{P}_{L}}{{u}_{r}}(k)\right.$
(46)
$\displaystyle\left.-{{{\bar{v}}}_{r}}(k)\left(\frac{1}{4}{{\gamma}^{a}}{{\Gamma}_{a}}(q)+\frac{1}{2}{{h}^{\mu}}_{a}(q){{\gamma}^{a}}{{{{k}^{\prime}}}_{\mu}}\right){{P}_{R}}{{v}_{s}}({k}^{\prime})\right]\;.$

However, a similar derivation as in (44) shows that the relation
$M_{Majorana}^{N}=M_{Dirac}^{N}\equiv M_{N}$ still holds.

Formula (45) shows that the scattering amplitude is non-zero only when the
helicities of the initial and final states of neutrino are the same left
handed one. For neutrino scattering, the summation over final spin states and
average over the initial states are trivial,

$|{{M}_{N}}{{|}^{2}}={{\sum\limits_{spins}{|{{{{M}}}_{N}}|}}^{2}}\;.$ (47)

The spinor calculation gives

$\displaystyle\sum\limits_{spins}{|{{{\bar{u}}}_{s}}({k}^{\prime})\tilde{\not{k}}{{P}_{L}}{{u}_{r}}(k)}{{|}^{2}}$
(48)
$\displaystyle=\sum\limits_{spins}{{{{\bar{u}}}_{s}}({k}^{\prime})\tilde{\not{k}}{{P}_{L}}{{u}_{r}}(k)}{{{\bar{u}}}_{r}}(k)\tilde{\not{k}}{{P}_{L}}{{u}_{s}}({k}^{\prime})$
$\displaystyle=\frac{1}{4}Tr\left[\not{{k}^{\prime}}\tilde{\not{k}}{{P}_{L}}\not{k}\tilde{\not{k}}{{P}_{L}}\right]$
$\displaystyle=\frac{1}{8}Tr\left[\not{{k}^{\prime}}\tilde{\not{k}}\not{k}\tilde{\not{k}}(1-{{\gamma}_{5}})\right]$
$\displaystyle=\frac{1}{2}\left[2\tilde{k}\cdot{k}^{\prime}k\cdot\tilde{k}+i{{\varepsilon}^{\mu\nu\rho\lambda}}{{{{k}^{\prime}}}_{\mu}}{{{\tilde{k}}}_{\nu}}{{k}_{\rho}}{{{\tilde{k}}}_{\lambda}}\right]$
$\displaystyle=4{{E}^{4}}(1-{{\sin}^{2}}\frac{\theta}{2})\;.$

The differential cross section can be given by

$\frac{d{{\sigma}_{N}}}{d\Omega}=\frac{{{G}^{2}}{{M}^{2}}}{{{\sin}^{4}}(\theta/2)}\left[1-{{\sin}^{2}}\frac{\theta}{2}\right]\;,$
(49)

which is just the $v\approx c$ limit in (34).

In conclusion, up to the lowest order expansion in the interaction of neutrino
and gravity, the scattering behavior of neutrino in the spherically symmetric
gravitational field background can not distinguish the type of neutrino from
Dirac one to Majorana one.

###  3.4 The Effect of Torsion on Scattering

Torsion may be an essential part of gravity from a point of view of either
Einstein-Cartan extension of general relativity or the equivalent description
of gravity by both curvature and torsion. To get the torsion involved in the
scattering of fermion by gravity, one needs a detailed investigation of the
interaction between torsion and fermion. The torsion tensor $T^{\mu\nu\rho}$
can be decomposed into three parts according to the irreducible representation
of Lorentz group, i.e. the vector torsion $V^{\mu}$,the axial vector torsion
$A^{\mu}$ and the pure tensor torsion $F^{\mu\nu\rho}$[14],

$\displaystyle{V^{\mu}}$ $\displaystyle={T^{\mu\nu}}_{\nu}$ (50)
$\displaystyle{A^{\mu}}$
$\displaystyle=\frac{1}{{3!}}{\varepsilon^{\mu\nu\rho\lambda}}{T_{\nu\rho\lambda}}$
$\displaystyle{F^{\mu\nu\rho}}$
$\displaystyle={T^{\mu(\nu\rho)}}+\frac{1}{3}{g^{\mu(\nu}}{V^{\rho)}}-\frac{1}{3}{g^{\nu\rho}}{V^{\mu}}\;.$

The pure tensor torsion does not involve in the coupling between torsion and
fermion and the covariant derivative of spinor field can be written as

${\gamma^{\mu}}{D_{\mu}}=\not\partial-\frac{i}{4}{\gamma^{\mu}}{\tilde{\Gamma}_{\mu
ab}}{\sigma^{ab}}+\frac{1}{2}\eta_{1}{\gamma^{\mu}}{V_{\mu}}+\frac{{3i}}{4}{\gamma^{\mu}}{A_{\mu}}{\gamma_{5}}$
(51)

where $\eta_{1}=0$ in the case of minimal coupling of torsion with fermion and
$\eta_{1}\neq 0$ in the non-minimal coupling scheme as well as the framework
of equivalent description of gravity by torsion, e.g. $\eta_{1}=1$ in the
framework of teleparallel gravity. For simplicity, we keep $\eta_{1}$ as te
coupling constant of the vector torsion term in the following of the paper to
account for the latter two cases.

The scattering amplitude of Dirac fermions under the presence of torsion is

${M^{T}_{Dirac}}=\frac{{2\pi
i}}{{VE}}\delta(E-E^{\prime}){\bar{u}_{s}}(k^{\prime})\left({\frac{1}{4}{\gamma^{a}}{{\tilde{\Gamma}}_{a}}(q)+\frac{i}{2}\eta_{1}\not
V(q)-\frac{3}{4}\not
A(q){\gamma_{5}}-\frac{1}{2}{h^{\mu}}_{a}(q){\gamma^{a}}{k_{\mu}}}\right){u_{r}}(k)\;.$
(52)

The vector part contributes a vector current coupling, while the axial vector
part contributes an axial vector current coupling. Similar to the Majorana
fermion scattering amplitude formula (43) by gravity , the scattering
amplitude of the Majorana fermion under the presence of torsion can be written
as

$\begin{array}[]{l}{M^{T}_{Majorana}}=\frac{{\pi
i}}{{VE}}\delta(E-E^{\prime}){{\bar{u}}_{s}}(k^{\prime})\left({\frac{1}{4}{\gamma^{a}}{{\tilde{\Gamma}}_{a}}(q)+\frac{i}{2}\eta_{1}\not
V(q)-\frac{3}{4}\not
A(q){\gamma_{5}}-\frac{1}{2}{h^{\mu}}_{a}(q){\gamma^{a}}{k_{\mu}}}\right){u_{r}}(k)\\\
-\frac{{\pi
i}}{{VE}}\delta(E-E^{\prime}){{\bar{v}}_{r}}(k)\left({\frac{1}{4}{\gamma^{a}}{{\tilde{\Gamma}}_{a}}(q)+\frac{i}{2}\eta_{1}\not
V(q)-\frac{3}{4}\not
A(q){\gamma_{5}}+\frac{1}{2}{h^{\mu}}_{a}(q){\gamma^{a}}{{k^{\prime}}_{\mu}}}\right){v_{s}}(k^{\prime})\;.\end{array}$
(53)

The first term in (53) is one half of the amplitude of Dirac fermion
scattering (52) again just as torsion free case in (43). The effect of torsion
can be separated from the full expressions (52) and (53). As in the torsion
free derivation (44), the second part of torsion terms from the Majorana
particle scattering (53) can be converted to

${\bar{v}_{r}}(k)\left({\frac{i}{2}\eta_{1}\not V(q)-\frac{3}{4}\not
A(q){\gamma_{5}}}\right){v_{s}}(k^{\prime})={\bar{u}_{s}}(k^{\prime})\left({\frac{i}{2}\eta_{1}\not
V(q)+\frac{3}{4}\not A(q){\gamma_{5}}}\right){u_{r}}(k)\;,$ (54)

which reduces the full expression (53) to

$\begin{array}[]{l}{M^{T}_{Majorana}}=\frac{{\pi
i}}{{VE}}\delta(E-E^{\prime}){{\bar{u}}_{s}}(k^{\prime})\left({\frac{1}{4}{\gamma^{a}}{{\tilde{\Gamma}}_{a}}(q)-\frac{1}{2}{h^{\mu}}_{a}(q){\gamma^{a}}{k_{\mu}}}\right){u_{r}}(k)\\\
-\frac{{\pi
i}}{{VE}}\delta(E-E^{\prime}){{\bar{v}}_{r}}(k)\left({\frac{1}{4}{\gamma^{a}}{{\tilde{\Gamma}}_{a}}(q)+\frac{1}{2}{h^{\mu}}_{a}(q){\gamma^{a}}{{k^{\prime}}_{\mu}}}\right){v_{s}}(k^{\prime})\\\
+\frac{{2\pi
i}}{{VE}}\delta(E-E^{\prime}){{\bar{u}}_{s}}(k^{\prime})\left({-\frac{3}{4}\not
A(q){\gamma_{5}}}\right){u_{r}}(k)\;.\end{array}$ (55)

The expression (55) reveals that the vector torsion does not play any role in
the scattering amplitude of the Majorana fermion, while the axial vector
torsion contributes the same term in Majorana case just as in the Dirac
fermion case up to the first order of perturbation. The result is independent
of the spacetime background actually although it is derived in a weak field
approximation of the spherically symmetric background. The conclusion may be a
reflection of the fact that the Majorana fermion has neither a $U(1)$ symmetry
hence nor a conserved $U(1)$ current which can couple to the vector torsion to
contribute to the scattering amplitude. The possible different response to the
vector torsion between Majorana and Dirac particle can distinguish the fermion
type of neutrino from the scattering behavior by vector torsion. It can also
serve as a prediction of the non-minimal coupling between torsion and fermion
or a verification of the equivalence between the curvature and the torsion
descriptions of gravity.

##  4 Summary and Outlook

The conclusion that Dirac and Majorana fermion are indistinguishable by the
torsionless spherically symmetric gravitational field can not be regarded as a
counter proof of the conclusion of Papini et.al. in literature [6] because the
Lense-Thirring background considered in [6] is not spherically symmetric and
is generated by the rotating mass with angular momentum. The spin-gravity
interaction is generated by the gravitomagenetic frame dragging effect of the
rotating source. On the other hand, our conclusion is only based on the first
order perturbation of the interaction between gravitational field and the
fermion field. Results beyond first order perturbation and in the case of
gravitational field without spherical symmetry are expected in the subsequent
works. However, the response difference on vector torsion coupling between
Dirac in (52) and Majorana fields in (55) is expected to hold beyond the
lowest order perturbation.

Torsion can be generated from the polarized distribution of spinor field which
may happen in some high density astronomical objects such as neutron star or
near the blackhole[6, 25]. The possible detection of torsion-spinor
interaction can be the scattering behavior of neutrino beam passing by these
kind of objects.

In view point of effective field theory, the connection may be an active
degree of freedom, which can be excited at high energy scale[11, 30]. In low
energy effective theory of gravity, the torsion free Levi-Civita connection
may be a confined or condensed phase of the high energy connection degree of
freedom. There remains the possibility that the low energy connection phase
deviates from Levi-Civita one at the large scale as in the large scale Lorentz
violation scheme[16, 17, 18, 19]. The large scale Lorentz violation approach
to the gravitational theory at cosmic scale and cosmology predicts a torsion
distribution at the cosmic scale[20, 21, 22, 23, 24]. The prediction and
detection on the behavior of neutrino propagation at cosmic scale can verify
the predicted cosmic torsion by the large scale Lorentz violation model.
However the technique developed in this paper needs modification and extension
to the cosmological torsion case. The geometro-hydrodynamics approach may be
employed to the prediction of neutrino propagation path in the FRW background
with cosmic torsion distribution rather than the scattering of free neutrino
by a localized torsion and curvature distribution[31, 32, 15].

In the framework of equivalent description of gravity by torsion rather than
curvature or both of them, the coupling of vector torsion with fermion field
does not vanish[14, 12, 13]. Its correspondence in the Hilbert——Einstein
framework needs to be clarified. The possible relation between different
response to the vector torsion by Dirac and Majorana fermion in the frame work
of teleparallel gravity and spin-gravity coupling induced different spin-flip
matrix element between Dirac and Majorana fermion in a rotating gravitational
field background[6, 10, 9] needs to be investigated further.

## Acknowledgment

This work is supported by the National Natural Science Foundation of China,
under Grant No. 11775080 and Grant No. 11865016.

## References

  * [1] E. Majorana, “Teoria simmetrica dell’elettrone e del positrone,” Nuovo Cim. 14 (1937) 171–184.
  * [2] W. H. Furry, “On transition probabilities in double beta-disintegration,” Phys. Rev. 56 (1939) 1184–1193.
  * [3] L. Oberauer, A. Ianni, and A. Serenelli, Solar Neutrino Physics: The Interplay between Particle Physics and Astronomy.  Wiley-VCH, 2020.
  * [4] K. L. Ng, “Gravitational form-factors of the neutrino,” Phys. Rev. D 47 (1993) 5187–5190, arXiv:gr-qc/9305002.
  * [5] K. L. Ng, “Equivalence principle and discrete symmetries of Dirac and Majorana fermions in gravitational field,” Nuovo Cim. B 109 (1994) 1143–1146.
  * [6] D. Singh, N. Mobed, and G. Papini, “Can gravity distinguish between Dirac and Majorana neutrinos?,” Phys. Rev. Lett. 97 (2006) 041101, arXiv:gr-qc/0605153.
  * [7] A. Menon and A. M. Thalapillil, “Interaction of Dirac and Majorana Neutrinos with Weak Gravitational Fields,” Phys. Rev. D 78 (2008) 113003, arXiv:0804.3833 [hep-ph].
  * [8] S. A. Alavi and A. Abbasnezhad, “Can gravity distinguish between Dirac and Majorana Neutrinos?,” Grav. Cosmol. 22 no. 3, (2016) 288–298, arXiv:1201.4741 [hep-ph].
  * [9] D. Singh, N. Mobed, and G. Papini, “Reply to comment on ‘Can gravity distinguish between Dirac and Majorana neutrinos?’,” Phys. Rev. Lett. 98 (2007) 069002, arXiv:gr-qc/0611016.
  * [10] J. F. Nieves and P. B. Pal, “Comment on ‘Can gravity distinguish between Dirac and Majorana neutrinos?’,” Phys. Rev. Lett. 98 (2007) 069001, arXiv:gr-qc/0610098.
  * [11] J. F. Donoghue, “Is the spin connection confined or condensed?,” Phys. Rev. D 96 no. 4, (2017) 044003, arXiv:1609.03523 [hep-th].
  * [12] H. I. Arcos, V. C. de Andrade, and J. G. Pereira, “Torsion and gravitation: A New view,” Int. J. Mod. Phys. D 13 (2004) 807–818, arXiv:gr-qc/0403074.
  * [13] R. T. Hammond, “Torsion gravity,” Rept. Prog. Phys. 65 (2002) 599–649.
  * [14] R. Aldrovandi and J. G. Pereira, Teleparallel Gravity: An Introduction.  Springer, 2013.
  * [15] I. L. Shapiro, “Physical aspects of the space-time torsion,” Phys. Rept. 357 (2002) 113, arXiv:hep-th/0103093.
  * [16] Y. Wu, X. Xue, L. Yang, and T.-C. Yuan, “The Effective Gravitational Theory at Large Scale with Lorentz Violation,” (2015) , arXiv:1510.00814 [gr-qc].
  * [17] Y. Wu and X. Xue, “${SIM(2)}$ gravity gauge theory,” Journal of East China Normal University (Natural Science) 2016 no. 3, (2016) 76–83.
  * [18] L. Yang, Y. Wu, W. Wei, X. Xue, and T.-C. Yuan, “The effective gravitation theory at large scale with Lorentz violation,” Chinese Science Bulletin no. 9, (2017) 944–950.
  * [19] W. Wei, J. Shen, Y. Wu, L. Yang, X. Xue, and T.-C. Yuan, “${E(2)}$ gauge theory model of effective gravitational theory at large scale,” Acta. Phys. Sin. 66 no. 13, (2017) 130301–130301.
  * [20] J. Shen and X. Xue, “Large Scale Lorentz Violation Gravity and Dark Energy,” in Proceedings, 28th International Symposium on Lepton Photon Interactions at High Energies (LP17): Guangzhou (Guangdong), China, August 7-12, 2017, W. Wang and Z.-z. Xing, eds., pp. 459–475.  2020\.  arXiv:1802.03502 [gr-qc].
  * [21] H. Zhai, J. Shen, and X. Xue, “Effective quintessence from string landscape,” Acta Phys. Sin. 68 no. 13, (2019) 139501.
  * [22] H. Zhai, J. Shen, and X. Xue, “Uplifting of AdS type to quintessence-like potential induced by frozen large-scale Lorentz violation,” Chin. Phys. C 44 no. 8, (2020) 085101, arXiv:1906.11860 [hep-th].
  * [23] Q. Li, J. Li, Y. Zhou, and X. Xue, “The Effective Potential Originating from Swampland and the Non-trivial Brans-Dicke Coupling,” Chin. Phys. C 44 no. 10, (2020) 105108, arXiv:2003.09121 [gr-qc].
  * [24] J. Li, Y. Zhou, and X. Xue, “Spatial Curvature and Large Scale Lorentz Violation,” arXiv:2112.02364 [gr-qc].
  * [25] F. W. Hehl, P. Von Der Heyde, G. D. Kerlick, and J. M. Nester, “General Relativity with Spin and Torsion: Foundations and Prospects,” Rev. Mod. Phys. 48 (1976) 393–416.
  * [26] G. Lambiase, G. Papini, R. Punzi, and G. Scarpetta, “Neutrino optics and oscillations in gravitational fields,” Physics 71 no. 7, (2005) 118–120.
  * [27] G. Lambiase, G. Papini, R. Punzi, and G. Scarpetta, “Neutrino optics and oscillations in gravitational fields,” Phys. Rev. D 71 (2005) 073011, arXiv:gr-qc/0503027.
  * [28] V. P. Nair, Quantum field theory: A modern perspective.  2005\.
  * [29] M. Srednicki, “Quantum Field Theory: Spin One Half,” Quantum Field Theory 231 no. 1186, (2004) 664.
  * [30] J. F. Donoghue, “Conformal model of gravitons,” Phys. Rev. D96 no. 4, (2017) 044006, arXiv:1609.03524 [hep-th].
  * [31] M. Trukhanova, Iv., “The geometro-hydrodynamical formalism of the quantum spinning particle,” PTEP 2018 no. 12, (2018) 123A01, arXiv:1810.04716 [gr-qc].
  * [32] M. I. Trukhanova, “The Geometro-Hydrodynamical Representation of the Torsion Field,” Phys. Lett. A 381 (2017) 2887–2892, arXiv:1702.02168 [gr-qc].