# Effects of gravitational lensing on neutrino oscillation in
$\gamma$-spacetime

Hrishikesh Chakrabarty<EMAIL_ADDRESS>School of Astronomy and
Space Sciences, University of Chinese Academy of Sciences, Beijing, China
Debasish Borah<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology Guwahati, Assam 781039, India  Ahmadjon Abdujabbarov
<EMAIL_ADDRESS>Shanghai Astronomical Observatory, 80 Nandan Road, Shanghai
200030, P. R. China  Ulugh Beg Astronomical Institute, Astronomy St 33,
Tashkent 100052, Uzbekistan  National University of Uzbekistan, Tashkent
100174, Uzbekistan  Tashkent Institute of Irrigation and Agricultural
Mechanization Engineers,39 Kori Niyoziy st., Tashkent 100000, Uzbekistan
Daniele Malafarina<EMAIL_ADDRESS>Department of Physics,
Nazarbayev University, 53 Kabanbay Batyr avenue, 010000 Astana, Kazakhstan
Bobomurat Ahmedov<EMAIL_ADDRESS>Ulugh Beg Astronomical Institute,
Astronomy St 33, Tashkent 100052, Uzbekistan  National University of
Uzbekistan, Tashkent 100174, Uzbekistan  Tashkent Institute of Irrigation and
Agricultural Mechanization Engineers,39 Kori Niyoziy st., Tashkent 100000,
Uzbekistan

###### Abstract

We study the effects of gravitational lensing on neutrino oscillations in the
$\gamma$-spacetime which describes a static, axially-symmetric and
asymptotically flat solution of the Einstein’s field equations in vacuum.
Using the quantum-mechanical treatment for relativistic neutrinos, we
calculate the phase of neutrino oscillations in this spacetime by considering
both radial and non-radial propagation. We show the dependence of the
oscillation probability on the absolute neutrino masses, which in the two-
flavour case also depends upon the sign of mass squared difference, in sharp
contrast with the well-known results of vacuum oscillation in flat spacetime.
We also show the effects of the deformation parameter $\gamma$ on neutrino
oscillations and reproduce previously known results for the Schwarzschild
metric. We then extend these to a more realistic three flavours neutrino
scenario and study the effects of the parameter $\gamma$ and the lightest
neutrino mass while using best fit values of neutrino oscillation parameters.

##  I Introduction

The fact that neutrinos have tiny masses and large mixing has been well-
established by now, thanks to several experimental evidences gathered over a
long periods of time Zyla et al. (2020). For a latest update on neutrino
parameter values, one may have a look at the latest global fit of neutrino
data Esteban et al. (2020); de Salas et al. (2021). In spite of this great
success story, there are still several unanswered questions related to
neutrinos. In addition to the fundamental question of the origin of light
neutrino mass, we do not know a few things from experimental point of view.
For example, we still do not know whether the atmospheric mixing angle lies in
the first or second octant, even though the best fit value prefers the second
octant. Also, since only two mass squared differences are experimentally
measured, we do not have any idea about the lightest neutrino mass leading to
the possibility of both normal ordering (NO): $m_{1}<m_{2}<m_{3}$ as well as
inverted ordering (IO): $m_{3}<m_{1}<m_{2}$. The global fit data, however,
suggests a normal mass ordering. While neutrino oscillation experiments are
sensitive to mass squared differences only, there are other ways to measure
the absolute neutrino mass scale. For example, the tritium beta decay
experiment KATRIN has recently measured an upper limit of 1.1 eV on absolute
neutrino mass scale Aker et al. (2019) which is however, weaker than the upper
limits from cosmology. Latest data from Planck collaboration constrains the
sum of absolute neutrino masses as $\sum_{i}\lvert m_{i}\rvert<0.12$ eV
Aghanim et al. (2018). In spite of recent advances in determining the leptonic
Dirac CP phase Abe et al. (2020) suggesting a maximal CP violation, the
measurement of leptonic Dirac CP phase is not statistically significant enough
to be considered as a discovery.

The results of various neutrino oscillation experiments are interpreted by
using the quantum mechanical neutrino oscillation probability derived in flat
Minkowski spacetime. While such descriptions are valid for the laboratory
based oscillation experiments, in extreme astrophysical environments or
cosmology in general, the effects of non-trivial spacetime background may be
important. This has led to several theoretical studies on the effects of
curved spacetime on neutrino oscillations Wudka (1991); Grossman and Lipkin
(1997); Cardall and Fuller (1997); Piriz et al. (1996); Fornengo et al.
(1997); Ahluwalia and Burgard (1996); Bhattacharya et al. (1999); Pereira and
Zhang (2000); Crocker et al. (2004); Lambiase et al. (2005); Godunov and
Pastukhov (2011); Ren and Zhang (2010); Geralico and Luongo (2012);
Chakraborty (2014); Visinelli (2015); Chakraborty (2015); Zhang and Li (2016);
Alexandre and Clough (2018); Blasone et al. (2020); Buoninfante et al. (2020);
Boshkayev et al. (2020); Mandal (2021); Koutsoumbas and Metaxas (2020); Swami
et al. (2020); Capolupo et al. (2020). The general conclusion reached in most
of these works is the increase in effective neutrino oscillation length due to
gravitational redshift in such curved spacetime. In this work, we consider the
effects of gravitational lensing on neutrino oscillations. In the presence of
lensing around a massive astrophysical object, the oscillation probability at
a particular point is calculated by considering all possible trajectories of
neutrinos which get focused at that point due to the lensing phenomenon. While
this was first discussed in Fornengo et al. (1997), some follow-up works have
also appeared in the literature Crocker et al. (2004); Alexandre and Clough
(2018); Swami et al. (2020). While these works have focused on the
Schwarzschild spacetime for such studies, we consider here a different and
well-motivated spacetime metric.

In this paper, we study gravitational lensing of neutrinos in the so-called
$\gamma$-metric, also known as the ‘Zipoy-Voorhees’ spacetime Zipoy (1966);
Voorhees (1970). The $\gamma$-metric is a static, axially-symmetric and
asymptotically flat solution of Einstein’s equations in vacuum and it belongs
to the Weyl class of exact solutions. The metric is characterized by only two
parameters: $M$ related to the mass of the source and $\gamma$ related to the
static mass quadrupole moment of the source. The parameter $\gamma$ describes
deviations from the spherical symmetry in such a way that $\gamma=1$
reproduces the Schwarzschild metric while $\gamma>1$ ($\gamma<1$) describes an
oblate (prolate) source. The structure of the metric at small radii is also
very interesting since when $\gamma\neq 0$ the horizon at $2M$ is replaced by
a true curvature singularity. Therefore, the singularity in the
$\gamma$-metric is naked and Virbhadra in Virbhadra (1996) showed that the
singularity has ‘directional nakedness’ which depends on the angle of approach
Kodama and Hikida (2003). For $\gamma<1$, this singularity is globally visible
along polar as well as equatorial direction. For $\gamma>1$, the singularity
is globally naked along equatorial direction, but it is not visible (even
locally) along the polar direction. In Chakrabarty et al. (2018), it was shown
that the same singularity at $2M$ can be resolved in conformal gravity
theories. Also it is worth noticing that, the ‘Zipoy-Voorhees’ metric is non-
integrable, as shown in Lukes-Gerakopoulos (2012), which leads to interesting
and possibly chaotic behavior for the general motion of test particles.

Goedesic motion of massive and massless particles in the $\gamma$-metric was
studied in Herrera et al. (1999); Chowdhury et al. (2012); Boshkayev et al.
(2016); Toshmatov and Malafarina (2019). Shadows and optical properties was
calculated in Abdikamalov et al. (2019) and it was shown that this spacetime
can behave as a black hole mimicker for some allowed values of $\gamma$.
Therefore, the $\gamma$-metric provides an excellent and simple candidate to
study toy models of astrophysical scenarios where the exterior region of a
massive compact object is not given by a black hole line element. We calculate
the quantum mechanical phase of neutrino oscillations in this metric by
considering both radial and non-radial propagation of neutrinos. Using the
weak-field approximation, we also arrive at simple expressions for this phase
and use it to calculate the neutrino oscillation probabilities. We then use
the effects of gravitational lensing and calculate the oscillation
probabilities for both two-flavor and three-flavor regimes.

The article is structured as follows: In Sec. II, we give a very brief review
of neutrino oscillation in flat and curved spacetimes and in Sec. III, we
present a brief overview of $\gamma$-metric. In Sec. IV, we calculate the
phase for neutrino oscillations in the $\gamma$-metric for radial and non-
radial propagation of neutrinos and in Sec. V, we calculate the oscillation
probability of neutrinos moving in a geometry described by $\gamma$-metric. In
Sec. VI, we discuss gravitational lensing of neutrinos and present the
numerical results of both two and three flavor neutrino oscillation
probability in this spacetime. Finally, in Sec. VII we summarize and comment
on our results. Throughout the paper we employ $(-,+,+,+)$ signature for the
line element and adopt geometrical units setting $c=G=\hbar=1$.

##  II Neutrino oscillations

Neutrinos are produced and detected in different flavor eigenstates denoted by
$\ket{\nu_{\alpha}}$, and the flavor eigenstates are superposition of mass
eigenstates denoted by $\ket{\nu_{i}}$. So a flavor eigenstate can be written
in terms of mass eigenstates as

$\ket{\nu_{\alpha}}=\sum U^{*}_{\alpha i}\ket{\nu_{i}},$ (1)

where $\alpha=e,\mu,\tau$ and $i=1,2,3$. $U$ is a $3\times 3$ unitary mixing
matrix. For three flavor neutrino oscillation, this is known as the
Pontecorvo-Maki-Nakagawa-Sakata (PMNS) leptonic mixing matrix Pontecorvo
(1957); Maki et al. (1962); Pontecorvo (1967).

We assume that the neutrino wave-function is a plane wave as considered
originally in Pontecorvo (1957); Maki et al. (1962); Pontecorvo (1967) and it
propagates from a source S located at a spacetime event $(t_{S},x_{S})$ to a
detector D located at a spacetime event $(t_{D},x_{D})$. So the wave-function
at the detector point is given by

$\ket{\nu_{i}(t_{D},x_{D})}=\exp\left(-i\Phi_{i}\right)\ket{\nu_{i}(t_{S},x_{S})}\
.$ (2)

Neutrinos are expected to be produced initially in the flavour eigenstate
$\ket{\nu_{\alpha}}$ at S and then travel to the detector D. In that case, the
probability of the change in neutrino flavour from $\nu_{\alpha}$ to
$\nu_{\beta}$ at D is given by

$\displaystyle\mathcal{P}_{\alpha\beta}$
$\displaystyle=|\braket{\nu_{\beta}}{\nu_{\alpha}(t_{D},x_{D})}|^{2}=$ (3)
$\displaystyle=\sum_{i,j}U_{\beta i}U_{\beta j}^{*}U_{\alpha j}U_{\beta
i}^{*}\exp\left[-i(\Phi_{i}-\Phi_{j})\right].$

In flat spacetime, the phase $\Phi_{i}$ is given by

$\Phi_{i}=E_{i}(t_{D}-t_{S})-p_{i}(\vec{x}_{D}-\vec{x}_{S}).$ (4)

It is typically assumed that all the mass eigenstates in a flavour eigenstate
initially produced at the source have equal momentum or energy Bilenky and
Pontecorvo (1978). Either of these assumptions together with
$(t_{D}-t_{S})\simeq|(\vec{x}_{D}-\vec{x}_{S})|$ for relativistic neutrinos
($E_{i}\gg m_{i}$) lead to

$\Delta\Phi_{ij}\equiv\Phi_{i}-\Phi_{j}\simeq\frac{\Delta
m_{ij}^{2}}{2E_{0}}(\vec{x}_{D}-\vec{x}_{S}),$ (5)

where $\Delta m_{ij}^{2}=m_{i}^{2}-m_{j}^{2}$, and $E_{0}$ is the average
energy of the relativistic neutrinos produced at the source.

To generalize the expression of the phase $\Phi_{i}$ for neutrino propagation
in curved spacetime, Eq. (4) is written in the covariant form Fornengo et al.
(1997)

$\Phi_{i}=\int^{D}_{S}p_{\mu}^{(i)}dx^{\mu},$ (6)

where

$p_{\mu}^{(i)}=m_{i}g_{\mu\nu}\frac{dx^{\mu}}{ds}$ (7)

is the canonical conjugate momentum to the coordinates $x^{\mu}$ and
$g_{\mu\nu}$ and $ds$ are the metric tensor and line element of the curved
spacetime, respectively.

##  III  The $\gamma$-metric

The $\gamma$-metric is an exact solution of the Einstein’s field equations in
vacuum belonging to the class of Weyl metrics. The line element for
$\gamma$-metric can be written as Zipoy (1966); Voorhees (1970)

$\displaystyle ds^{2}=$
$\displaystyle-f^{\gamma}dt^{2}+f^{\gamma^{2}-\gamma}g^{1-\gamma^{2}}\left(\frac{dr^{2}}{f}+r^{2}d\theta^{2}\right)+$
(8) $\displaystyle+f^{1-\gamma}r^{2}\sin^{2}\theta d\phi^{2},$

where

$\displaystyle f=1-\frac{2M}{r},$ (9) $\displaystyle
g=1-\frac{2M}{r}+\frac{M^{2}\sin^{2}\theta}{r^{2}}.$

There are two parameters characterizing this line element: $M>0$ is related to
the mass of the source, and $\gamma>0$ quantifies the departure of the
geometry from spherical symmetry. For $\gamma=1$, the spacetime becomes
spherically symmetric and from Eq. 8 we recover the Schwarzschild solution in
Schwarzschild coordinates. For $\gamma>1$ ($\gamma<1$), the spacetime posses
oblate (prolate) structure, which may be used to describe the exterior of a
compact object. Interiors for the $\gamma$-spacetime were obtained in Stewart
et al. (1982); Herrera et al. (2005); Hernandez-Pastora et al. (2016). In
order to have a better understanding of the interpretation of the parameters,
it is useful to consider the asymptotic expansion in multipoles of the
gravitational potential. We see that the total mass as measured by an observer
at infinity $M_{\rm tot}$, namely the monopole moment, and the quardupole
moment $Q$ are given by Hernández-Pastora and Martín (1994)

$\displaystyle M_{\rm tot}$ $\displaystyle=M_{\rm ADM}=\gamma M\ ,$ (10)
$\displaystyle Q$ $\displaystyle=\frac{\gamma M^{3}(1-\gamma^{2})}{3}.$

For $\gamma\neq 1$, the spacetime has a naked singularity at the radial
coordinate $r_{\rm sing}=2M$ as it can be seen by calculating the Kretschmann
invariant. So the spacetime is geodesically incomplete with the radial
coordinate range $r\in(2M,\infty)$. In the following discussions, we consider
the exterior of a massive object to be approximated by the $\gamma$-metric in
the mentioned radial coordinate range.

##  IV Phase of neutrino oscillation

Let us now calculate the phase of neutrino oscillation probability in the
$\gamma$-spacetime. First, for simplicity, we shall rename the metric
coefficients as follows.

$\displaystyle g_{tt}$ $\displaystyle=-A,\ \ \ \ g^{tt}=-1/A,$ (11)
$\displaystyle g_{rr}$ $\displaystyle=B,\ \ \ \ g^{rr}=1/B,$ $\displaystyle
g_{\theta\theta}$ $\displaystyle=C,\ \ \ \ g^{\theta\theta}=1/C,$
$\displaystyle g_{\phi\phi}$ $\displaystyle=D,\ \ \ \ g_{\phi\phi}=1/D.$

We can write the components of canonical momenta $p_{\mu}^{(k)}$ for test
particles moving on the equatorial plane $\theta=\pi/2$ in the
$\gamma$-spacetime as

$\displaystyle p_{t}^{(k)}=m_{k}g_{tt}\frac{dt}{ds},$ (12) $\displaystyle
p_{r}^{(k)}=m_{k}g_{rr}\frac{dr}{ds},$ $\displaystyle
p_{\phi}^{(k)}=m_{k}g_{\phi\phi}\frac{d\phi}{ds},$

where $m_{k}$ is the particle’s mass, i.e. the neutrino mass eigenstate, ad
$s$ is the trajectory’s affine parameter. These momenta are related to each
other and the mass of the $k$-th eigenstate by the mass-shell relation

$-m_{k}^{2}=g^{tt}p_{t}^{2}+g^{rr}p_{r}^{2}+g^{\phi\phi}p_{\phi}^{2}.$ (13)

To avoid confusion, we have dropped the prefix ‘$k$’ in the momenta.

###  IV.1 Radial case

We first investigate the case of radial propagation of neutrinos. For radial
propagation, $d\phi=0$. Therefore, from Eq. (12), we have

$\frac{dt}{ds}=\frac{p_{t}}{m_{k}g_{tt}},\ \ \ \
\frac{dr}{ds}=\frac{p_{r}}{m_{k}g_{rr}}.$ (14)

Thanks to the spacetime being static we know that the energy $E$ of test
particles in conserved and we can define, $p_{t}=-E_{k}$ and $p_{r}=p_{k}(r)$,
where the subscript $k$ refers to the neutrino flavor. With this notation, the
phase of a neutrino propagating radially in a light-ray trajectory becomes

$\Phi_{k}=\int^{D}_{S}\left[-E_{k}\left(\frac{dt}{dr}\right)_{0}+p_{k}(r)\right]dr,$
(15)

where, $S$ and $D$ denote the source of the neutrino and detector,
respectively. We again emphasize that the phase in Eq. (15) is not the phase
on a classical trajectory for the mass eigenstates but the phase calculated on
a light-ray trajectory Fornengo et al. (1997). Now, from Eq. (14), we get the
light-ray differential as

$\left(\frac{dt}{dr}\right)_{0}=\frac{E_{0}}{p_{0}(r)}\frac{B}{A},$ (16)

where $E_{0}$ is the energy of a massless particle at infinity. The mass-shell
relation gives

$\displaystyle p_{0}(r)=\pm E_{0}\sqrt{\frac{B}{A}},$ (17) $\displaystyle
p_{k}(r)=\pm\sqrt{\frac{BE_{k}^{2}}{A}-Bm_{k}^{2}}.$

Using Eq. (16) and (17) in the expression for phase, we obtain

$\displaystyle\Phi_{k}$
$\displaystyle=\pm\int^{D}_{S}E_{k}\sqrt{\frac{B}{A}}\left[-1+\sqrt{1-\frac{m_{k}^{2}A}{E_{k}^{2}}}\right]dr.$
(18)

Now, using the fact that $0<A<1$, we expand the square-root under the bracket
to write

$\Phi_{k}=\pm\int^{D}_{S}\sqrt{AB}E_{k}\frac{m_{k}^{2}}{2E_{k}^{2}}\ dr.$ (19)

In relativistic approximation ($m_{k}<<E_{k}$), the following relation holds
Fornengo et al. (1997)

$\displaystyle E_{k}$ $\displaystyle\simeq
E_{0}+\mathcal{O}\left(\frac{m_{k}^{2}}{2E_{0}}\right)\ ,$ (20) $\displaystyle
E_{k}\frac{m_{k}^{2}}{2E_{k}^{2}}$ $\displaystyle\simeq
E_{0}\frac{m_{k}^{2}}{2E_{0}^{2}}.$

Hence the phase becomes

$\Phi_{k}=\pm\frac{m_{k}^{2}}{2E_{0}}\int^{D}_{S}\sqrt{AB}\ dr.$ (21)

Integrating the above integral from $r_{S}$ to $r_{D}$, we get

$\Phi_{k}\simeq\pm\frac{m_{k}^{2}}{2E_{0}}(r_{D}-r_{S})\left[1-\frac{M^{2}(\gamma^{2}-1)}{2r_{S}r_{D}}\right].$
(22)

Here the phase is calculated up to order $\mathcal{O}\left(M^{2}/r^{2}\right)$
to obtain some information on its dependence on $\gamma$. In fact, we can see
that $\gamma$ has no effect in the order $\mathcal{O}\left(M/r\right)$. So in
the lowest order in ($M/r$), the phase is same as that in Schwarzschild
spacetime. This behaviour is interesting as the phase of the radially
propagating neutrinos in the equatorial plane and in the weak field limit is
not affected by the deformed geometry of the spacetime. The reason behind this
is that, in the weak field limit $g_{tt}g_{rr}\simeq 1$ for the
$\gamma$-metric.

###  IV.2 Non-radial case

Now we would like to concentrate on non-radial motion of neutrinos. Let us
write the phase for a neutrino travelling along a light ray using Eq. (12) and
Eq. (17) as

$\Phi_{k}=\int^{D}_{S}\left[-E_{k}\left(\frac{dt}{dr}\right)_{0}+p_{r}+J_{k}\left(\frac{d\phi}{dr}\right)_{0}\right]dr,$
(23)

where the integral is taken along the light-ray trajectory. For simplicity, we
have ignored the indices and used conservation of energy and angular momentum
which give $p_{t}^{(k)}=-E_{k}$, $p_{\phi}^{(k)}=J_{k}$, where $J$ is the
angular momentum of the test particle. Now from the equations for canonical
momenta we can obtain

$\frac{dt}{dr}=\frac{-E_{k}}{p_{r}}\frac{B}{(-A)},\ \ \ \
\frac{d\phi}{dr}=\frac{J_{k}}{p_{r}}\frac{B}{D}.$ (24)

And along the light-ray paths, these equations become

$\left(\frac{dt}{dr}\right)_{0}=\frac{E_{0}}{p_{0}(r)}\frac{B}{A},\ \ \ \
\left(\frac{d\phi}{dr}\right)_{0}=\frac{J_{0}}{p_{0}(r)}\frac{B}{D}.$ (25)

It is convenient to express the angular momentum $J_{k}$ as a function of the
energy $E_{k}$, the impact parameter $b$ and the velocity at infinity
$v_{k}^{(\infty)}$ as

$J_{k}=E_{k}bv_{k}^{(\infty)}.$ (26)

Asymptotically, the metric is flat. So we can write

$\displaystyle v_{k}^{(\infty)}$
$\displaystyle=\frac{\sqrt{E_{k}^{2}-m_{k}^{2}}}{E_{k}}\simeq
1-\frac{m_{k}^{2}}{2E_{k}^{2}},$ (27) $\displaystyle J_{k}$
$\displaystyle\simeq E_{k}b\left(1-\frac{m_{k}^{2}}{2E_{k}^{2}}\right)\ ,$

where in the last equality we used the relativistic approximation. For a
massless particle, the angular momentum is

$J_{0}=E_{0}b.$ (28)

Using Eqs. (25), (26), (27) and (28), the phase can be written as

$\Phi_{k}=\int^{d}_{S}\frac{E_{0}E_{k}B}{p_{0}(r)}\left[-\frac{1}{A}+\frac{p_{0}p_{k}(r)}{E_{0}E_{k}B}+\frac{b^{2}}{D}\left(1-\frac{m_{k}^{2}}{2E_{k}^{2}}\right)\right].$
(29)

Now from the mass-shell relation, we can find the following two relations

$\displaystyle\frac{p_{0}(r)}{E_{0}B}$
$\displaystyle=\pm\sqrt{\frac{1}{AB}-\frac{b^{2}}{BD}},$ (30)
$\displaystyle\frac{p_{0}(r)p_{k}(r)}{E_{0}E_{k}B}$
$\displaystyle=\frac{1}{A}-\frac{b^{2}}{D}-\frac{m_{k}^{2}}{2E_{k}^{2}}.$

Using Eq. (30) in the integral for the phase we can again write

$\Phi_{k}=-\int^{D}_{S}\frac{E_{0}B}{p_{0}(r)}E_{k}\frac{m_{k}^{2}}{2E_{k}^{2}}dr.$
(31)

Using the second equation in relativistic approximation given in Eqs. (20),
the phase becomes

$\displaystyle\Phi_{k}$
$\displaystyle=-\frac{m_{k}^{2}}{2E_{0}}\int^{D}_{S}\frac{E_{0}B}{p_{0}(r)}dr=$
(32)
$\displaystyle=\pm\frac{m_{k}^{2}}{2E_{0}}\int^{D}_{S}\sqrt{AB}\left(1-\frac{b^{2}A}{D}\right)^{-1/2}dr.$

For $b=0$ and $\gamma=1$, we recover the well-known result of the phase for a
radially propagating neutrino in Schwarzschild spacetime

$\Phi_{k}\simeq\frac{m_{k}^{2}}{2E_{0}}|r_{D}-r_{S}|\ .$ (33)

Now we would like to integrate the integral in Eq. (32) for two cases: (a) in
the first, we shall consider the case where a neutrino is produced in a
gravitational field and then propagates outward in the potential non-radially;
(b) in the second, we shall consider the case where a neutrino is produced at
a source $S$, travels outward, moves around the massive object (which is
effectively the gravitational lens), crossing the closest approach point at
$r=r_{0}$ and reaches the detector $D$. The second situation will be used for
calculating the oscillation probability in a setting of gravitational lensing
of neutrinos.

Let us first look at case (a). In the weak field approximation, we expand the
quantity under integral in Eq. (32) as

$\Phi_{k}\simeq\pm\frac{m_{k}^{2}}{2E_{0}}\int^{r_{D}}_{r_{S}}\left[\frac{1}{\sqrt{1-\frac{b^{2}}{r^{2}}}}+\frac{b^{2}(1-2\gamma)M/r}{(r^{2}-b^{2})\sqrt{1-\frac{b^{2}}{r^{2}}}}\right]dr.$
(34)

Now this can be easily integrated and the result is

$\displaystyle\Phi_{k}$
$\displaystyle=\frac{m_{k}^{2}}{2E_{0}}\Bigg{[}\sqrt{r_{D}^{2}-b^{2}}-\sqrt{r_{S}^{2}-b^{2}}+$
(35)
$\displaystyle+(2\gamma-1)M\left(\frac{r_{D}}{\sqrt{r_{D}^{2}-b^{2}}}-\frac{r_{S}}{\sqrt{r_{S}^{2}-b^{2}}}\right)\Bigg{]}.$

It can be seen from this expression that when $b=0$, we recover the phase for
the case of radial motion of neutrinos.

In case (b) a propagating neutrino moves in the vicinity of the lensing object
passing the closest point of approach at $r=r_{0}$. In this case, the integral
for the phase can be written as the sum of two parts taking into account the
sign of the momentum,

$\displaystyle\Phi_{k}(r_{S}\rightarrow r_{0}\rightarrow r_{D})$
$\displaystyle=\frac{m_{k}^{2}}{2E_{0}}\int^{r_{S}}_{r_{0}}\sqrt{\frac{AB}{\left(1-\frac{b^{2}A}{D}\right)}}dr+$
(36)
$\displaystyle+\frac{m_{k}^{2}}{2E_{0}}\int^{r_{D}}_{r_{0}}\sqrt{\frac{AB}{\left(1-\frac{b^{2}A}{D}\right)}}dr.$

The position of closest approach of the neutrinos can be obtained from the
equation

$\left(\frac{dr}{d\phi}\right)_{0}=\frac{p_{0}(r_{0})D}{J_{0}B}=0.$ (37)

In the weak-field approximation, we solve this equation to obtain the point of
closest approach as

$r_{0}\simeq b\left[1-(2\gamma-1)\frac{M}{b}\right].$ (38)

Now integrating Eq. (36) using the point of closest approach we get

$\displaystyle\Phi_{k}$
$\displaystyle=\frac{m_{k}^{2}}{2E_{0}}\Bigg{[}\sqrt{r_{D}^{2}-r_{0}^{2}}+\sqrt{r_{S}^{2}-r_{0}^{2}}+$
(39)
$\displaystyle+(2\gamma-1)M\left(\sqrt{\frac{r_{D}-r_{0}}{r_{D}^{2}+r_{0}^{2}}}+\sqrt{\frac{r_{S}-r_{0}}{r_{S}^{2}+r_{0}^{2}}}\right)\Bigg{]},$

or

$\displaystyle\Phi_{k}$
$\displaystyle\simeq\frac{m_{k}^{2}}{2E_{0}}\Bigg{[}\sqrt{r_{D}^{2}-b^{2}}+\sqrt{r_{S}^{2}-b^{2}}+$
(40)
$\displaystyle+(2\gamma-1)M\Big{(}\frac{b}{\sqrt{r_{S}^{2}-b^{2}}}+\frac{b}{\sqrt{r_{D}^{2}-b^{2}}}+$
$\displaystyle+\frac{r_{S}-b}{\sqrt{r_{S}+b}}++\frac{r_{S}-b}{\sqrt{r_{S}+b}}\Big{)}\Bigg{]}\
.$

Now, for $b<<r_{S,D}$ expanding the above quantity over $b/r_{S}$ and
$b/r_{D}$ and keeping terms of the order up to $(b^{2}/r_{S,D}^{2})$ we obtain

$\Phi_{k}=\frac{m_{k}^{2}}{2E_{0}}(r_{S}+r_{D})\left[1-\frac{b^{2}}{2r_{S}r_{D}}+\frac{2(2\gamma-1)M}{r_{S}+r_{D}}\right].$
(41)

As we can see, the phase now depends on $\gamma$ in the first order in
($M/r_{S,D}$). This is because, all the particles travelling non-radially will
be affected by the deformed geometry. This expression will be used to discuss
gravitational lensing of neutrinos in $\gamma$-spacetime in the following
section.

Figure 1: Schematic diagram for weak lensing of neutrinos in the
$\gamma$-spacetime. Neutrinos propagate from the source S to detector D in the
exterior of a static and non spherical massive object which is described by
the $\gamma$-metric.

Figure 2:  Top panel: Neutrino oscillation probability for $\gamma=1.0$ when
the mixing angle is $\alpha=\pi/6$. Bottom panel: Neutrino oscillation
probability for $\gamma=1.5$ when the mixing angle is $\alpha=\pi/6$. The
solid and the dashed curves represent normal hierarchy and inverted hierarchy
respectively. Values of the other parameters are as follows: $M_{\rm
ADM}=1M_{\odot}$, $\Delta m^{2}=10^{-3}\;{\rm eV}^{2}$, and the lightest
neutrino here is considered to be massless.

Figure 3: Top panel: Neutrino oscillation probability as a function of
azimuthal angle $\phi$ for $\gamma=1.0$ (solid line),$\gamma=0.5$ (dashed
line) and $\gamma=1.5$ (dotted line) for normal hierarchy $\Delta m^{2}>0$.
Bottom panel: Neutrino oscillation probability for $\gamma=1.0$ (solid
line),$\gamma=0.5$ (dashed line) and $\gamma=1.5$ (dotted line) for inverted
hierarchy $\Delta m^{2}<0$. The mixing angle here is $\alpha=\pi/6$. Values of
the other parameters are as follows: $M_{\rm ADM}=1M_{\odot}$, $\Delta
m^{2}=10^{-3}\;{\rm eV}^{2}$, and the lightest neutrino here is considered to
be massless.

##  V Neutrino oscillation probability

Let us now consider neutrinos with mass eigentaste $\nu_{i}$ travelling in the
$\gamma$-spacetime, through different classical paths, from a source S to meet
at a detector D. The neutrino flavor eigenstate $\nu_{\alpha}$, propagated
from the source to the detector through a path denoted by $p$, is given by,

$\ket{\nu_{\alpha}(t_{D},x_{D})}=N\sum_{i}U^{*}_{\alpha
i}\sum_{p}\exp\left(-i\Phi_{i}^{p}\ket{\nu_{i}(t_{S},x_{S})}\right).$ (42)

Here $\Phi_{i}^{p}$ is given by Eq. (41) with $b_{p}$ as the path dependent
impact parameter. If a neutrino is produced in $\alpha$ flavor eigenstate at
the source $S$ then the probability of it being detected in $\beta$ flavor at
the detector location, is given by Swami et al. (2020)

$\displaystyle\mathcal{P}_{\alpha\beta}^{\rm lens}$
$\displaystyle=|\braket{\nu_{\beta}}{\nu_{\alpha}(t_{D},x_{D})}|^{2}=$ (43)
$\displaystyle=|N|^{2}\sum_{i,j}U_{\beta i}U_{\beta j}^{*}U_{\alpha
j}U_{\alpha i}^{*}\sum_{p,q}\exp\left(-i\Delta\Phi_{ij}^{pq}\right),$

where

$|N|^{2}=\left(\sum_{i}|U_{\alpha
i}|^{2}\sum_{p,q}\left(-i\Delta\Phi_{ii}^{pq}\right)\right)$ (44)

is the normalization constant and $\Delta\Phi_{ij}^{pq}$ can be given as

$\displaystyle\Delta\Phi_{ij}^{pq}$
$\displaystyle=\Phi_{p}^{i}-\Phi_{q}^{j}=\Delta m_{ij}^{2}{\rm A}_{pq}+\Delta
b_{pq}^{2}{\rm B}_{ij},$ (45)

where

$\displaystyle{\rm A}_{pq}$
$\displaystyle=\frac{r_{S}+r_{D}}{2E_{0}}\Bigg{(}1+\frac{2(2\gamma-1)M}{2r_{S}r_{D}}-\frac{\Sigma
b_{pq}^{2}}{4r_{S}r_{D}}\Bigg{)},$ (46) $\displaystyle{\rm B}_{ij}$
$\displaystyle=-\frac{\Sigma
m_{ij}^{2}}{8E_{0}}\left(\frac{1}{r_{S}}+\frac{1}{r_{D}}\right).$

In the above equations, the quantities $\Delta m_{ij}^{2}$, $\sum m_{ij}^{2}$,
$\Delta b_{pq}^{2}$ and $\sum b_{pq}^{2}$ are given by

$\displaystyle\Delta m_{ij}^{2}=m_{i}^{2}-m_{j}^{2},\ \ \ \sum
m_{ij}^{2}=m_{i}^{2}+m_{j}^{2},$ (47) $\displaystyle\Delta
b_{pq}^{2}=b_{p}^{2}-b_{q}^{2},\ \ \ \sum b_{pq}^{2}=b_{p}^{2}+b_{q}^{2},$

respectively. Based upon these, we can make the following observations.

  * •

${\rm A}_{pq}$ and ${\rm B}_{ij}$ are symmetric under interchange of their
indices, ${\rm A}_{pq}={\rm A}_{qp}$ and ${\rm B}_{ij}={\rm B}_{ji}$.

  * •

The oscillation probability depends on the sum of individual mass squared of
neutrinos $\sum m_{ij}^{2}$ through $\Delta b_{pq}^{2}$.

  * •

For those paths for which $\Delta b_{pq}^{2}$ vanishes,
$\mathcal{P}_{\alpha\beta}^{\rm lens}$ is invariant under the symmetry
$m_{i}^{2}\rightarrow m_{i}^{2}+C$.

  * •

For the paths for which $\Delta b_{pq}^{2}$ does not vanish, the symmetry
$m_{i}^{2}\rightarrow m_{i}^{2}+C$ is broken. The shift implies ${\rm
B}_{ij}\rightarrow{\rm B}_{ij}+2C$. Therefore, a generic configuration of
source S, lens and the detector D will retain the information about $\sum
m_{ij}^{2}$. Hence the observation of lensed neutrinos would carry information
about their absolute masses Swami et al. (2020).

##  VI Gravitational lensing of neutrinos in $\gamma$-spacetime

We can see clearly from the expression for phase difference derived in
previous section that the probability of neutrino oscillation depends on
$\Sigma m_{ij}^{2}$ through a path dependent term $\Delta b_{pq}^{2}$. This
was also shown in two and three flavour neutrino oscillation case in Swami et
al. (2020) for Schwarzschild. Similarly, in the $\gamma$-spacetime, it appears
that the absolute neutrino mass dependent effect is same as in Schwarzschild,
since ${\rm B}_{ij}$ does not depend on $\gamma$. However, we shall see later
that the impact parameters $b_{p}$ depend on $\gamma$ and hence we should be
able to see the effect of the quadruple deformation parameter $\gamma$ in both
terms of the phase difference (45).

Now we can simplify the expression for the oscillation probability in order to
understand the effect of $\gamma$ more clearly. For simplicity, first we shall
consider a toy model with two neutrino flavors. Specifically, we shall
evaluate the two flavor neutrino oscillation probability at a generic point in
a plane connecting the source, the lens and the detector in the weak field
limit. Therefore we substitute equation (45) in equation (43) to get

$\mathcal{P}^{\rm lens}_{\alpha\beta}=|N|^{2}\left[\sum_{i,j}U_{\beta
i}U^{*}_{\beta j}U_{\alpha j}U^{*}_{\alpha
i}\Bigg{(}\sum_{p=q}\exp\left(-i\Delta m^{2}_{ij}{\rm
A}_{pp}\right)+2\sum_{p>q}\cos\left(\Delta b_{pq}^{2}\,{\rm
B}_{ij}\right)\exp\left(-i\Delta m^{2}_{ij}\,{\rm
A}_{pq}\right)\Bigg{)}\right]\ ,$ (48)

where $|N|^{2}$ is given by

$|N|^{2}=\left(N_{\rm path}+\sum_{i}|U_{\alpha i}|^{2}\sum_{q>p}2\cos(\Delta
b_{pq}^{2}{\rm B}_{ii})\right)^{-1}.$ (49)

Here, for simplicity, we consider neutrino propagation in the equatorial
($\theta=\pi/2$) plane and in this case, $N_{\rm path}=2$. Now the general
expression of probability reduces to a simpler form

$\displaystyle\mathcal{P}_{\alpha\beta}^{\rm lens}$
$\displaystyle=|N|^{2}\Bigg{[}2\sum_{i}|U_{\beta i}|^{2}|U_{\alpha
i}|^{2}(1+\cos(\Delta b^{2}{\rm B}_{ii}))+$ (50) $\displaystyle+\sum_{i,j\neq
i}U_{\beta i}U_{\beta j}^{*}U_{\alpha j}U_{\alpha i}^{*}\big{[}\exp(-i\Delta
m_{ij}^{2}{\rm A}_{11})+$ $\displaystyle+\exp(-i\Delta m_{ij}^{2}{\rm
A}_{22})+$ $\displaystyle+2\cos(\Delta b^{2}{\rm B}_{ij})\exp(-i\Delta
m_{ij}^{2}{\rm A}_{12})\big{]}\Bigg{]}\ ,$

with ${\rm A}_{pq}$ and ${\rm B}_{ij}$ being given by equation (46) and the
normalization factor reduces to

$|N|^{2}=\left(2+2\sum_{i}|U_{\alpha i}|^{2}\cos(\Delta b^{2}{\rm
B}_{ii})\right)^{-1}\ ,$ (51)

where $\Delta b^{2}=\Delta b_{12}^{2}$. Now for the qualitative and
quantitative understanding of neutrino oscillation in the $\gamma$-spacetime,
we consider the simplest two neutrino flavors oscillation,
$\nu_{e}\rightarrow\nu_{\mu}$. In this case, the oscillation probability
becomes

$\displaystyle\mathcal{P}_{e\mu}^{\rm lens}$
$\displaystyle=|N|^{2}\sin^{2}2\alpha\Bigg{[}\sin^{2}\left(\Delta
m^{2}\frac{{\rm A}_{11}}{2}\right)+\sin^{2}\left(\Delta m^{2}\frac{{\rm
A}_{22}}{2}\right)+$ (52) $\displaystyle-\cos(\Delta b^{2}{\rm
B}_{12})\cos(\Delta m^{2}{\rm A}_{12})+\frac{1}{2}\cos(\Delta b^{2}{\rm
B}_{11})+$ $\displaystyle+\frac{1}{2}\cos(\Delta b^{2}{\rm B}_{22})\Bigg{]},$

where

$|N|^{2}=\left[2(1+\cos^{2}\alpha\cos(\Delta b^{2}{\rm
B}_{11})\sin^{2}\alpha\cos(\Delta b^{2}{\rm B}_{22}))\right]^{-1},$ (53)

and $\Delta m^{2}=\Delta m_{21}^{2}$.

###  VI.1 Two flavor case: numerical results

Figure 4: Left panel: Neutrino oscillation probability as a function of
azimuthal angle $\phi$ for $\gamma=1.0$, $0.5$ and $1.5$ with $\Delta m^{2}>0$
(normal hierarchy). Here the solid line represents $m_{1}=0$ eV and
$\gamma=1.0$, the dashed line represents $m_{1}=0.03$ eV and $\gamma=1.0$, the
dotted line represents $m_{1}=0$ eV and $\gamma=1.5$, the dotted-dashed line
represents $m_{1}=0.03$ eV and $\gamma=1.5$, the thick solid line represents
$m_{1}=0$ eV and $\gamma=0.5$, and the triple-dotted dashed line represents
$m_{1}=0.03$ eV and $\gamma=0.5$. Right panel: Neutrino oscillation
probability as a function of azimuthal angle $\phi$ for $\gamma=1.0$, $0.5$
and $1.5$ with $\Delta m^{2}<0$ (inverted hierarchy). Here the solid line
represents $m_{1}=0$ eV and $\gamma=1.0$, the dashed line represents
$m_{1}=0.03$ eV and $\gamma=1.0$, the dotted line represents $m_{1}=0$ eV and
$\gamma=1.5$, the dotted-dashed line represents $m_{1}=0.03$ eV and
$\gamma=1.5$, the thick solid line represents $m_{1}=0$ eV and $\gamma=0.5$,
and the triple-dotted dashed line represents $m_{1}=0.03$ eV and $\gamma=0.5$.
The other parameters are as follows: $M_{\rm ADM}=1M_{\odot}$, $\Delta
m^{2}=10^{-3}{\rm eV}^{2}$.

Figure 5: Zoomed-in version of the $\gamma$ dependence for the neutrino
oscillation probability. Left panel: Neutrino oscillation probability as a
function of azimuthal angle $\phi$ for $\gamma=0.5$ (thin solid line),
$\gamma=0.75$ (dashed line), $\gamma=1.0$ (dotted line), $\gamma=1.25$
(dotted-dashed line), $\gamma=1.5$ (thick solid line) in case of normal
hierarchy. Right panel: Neutrino oscillation probability for $\gamma=0.5$
(thin solid line), $\gamma=0.75$ (dashed line), $\gamma=1.0$ (dotted line),
$\gamma=1.25$ (dotted-dashed line), $\gamma=1.5$ (thick solid line) in case of
inverted hierarchy.

Now for a quantitative understanding, we would like to see how the probability
changes over some change in the lensing parameters. Therefore it is useful
express the impact parameter $b_{p}$ in terms of the geometric quantities of
the system. To this aim we can refer to Figure 1, which shows a schematic
representation of a weak lensing event in $\gamma$-spacetime. Here we have the
source (S) of neutrinos, the gravitational lens characterized by the
$\gamma$-metric and the detector (D) on a plane. Physical distances from the
source to the lens and from the lens to the detector are $r_{S}$ and $r_{D}$,
respectively, in the $(x,y)$ coordinate system. We can also consider another
coordinate system $(x^{\prime},y^{\prime})$ obtained by rotating the original
coordinate system $(x,y)$ by an angle $\varphi$ such that
$x^{\prime}=x\cos\varphi+y\sin\varphi$ and
$y^{\prime}=-x\sin\varphi+y\cos\varphi$. In the rotated frame, the angle of
deflection of the neutrino, $\delta$ from its original path is related to the
impact parameter $b$ by the following relation

$\delta\sim\frac{y_{D}^{\prime}-b}{x_{D}^{\prime}}=-\frac{4\gamma
M}{b}=-\frac{2R_{x}}{b},$ (54)

where $(x_{D}^{\prime},y_{D}^{\prime})$ is the location of the detector,
$R_{x}=2\gamma M$ and we have used the expression for $\delta$ in the weak
lensing limit. Now using the identity $\sin\varphi=b/r_{S}$ from Figure 1, the
previous equation Eq. (54) can be written as

$(2R_{x}x_{D}+by_{D})\sqrt{1-\frac{b^{2}}{r_{S}^{2}}}=b^{2}\left(\frac{x_{D}}{r_{S}}+1\right)-\frac{2R_{x}by_{D}}{r_{S}}.$
(55)

Solution of this equation gives the impact parameters in terms of $r_{S}$,
$R_{x}$ and the lensing location $(x_{D},y_{D})$. As a numerical exercise we
can consider the Sun-Earth system with typical values of the geometrical
quantities and assume that the gravitational field of the Sun is represented
by the $\gamma$-metric while the Earth is taken to be the detector. The source
is considered to be situated behind the Sun and it emits relativistic
neutrinos with typical energy $E_{0}=10\ \text{MeV}$. Now assuming a circular
trajectory of the detector around the Sun such that $x_{D}=r_{D}\cos\phi$ and
$y_{D}=r_{D}\sin\phi$, we can numerically solve the quartic polynomial given
by Eq. (55) in the equatorial plane and obtain two positive real roots $b_{1}$
and $b_{2}$ for every $\phi$. In this system, $r_{D}=10^{8}\ \text{km}$ and
$r_{S}=10^{5}r_{D}$. In this numerical exercise, neutrino oscillation
probability is calculated only for those value of $b_{p}$ for which $R_{x}\ll
b_{p}\ll r_{D}$. Other relevant parameters are $M_{\rm ADM}=\gamma
M=1M_{\odot}$ (as defined in (10)) and $|\Delta m^{2}|=10^{-3}\;{\rm eV}^{2}$.
Note that, these numbers are for illustrative purposes only and in a realistic
scenario proper numerical values of the geometric parameters of the model have
to be selected.

Oscillation probabilities for the two-flavor toy model of neutrinos are shown
in Figure 2, 3, 4, and 5 as a function of the azimuthal angle $\phi$. Note
that our main aim is to investigate the dependence of the oscillation
probability on the deformation parameter $\gamma$. Figure 2 shows the neutrino
oscillation probability $\nu_{e}\rightarrow\nu_{\mu}$ for $\gamma=1.0$ (top
panel) and $\gamma=1.5$ (bottom panel)111It is worth mentioning that the
actual realistic values of the $\gamma$ parameter are expected to be much
closer to $1$ both for the solar system as well as for any astrophysical
compact object. Here, we choose large values of $\gamma$ for illustrative
purposes only.. In each of these two plots shown in Figure 2, the solid line
corresponds to normal hierarchy ($\Delta m^{2}>0$) and the dashed line
corresponds to inverted hierarchy ($\Delta m^{2}<0$). The mixing angle here is
$\alpha=\pi/6$. The first thing we notice from these three plots is that the
oscillation probability is different for normal and inverted mass ordering
(except for some values of $\phi$) confirming the absolute mass dependent
effects in gravitational lensing Swami et al. (2020). This is in sharp
contrast with two flavour vacuum oscillation probability in flat spacetime
where the sign of $\Delta m^{2}$ does not make any distinction. We also notice
an increase in frequency for higher values of $\phi$. Figure 3 shows the same
results in order to better highlight the dependence of the probability on
$\gamma$.

The top panel in Figure 3 shows the conversion probability for $\gamma=1$
(solid line), $\gamma=0.5$ (dashed line) and $\gamma=1.5$ (dotted line) for
normal hierarchy ($\Delta m^{2}>0$). The bottom panel shows the conversion
probability for $\gamma=1$ (solid line), $\gamma=0.5$ (dashed line) and
$\gamma=1.5$ (dotted line) for inverted hierarchy ($\Delta m^{2}<0$). It can
be clearly seen from the plots that the conversion probability is different
for different values of the $\gamma$ parameter even though the observable mass
$M_{\rm ADM}=\gamma M$ is same. In standard optical observations, such as
deflection of light in weak-field limit, it is not possible to constrain the
parameter $\gamma$ independently as it appears together with the mass
parameter $M$, creating a degeneracy. For example, the deflection angle in the
weak-field limit is $\delta=-4\gamma M/b=-4M_{\rm ADM}/b$ which remains
unchanged if we change $\gamma$ and $M$ is such a way that $M_{\rm ADM}$
remains the same. However, this degeneracy is broken in the phase difference
of neutrino oscillation and appears in ${\rm A}_{pq}$ as
$2(2\gamma-1)M/2r_{S}r_{D}$ in the coefficients of neutrino mass squared
difference ($\Delta m_{ij}^{2}$). Note that, ${\rm B}_{ij}$ is independent of
$\gamma$ however, $b_{p,q}$ and hence $\Delta b_{pq}^{2}$ and $\Sigma
b_{pq}^{2}$ are $\gamma$ dependent as the paths will be affected by the
deviations from the spherical symmetry (see Eq. (55)). Therefore, the
deviations of the geometry from spherical symmetry can in principle be
constrained from neutrino observations. For values of $\gamma<1$, we see a
larger change in conversion probability than for values of $\gamma>1$. With
decrease in $\gamma$, the conversion probability increases for both normal and
inverted hierarchy as seen from Figure 5. Also it is important to note that in
Figures 2, 3 and 5, the lightest neutrino is assumed to be massless.

In Figure 4, we compare the conversion probability of the massless case for
lightest neutrino with the massive case. The left panel in Figure 4 shows
conversion probabilities for the case of normal hierarchy for
$\gamma=0.5,1\text{ and},1.5$ when the lightest neutrino mass $m_{1}=0\ {\rm
eV}$ and $0.03\ {\rm eV}$. Similarly the left panel of Figure 4 shows the case
of inverted hierarchy for $\gamma=0.5,1\text{ and }1.5$ when the lightest
neutrino mass is $m_{1}=0\ {\rm eV}$ and $0.03\ {\rm eV}$. We notice similar
effects of $\gamma$ on the conversion probability in this case as well. This,
in principle, confirms that neutrino probability observation can constrain
deviations from spherical symmetry of a gravitational system.

Figure 6: Probability of $\nu_{\alpha}\rightarrow\nu_{\beta}$ conversion, for
different $\alpha$ and $\beta$, as function of the azimuthal angle $\phi$ in
the three flavor case for $\gamma=1$ (spherical symmetry). The solid, dashed
and the dotted lines represent the $\nu_{e}\rightarrow\nu_{\mu}$,
$\nu_{e}\rightarrow\nu_{\tau}$ and $\nu_{\mu}\rightarrow\nu_{\tau}$ conversion
probabilities, respectively. We take $r_{D}=10^{8}$ km, $r_{S}=10^{5}r_{D}$,
$M_{\rm ADM}=1M_{\odot}$ and $E_{0}=10$ MeV. Neutrino mass squared
differences, mixing angles and the Dirac CP phase are taken from the latest
global fit Esteban et al. (2020).

Figure 7: Oscillation probability for the three flavor case in normal
hierarchy with $m_{1}=0$. Top panel: Probability of
$\nu_{e}\rightarrow\nu_{\mu}$ conversion for $\gamma=1.0$ (solid line),
$\gamma=0.5$ (dashed line), $\gamma=1.5$ (dotted line). Mid panel: Probability
of $\nu_{e}\rightarrow\nu_{\tau}$ conversion for $\gamma=1.0$ (solid line),
$\gamma=0.5$ (dashed line), $\gamma=1.5$ (dotted line). Bottom panel:
Probability of $\nu_{\mu}\rightarrow\nu_{\tau}$ conversion for $\gamma=1.0$
(solid line), $\gamma=0.5$ (dashed line), $\gamma=1.5$ (dotted line). We take
$r_{D}=10^{8}$ km, $r_{S}=10^{5}r_{D}$, $M_{\rm ADM}=1M_{\odot}$ and
$E_{0}=10$ MeV. Neutrino mass squared differences, mixing angles and the Dirac
CP phase are taken from the latest global fit Esteban et al. (2020).

Figure 8: Oscillation probability for the three flavor case in inverted
hierarchy with $m_{3}=0$. Top panel: Probability of
$\nu_{e}\rightarrow\nu_{\mu}$ conversion for $\gamma=1.0$ (solid line),
$\gamma=0.5$ (dashed line), $\gamma=1.5$ (dotted line). Mid panel: Probability
of $\nu_{e}\rightarrow\nu_{\tau}$ conversion for $\gamma=1.0$ (solid line),
$\gamma=0.5$ (dashed line), $\gamma=1.5$ (dotted line). Bottom panel:
Probability of $\nu_{\mu}\rightarrow\nu_{\tau}$ conversion for $\gamma=1.0$
(solid line), $\gamma=0.5$ (dashed line), $\gamma=1.5$ (dotted line). We take
$r_{D}=10^{8}$ km, $r_{S}=10^{5}r_{D}$, $M_{\rm ADM}=1M_{\odot}$ and
$E_{0}=10$ MeV. Neutrino mass squared differences, mixing angles and the Dirac
CP phase are taken from the latest global fit Esteban et al. (2020).

Figure 9: Oscillation probability for the three flavor case in normal
hierarchy with $m_{1}=0.08$ eV. Top panel: Probability of
$\nu_{e}\rightarrow\nu_{\mu}$ conversion for $\gamma=1.0$ (solid line),
$\gamma=0.5$ (dashed line), $\gamma=1.5$ (dotted line). Mid panel: Probability
of $\nu_{e}\rightarrow\nu_{\tau}$ conversion for $\gamma=1.0$ (solid line),
$\gamma=0.5$ (dashed line), $\gamma=1.5$ (dotted line). Bottom panel:
Probability of $\nu_{\mu}\rightarrow\nu_{\tau}$ conversion for $\gamma=1.0$
(solid line), $\gamma=0.5$ (dashed line), $\gamma=1.5$ (dotted line). We take
$r_{D}=10^{8}$ km, $r_{S}=10^{5}r_{D}$, $M_{\rm ADM}=1M_{\odot}$ and
$E_{0}=10$ MeV. Neutrino mass squared differences, mixing angles and the Dirac
CP phase are taken from the latest (NuFIT 5.0 (2020)) global fit Esteban et
al. (2020).

###  VI.2 Three flavor case: numerical results

In this section, we shall apply the formalism developed earlier to the
realistic three flavor neutrino model. We start with Eq. (50) and write down
the oscillation probability expressions for $\nu_{e}\rightarrow\nu_{\mu}$,
$\nu_{e}\rightarrow\nu_{\tau}$ and $\nu_{\mu}\rightarrow\nu_{\tau}$
conversion. The $3\times 3$ matrix $U$ is the usual PMNS matrix parametrized
by three angles $\theta_{12}$, $\theta_{13}$, $\theta_{23}$ and the Dirac CP
phase $\delta_{\rm CP}$. Similar to the two flavor case, we solve the quartic
polynomial in Eq. (55) numerically to find the lensing locations and calculate
the different oscillation probabilities in these locations. The numerical
values of the parameters here are taken from the latest global best fit of
neutrino oscillation data Esteban et al. (2020). The parameters are
$\theta_{12}=33.44^{\circ}$ ($\theta_{12}=33.45^{\circ}$),
$\theta_{13}=8.57^{\circ}$ ($\theta_{13}=8.61^{\circ}$),
$\theta_{23}=49.0^{\circ}$ ($\theta_{23}=49.3^{\circ}$), $\Delta
m_{21}^{2}=7.42\times 10^{-5}\;{\rm eV}^{2}$, $\Delta m_{31}^{2}=2.514\times
10^{-3}\;{\rm eV}^{2}$ ($\Delta m_{32}^{2}=-2.497\times 10^{-3}\;{\rm
eV}^{2}$) and $\delta_{\rm CP}=222^{\circ}$ ($\delta_{\rm CP}=286^{\circ}$)
for normal (inverted) hierarchy.

The results for conversion probability in three flavor neutrino oscillation
are shown in Figures 6-9. $\mathcal{P}_{e\mu}$, $\mathcal{P}_{e\tau}$ and
$\mathcal{P}_{\mu\tau}$ for $\gamma=1$ and $m_{1}=0\;{\rm eV}$ is plotted in
Figure 6. Figure 7 shows the $\gamma$ dependence of the conversion
probabilities when the lightest neutrino is massless $m_{1}=0\;{\rm eV}$ for
normal hierarchy. The top, mid and the bottom panel shows
$\mathcal{P}_{e\mu}$, $\mathcal{P}_{e\tau}$ and $\mathcal{P}_{\mu\tau}$,
respectively, for $\gamma=1,0.5$ and $1.5$. Figure 8 shows the $\gamma$
dependence of the conversion probabilities when the lightest neutrino is
massless $m_{3}=0\;{\rm eV}$ for inverted hierarchy. The top, mid and the
bottom panels show $\mathcal{P}_{e\mu}$, $\mathcal{P}_{e\tau}$ and
$\mathcal{P}_{\mu\tau}$, respectively, for $\gamma=1,0.5$ and $1.5$.
Significant differences can be noticed between the normal and inverted
hierarchies from comparing Figures 7 and 8. Figure 9 shows the $\gamma$
dependence of the conversion probabilities when the lightest neutrino is
massive $m_{1}=0.08\;{\rm eV}$. Similarly to the two flavor case, we notice
that the presence of the parameter $\gamma$ induces significant effects in
three flavor case as well. This is especially noticeable when the lightest
neutrino is massive, when we observe the most significant difference in the
conversion probabilities for different values of $\gamma$.

##  VII Conclusion

Neutrino flavor oscillation is an interesting phenomenon which has the
potential to resolve several longstanding questions in fundamental physics. In
flat spacetime, with or without matter effects, the flavor oscillation depends
only on the difference of the squared masses and does not give us any
information about individual neutrino masses. However, neutrino flavor
oscillation in a curved spacetime depends on the absolute neutrino masses, as
shown in several earlier works, and may be relevant for neutrinos produced by
highly energetic astrophysical phenomena. In particular, the effects of
gravitational lensing on neutrino oscillations appear to be promising in
determining the neutrino masses Swami et al. (2020).

This dependency on the individual neutrino masses exists for a wide class of
trajectories, including of course geodesics, thus providing a novel potential
experimental method to determine masses through the measurement of flavor
transition probability. In this work, we extended this method to a larger
class of geometries which include mass quadrupole moment by considering a
static, axially symmetric and asymptotically flat solution of Einstein’s field
equations in vacuum known as $\gamma$-metric. The additional parameter
$\gamma$ in the line element describes deviations of the source from spherical
symmetry. We showed how deviations from the usual Schwarzschild metric for the
background geometry in which neutrinos travel affects the flavor transition
probability and thus potential observations. Conversely, if one knows the
transition probabilities, this method provides in principle a new and unique
tool to measure $\gamma$ and hence the quadrupole moment of a spacetime which
would otherwise be impossible by usual weak optical lensing experiments.

The numerical results of oscillation probability in the $\gamma$-spacetime for
two and three flavor neutrinos show that the probability increases for
decreasing $\gamma$ at certain locations of the detector on the orbit. This is
true for both normal and inverted hierarchies when the lightest neutrino is
massless, and changes in the value of $\gamma$ alter the peak of the
oscillation probability. Similar qualitative but different quantitative
behaviour is observed when the lightest neutrino is massive showing that the
effects due to the absolute neutrino mass depend also on the value of
$\gamma$. This suggest the possibility of measuring the neutrino masses if the
quadruple moment of the source is known or, vice versa, measuring the
quadrupole moment if the oscillation probabilities are known.

While the toy model presented here can be considered as a first step, a more
realistic approach that may be pursued in the future will require the use of a
quantum field theoretical treatment Grimus (2020); Capolupo et al. (2020) and
the choice of a geometry that describes the exterior of rotating compact
objects. Additionally, to have a complete picture of gravitational effects on
neutrino oscillations, one would need to consider effects of gravitational
decoherence Swami et al. (2021), gravitational spin-flip transitions of
neutrinos Sorge and Zilio (2007); Dvornikov (2019) etc. We would also like to
point out some earlier work Blasone1 (1995); Blasone2 (1998) where it was
shown using a QFT-based formula that neutrino oscillations (in vacuum-flat
spacetime) in the “non-relativistic" regime depends on individual neutrino
masses.

In the future the growth of observations of neutrinos from distant highly
energetic astrophysical sources will play an important role in our
understanding of compact objects in the universe as well as fundamental
physics. In particular, gravitational lensing of neutrinos can provide a
unique probe to test the nature of the source and the validity of modified
gravity theories. In this context, detecting neutrinos from extra-solar
sources lensed by compact objects is expected to provide constraints on the
geometry ond thus the nature of the lensing object. In future work we plan to
extend the framework to modified theories of gravity to determine how future
neutrino observations may be used to test alternative theories. The hope is
that the growth of multi-messenger astronomy Mészáros et al. (2019) and extra-
galactic neutrino astronomy Aartsen et al. (2018), will make it possible to
test the nature of the geometry surrounding compact objects, put constraints
on the validity of modified theories of gravity and explore fundamental
properties of neutrino physics such as their absolute masses, which remain
unknown with present oscillation experiments.

###### Acknowledgements.

The authors thank Dr. Yong Tang for careful reading of the manuscript and for
useful comments. HC acknowledges support from University of Chinese Academy of
Sciences through Special Research Assistant fellowship (SRA). DB acknowledges
the support from Early Career Research Award from DST-SERB, Government of
India (reference number: ECR/2017/001873). AA and BA acknowledge Grants of the
Uzbekistan Ministry for Innovative Development and Grant NT-01 of the Abdus
Salam International Centre of Theoretical Physics. AA also acknowledges the
support from Chinese academy of sciences through PIFI foundation.

## References

  * Zyla et al. (2020) P. A. Zyla et al. (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Esteban et al. (2020) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, JHEP 09, 178 (2020), eprint 2007.14792.
  * de Salas et al. (2021) P. F. de Salas, D. V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. A. Ternes, M. Tórtola, and J. W. F. Valle, JHEP 02, 071 (2021), eprint 2006.11237.
  * Aker et al. (2019) M. Aker et al. (KATRIN), Phys. Rev. Lett. 123, 221802 (2019), eprint 1909.06048.
  * Aghanim et al. (2018) N. Aghanim et al. (Planck) (2018), eprint 1807.06209.
  * Abe et al. (2020) K. Abe et al. (T2K), Nature 580, 339 (2020), [Erratum: Nature 583, E16 (2020)], eprint 1910.03887.
  * Wudka (1991) J. Wudka, Mod. Phys. Lett. A 6, 3291 (1991).
  * Grossman and Lipkin (1997) Y. Grossman and H. J. Lipkin, Phys. Rev. D 55, 2760 (1997), eprint hep-ph/9607201.
  * Cardall and Fuller (1997) C. Y. Cardall and G. M. Fuller, Phys. Rev. D 55, 7960 (1997), eprint hep-ph/9610494.
  * Piriz et al. (1996) D. Piriz, M. Roy, and J. Wudka, Phys. Rev. D 54, 1587 (1996), eprint hep-ph/9604403.
  * Fornengo et al. (1997) N. Fornengo, C. Giunti, C. W. Kim, and J. Song, Phys. Rev. D 56, 1895 (1997), eprint hep-ph/9611231.
  * Ahluwalia and Burgard (1996) D. V. Ahluwalia and C. Burgard, Gen. Rel. Grav. 28, 1161 (1996), eprint gr-qc/9603008.
  * Bhattacharya et al. (1999) T. Bhattacharya, S. Habib, and E. Mottola, Phys. Rev. D 59, 067301 (1999).
  * Pereira and Zhang (2000) J. G. Pereira and C. M. Zhang, Gen. Rel. Grav. 32, 1633 (2000), eprint gr-qc/0002066.
  * Crocker et al. (2004) R. M. Crocker, C. Giunti, and D. J. Mortlock, Phys. Rev. D 69, 063008 (2004), eprint hep-ph/0308168.
  * Lambiase et al. (2005) G. Lambiase, G. Papini, R. Punzi, and G. Scarpetta, Phys. Rev. D 71, 073011 (2005), eprint gr-qc/0503027.
  * Godunov and Pastukhov (2011) S. I. Godunov and G. S. Pastukhov, Phys. Atom. Nucl. 74, 302 (2011), eprint 0906.5556.
  * Ren and Zhang (2010) J. Ren and C.-M. Zhang, Class. Quant. Grav. 27, 065011 (2010), eprint 1002.0648.
  * Geralico and Luongo (2012) A. Geralico and O. Luongo, Phys. Lett. A 376, 1239 (2012), eprint 1202.5408.
  * Chakraborty (2014) S. Chakraborty, Class. Quant. Grav. 31, 055005 (2014), eprint 1309.0693.
  * Visinelli (2015) L. Visinelli, Gen. Rel. Grav. 47, 62 (2015), eprint 1410.1523.
  * Chakraborty (2015) S. Chakraborty, JCAP 10, 019 (2015), eprint 1506.02647.
  * Zhang and Li (2016) Y.-H. Zhang and X.-Q. Li, Nucl. Phys. B 911, 563 (2016), eprint 1606.05960.
  * Alexandre and Clough (2018) J. Alexandre and K. Clough, Phys. Rev. D 98, 043004 (2018), eprint 1805.01874.
  * Blasone et al. (2020) M. Blasone, G. Lambiase, G. G. Luciano, L. Petruzziello, and L. Smaldone, Class. Quant. Grav. 37, 155004 (2020), eprint 1904.05261.
  * Buoninfante et al. (2020) L. Buoninfante, G. G. Luciano, L. Petruzziello, and L. Smaldone, Phys. Rev. D 101, 024016 (2020), eprint 1906.03131.
  * Boshkayev et al. (2020) K. Boshkayev, O. Luongo, and M. Muccino, Eur. Phys. J. C 80, 964 (2020), eprint 2010.08254.
  * Mandal (2021) S. Mandal, Nucl. Phys. B 965, 115338 (2021).
  * Koutsoumbas and Metaxas (2020) G. Koutsoumbas and D. Metaxas, Gen. Rel. Grav. 52, 102 (2020), eprint 1909.02735.
  * Swami et al. (2020) H. Swami, K. Lochan, and K. M. Patel, Phys. Rev. D 102, 024043 (2020), eprint 2002.00977.
  * Capolupo et al. (2020) A. Capolupo, G. Lambiase, and A. Quaranta, Phys. Rev. D 101, 095022 (2020), eprint 2003.00516.
  * Zipoy (1966) D. M. Zipoy, Journal of Mathematical Physics 7, 1137 (1966), eprint https://doi.org/10.1063/1.1705005, URL https://doi.org/10.1063/1.1705005.
  * Voorhees (1970) B. H. Voorhees, Phys. Rev. D 2, 2119 (1970).
  * Virbhadra (1996) K. S. Virbhadra (1996), eprint gr-qc/9606004.
  * Kodama and Hikida (2003) H. Kodama and W. Hikida, Class. Quant. Grav. 20, 5121 (2003), eprint gr-qc/0304064.
  * Chakrabarty et al. (2018) H. Chakrabarty, C. A. Benavides-Gallego, C. Bambi, and L. Modesto, JHEP 03, 013 (2018), eprint 1711.07198.
  * Lukes-Gerakopoulos (2012) G. Lukes-Gerakopoulos, Phys. Rev. D 86, 044013 (2012), eprint 1206.0660.
  * Herrera et al. (1999) L. Herrera, F. M. Paiva, and N. O. Santos, J. Math. Phys. 40, 4064 (1999), eprint gr-qc/9810079.
  * Chowdhury et al. (2012) A. N. Chowdhury, M. Patil, D. Malafarina, and P. S. Joshi, Phys. Rev. D 85, 104031 (2012), eprint 1112.2522.
  * Boshkayev et al. (2016) K. Boshkayev, E. Gasperin, A. C. Gutierrez-Pineres, H. Quevedo, and S. Toktarbay, Phys. Rev. D 93, 024024 (2016), eprint 1509.03827.
  * Toshmatov and Malafarina (2019) B. Toshmatov and D. Malafarina, Phys. Rev. D 100, 104052 (2019), eprint 1910.11565.
  * Abdikamalov et al. (2019) A. B. Abdikamalov, A. A. Abdujabbarov, D. Ayzenberg, D. Malafarina, C. Bambi, and B. Ahmedov, Phys. Rev. D 100, 024014 (2019), eprint 1904.06207.
  * Pontecorvo (1957) B. Pontecorvo, Zh. Eksp. Teor. Fiz. 34, 247 (1957).
  * Maki et al. (1962) Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys. 28, 870 (1962).
  * Pontecorvo (1967) B. Pontecorvo, Zh. Eksp. Teor. Fiz. 53, 1717 (1967).
  * Bilenky and Pontecorvo (1978) S. M. Bilenky and B. Pontecorvo, Phys. Rept. 41, 225 (1978).
  * Stewart et al. (1982) B. W. Stewart, D. Papadopoulos, L. Witten, R. Berezdivin, and L. Herrera, Gen. Relat. Gravit. 14, 97 (1982).
  * Herrera et al. (2005) L. Herrera, G. Magli, and D. Malafarina, Gen. Rel. Grav. 37, 1371 (2005), eprint gr-qc/0407037.
  * Hernandez-Pastora et al. (2016) J. L. Hernandez-Pastora, L. Herrera, and J. Martin, Class. Quant. Grav. 33, 235005 (2016), eprint 1607.02315.
  * Hernández-Pastora and Martín (1994) J. L. Hernández-Pastora and J. Martín, Gen. Rel. Grav. 26, 877 (1994).
  * Grimus (2020) W. Grimus, J. Phys. G 47, 085004 (2020), eprint 1910.13776.
  * Swami et al. (2021) H. Swami, K. Lochan, and K. M. Patel (2021), eprint 2106.07671.
  * Sorge and Zilio (2007) F. Sorge and S. Zilio, Class. Quant. Grav. 24, 2653 (2007).
  * Dvornikov (2019) M. Dvornikov, Phys. Rev. D 99, 116021 (2019), eprint 1902.11285.
  * Blasone1 (1995) M. Blasone and G. Vitiello, Annals Phys. 244, 283-311 (1995), [erratum: Annals Phys. 249, 363-364 (1996)] eprint hep-ph/9501263.
  * Blasone2 (1998) M. Blasone, P. A. Henning and G. Vitiello, Phys. Lett. B 451, 140-145 (1999), eprint hep-th/9803157.
  * Mészáros et al. (2019) P. Mészáros, D. B. Fox, C. Hanna, and K. Murase, Nature Rev. Phys. 1, 585 (2019), eprint 1906.10212.
  * Aartsen et al. (2018) M. G. Aartsen et al. (IceCube), Science 361, 147 (2018), eprint 1807.08794.


# Single-mode lasing based on $\mathcal{PT}$-breaking of two-dimensional
photonic higher-order topological insulator

Bofeng Zhu  Division of Physics and Applied Physics, School of Physical and
Mathematical Sciences,
Nanyang Technological University, Singapore 637371, Singapore  School of
Electrical and Electronic Engineering,
Nanyang Technological University, Singapore 637371, Singapore  Qiang Wang
Division of Physics and Applied Physics, School of Physical and Mathematical
Sciences,
Nanyang Technological University, Singapore 637371, Singapore  Yongquan Zeng
Electronic Information School,
Wuhan University, Wuhan 430072, China.  Qi Jie Wang<EMAIL_ADDRESS>Division of Physics and Applied Physics, School of Physical and Mathematical
Sciences,
Nanyang Technological University, Singapore 637371, Singapore  School of
Electrical and Electronic Engineering,
Nanyang Technological University, Singapore 637371, Singapore  Centre for
Disruptive Photonic Technologies, Nanyang Technological University, Singapore
637371, Singapore  Y. D. Chong<EMAIL_ADDRESS>Division of Physics and
Applied Physics, School of Physical and Mathematical Sciences,
Nanyang Technological University, Singapore 637371, Singapore  Centre for
Disruptive Photonic Technologies, Nanyang Technological University, Singapore
637371, Singapore

###### Abstract

Topological lasers are a new class of lasers that seek to exploit the special
properties of topological states of light. A typical limiting factor in their
performance is the existence of non-topological states with quality factors
comparable to the desired topological states. We show theoretically that, by
special design of uniform gain and loss on two sublattices of a two-
dimensional higher-order topological insulator (HOTI) lattice, single-mode
lasing based on topological corner states can be sustained over a wide dynamic
range of laser. This behavior stems from the phenomenon of parity/time-
reversal breaking of the topological corner states, where the undesired edge
and bulk states are suppressed much more effectively than through localized
pumping of the domain corners. These promising results shed light on the
opportunities for exploiting non-Hermitian phenomena in the design of high
performance topological lasers with small footprint.

Photonic topological insulators (PTIs) are an emerging class of photonic
devices based on adapting the topological band insulator concept, drawn from
condensed matter physics, to light [1, 2, 3, 4, 5]. PTIs host novel photonic
modes with unusual properties, such as robustness against disorder-induced
localization, that may be useful for compact and low-loss waveguiding (e.g.,
in optical [6, 7], quantum optical [8], and terahertz chips [9]), signal
processing (e.g., amplification and harmonic generation [10, 11]), lasing [12,
13, 14, 15, 16, 17, 18, 20, 21, 22, 23, 24, 25, 26], and other applications
[27, 28]. Of these, PTI-based lasers have drawn particular attention. They are
conceptually intriguing since lasers are intrinsically non-Hermitian and
nonlinear, whereas the theory of topological band insulators is formulated in
the Hermitian and linear regime, raising the possibility of new phenomena
arising from the interplay between band topology and laser physics [29, 30,
31, 32, 33, 34, 35, 36]. The first PTI lasers were based on one-dimensional
(1D) Su-Schrieffer-Heeger (SSH) lattices, with SSH topological boundary states
promoted to laser modes by the application of gain [12, 14, 15, 16].
Subsequently, PTI lasers were devised into two-dimensional (2D) lattices with
1D topological edge states [17, 18, 20], whose backscattering-resistant
features may allow for single-mode lasing at high powers [17, 18, 19]. Lasers
based on higher-order topological insulator (HOTI) lattices have also been
studied [23, 24, 25, 26]. Unlike standard or “first-order” topological
insulators, HOTIs manifest their topological states on higher-order boundaries
such as the “zero-dimensional” corners of 2D domains [37, 38]. Photonic HOTIs
in general [39, 40, 41, 42, 43, 44, 45], and HOTI lasers in particular [23,
24, 25, 26], are a promising avenue of research since they can be relatively
easy to fabricate (e.g., in the form of photonic crystal slabs) and to
integrate into existing on-chip photonic platforms.

PTI lasers of all types face the challenge of ensuring that the topological
edge or boundary states lase preferentially. There is usually no a priori
reason for topological states to have Q-factors higher than the the bulk
states of the lattice, so simple pumping schemes tend to induce multimode
lasing, with much of the emission being due to (mutually incoherent) bulk
laser modes. Previous 2D PTI laser designs have sought to address this problem
by selectively pumping the spatial regions where the topological states reside
[17, 18, 20]; in HOTI lasers, for instance, corner regions can be selectively
pumped to induce lasing of topological corner states [23, 24, 25, 26]. Another
approach is to use frequency-selective gain [33], but this is harder to
implement in semiconductor laser lattice-type devices. A third approach, which
has been employed in 1D SSH lasers [14, 15], is to apply gain to only part of
each unit cell (e.g., one of the two sites in each SSH unit cell); even if
this pumping scheme is applied uniformly throughout the lattice, the
asymmetric gain distribution can supply more effective gain to the topological
state than the bulk states [14, 15, 46]. Thus far, this last type of pumping
scheme has not been explored beyond 1D lattices.

In this paper, we show that in a HOTI laser containing gain and loss in
different sublattices, the lasing of topological states can be protected by
parity/time-reversal ($\mathcal{PT}$) breaking [46, 47, 48, 49]. The 2D
photonic HOTI lattice exhibits stable single-mode lasing over a substantial
span of pump strengths. The laser mode, which is composed of a superposition
of topological corner states determined by mirror symmetry and $\mathcal{PT}$
symmetry breaking, experiences higher effective gain than all competing edge
and bulk states, even without localized pumping of the domain corners. By
contrast, a standard localized pumping scheme (i.e., pumping the domain
corners and without exploiting $\mathcal{PT}$ breaking) induces multimodal
lasing, with the topological corner states not even having the lowest
threshold. The use of $\mathcal{PT}$ breaking to stabilize non-topological
single-mode lasing has previously been demonstrated [51, 52], and the present
work shows that single-mode lasing based on topological HOTI corner states can
be protected in a similar way, through the suppression of competing edge and
bulk states. These results point to further opportunities in making use of the
interplay between non-Hermitian effects and band topology, a subject of
considerable recent theoretical interest [56], in the design of future
topological laser devices.

Figure 1: (a) Schematic of the lattice. A hexagonal inner domain (yellow) is
surrounded by an outer domain (white). In the Hermitian and linear regime, the
inner domain is a higher order topological insulator (HOTI) in the nontrivial
phase, and the outer domain is trivial. In the close-up view (right panel),
dashes indicate one unit cell, thin (thick) lines correspond to intracell
(intercell) couplings $t_{1}=1-\delta$ ($t_{2}=1+\delta$), sites with gain
(loss) are shown as filled red (unfilled blue) circles, and the reflection
symmetry operations $\mathcal{P}$ and $\mathcal{P}^{\prime}$ are indicated.
(b) Spectrum of the linear lattice, plotted against the gain/loss parameter
$\gamma$ for $\delta=0.6$. Left panel: $\mathrm{Re}[E]$ versus $\gamma$. At
$\gamma=0$, there are six mid-gap corner states, which undergo $\mathcal{PT}$
breaking as $|\gamma|$ is increased. Right panel: zoomed-in plot of
$\mathrm{Re}[E]$ (solid curves) and $\mathrm{Im}[E]$ (dashes) versus $\gamma$.
The branches labelled $\mathrm{I}$ and $\mathrm{II}$ are one- and two-fold
degenerate respectively. Square and star markers indicate the modes in Fig. 2.
The number of non-real eigenvalues in each $\gamma$ range is shown in
parentheses.

Our analysis is performed using a tight-binding model containing nearest-
neighbor couplings, with each site representing an optical resonator (a
similar approach has been successfully used to model real PTI lasers [14, 17,
24, 26]). The model is both non-Hermitian and nonlinear: the pumped lattice
sites are modelled with saturable gain dependent on the local site intensity,
while the remaining lattice sites are subject to linear fixed loss
(representing a combination of outcoupling loss and material absorption). The
above-threshold behavior of the laser is studied using both frequency-domain
eigenmode analysis and time-domain simulations.

Figure 2: (a) Phase diagram showing the number of $\mathcal{PT}$-broken
eigenstates versus $\gamma$ and $\delta$. The white dot indicates the
parameter choice $\gamma=0.012$ and $\delta=0.6$ used in the other subplots.
(b)–(c) Intensity and phase distribution of selected eigenstates,
corresponding to the markers in Fig. 1(b). The area of the circles is
proportional to the intensity, and the color indicates the phase. (b)
Eigenstate on the non-degenerate branch I, which is $\mathcal{PT}$ symmetric
(real $E$) and antisymmetric under $\mathcal{P}^{\prime}$. (c) Eigenstates on
the degenerate branch II, which are $\mathcal{PT}$ broken with
$\mathrm{Im}(E)>0$ (their damped $\mathcal{PT}$ partners are not shown), and
$\mathcal{P}^{\prime}$-symmetric (left) and
$\mathcal{P}^{\prime}$-antisymmetric (right) respectively. Figure 3: (a) Plot
of $\mathrm{Im}[E]$ versus pump parameter $\Gamma$ for the corner states, in
the linear regime (no gain saturation). Lasing occurs upon crossing
$\mathrm{Im}[E]=0$ (horizontal dashes); the above-threshold regime is marked
in yellow. Eigenvalues with gain and loss are respectively marked by red and
blue. In this and subsequent subplots, we take $\gamma_{0}=0.005$,
$\gamma_{B}=0.017$, and $\delta=0.6$. (b) Time dependence of the total
intensity $I_{\mathrm{tot}}=\sum_{n}|\psi_{n}|^{2}$, with $\Gamma=0.012$ and
gain saturation enabled. Each site is initialized to $(\alpha+i\beta)f_{0}$,
where $\alpha,\beta$ are drawn from the standard normal distribution and
$f_{0}$ is a scale factor; results are shown for three independent initial
conditions with different $f_{0}$, showing convergence to the same steady
state $I_{\mathrm{tot}}$. (c) Intensity distribution in the steady state. A
comparison with Fig. 2(c) shows that the lasing mode is a superposition of
$\mathcal{P}^{\prime}$-symmetric and $\mathcal{P}^{\prime}$-antisymmetric
modes with broken $\mathcal{PT}$. (d) Time-averaged total intensity
$\left\langle I_{\mathrm{tot}}\right\rangle$ versus $\Gamma$. The kink at
$\Gamma\approx 0.022$ corresponds to the onset of multimode lasing. The values
of $\left\langle I_{\mathrm{tot}}\right\rangle$ are obtained by averaging over
$t\in[9000,10000]$. (e) Time dependence of $I_{\mathrm{tot}}$ in the multimode
regime, $\Gamma=0.03$. (f) Snapshot of the intensity distribution in the
multimode regime, taken at $t=5000$ with the same parameters as in (e).

The model, depicted in Fig. 1(a), consists of a Wu-Hu model [57] with
gain/loss added to alternate sites. The lattice consists of six sites per unit
cell; the intracell (intercell) couplings are $t_{1}=1-\delta$
($t_{2}=1+\delta$), where $\delta$ is a real tuning parameter. In the inner
(outer) domain, we set $\delta>0$ ($\delta<0$). Hence, in the absence of
gain/loss, the inner domain constitutes a higher-order topological insulator
(HOTI) in the nontrivial topological phase [42, 45], surrounded by a
topologically trivial cladding, with topological corner states occurring at
each of the six corners of the domain (the higher-order topology in the model
is tied to the “filling anomaly” of the $C_{6v}$-symmetric topological
crystalline insulator lattice [42, 45]). Next, we subdivide the lattice sites
into two sublattices, A and B, and assign them gain and loss respectively.
These sites are shown as red and blue circles in the right panel of Fig. 1(a).
In the position representation, the Hamiltonian is

$H=\left\\{\sum_{r\in
A}\sum_{i=1}^{3}t_{r,s_{i}}a_{r}^{\dagger}b_{r+s_{i}}+\mathrm{h.c.}\right\\}\\\
+\sum_{r\in A}i\gamma_{A}(r)\;a_{r}^{\dagger}a_{r}-\sum_{r\in
B}i\gamma_{B}\,b_{r}^{\dagger}b_{r},$ (1)

where $a_{r}$ ($b_{r}$) denotes the annihilation operator on the A (B) site at
position $r$; $s_{1,2,3}$ are the nearest neighbor displacement vectors shown
in Fig. 1(b); $t_{r,s_{i}}\in[t_{1},t_{2}]$ are nearest neighbor hoppings
which are assigned as stated above; $\gamma_{A}(r)$ is a gain term on each A
site, which may be $r$-dependent due to gain saturation; and $\gamma_{B}$ is a
loss term for the B sites. Each eigenvalue of $H$, denoted by $E$, represents
a detuning relative to a principal frequency.

We first analyze the system in the absence of gain saturation and with
balanced gain/loss $\gamma_{A}=\gamma_{B}=\gamma$. The lattice is
$\mathcal{PT}$ symmetric [49], where $\mathcal{P}$ is a reflection around the
vertical axis as indicated in Fig. 1(a), and $\mathcal{T}$ is the complex
conjugation operation. Fig. 1(b) plots the eigenvalue spectrum against the
gain/loss parameter $\gamma$ for $\delta=0.6$. At $\gamma=0$, there are six
eigenstates (in two sub-bands) near $E=0$, which are the six HOTI corner
states; there is some hybridization-induced splitting due to the finite
separation between different corners. As $|\gamma|$ is increased, the corner
states successively undergo $\mathcal{PT}$ breaking [46, 47, 48, 49]: the
eigenvalues coalesce at $\mathrm{Re}[E]=0$ and then become purely imaginary,
as shown in the right panel of Fig. 2(b). The two “inner” branches, each two-
fold degnerate and labelled II, break $\mathcal{PT}$ first; the non-degenerate
“outer” branches, labelled I, break $\mathcal{PT}$ at a larger value of
$|\gamma|$.

Notably, the lattice’s other eigenstates (i.e., bulk and edge states) have
unbroken $\mathcal{PT}$ symmetry over the range of $\gamma$ shown here, and
their eigenvalues are real. In other words, the corner states undergo
$\mathcal{PT}$ breaking first. In the lasing model to be discussed below, this
is the key to producing laser modes from the topological corner states. The
$\mathcal{PT}$ breaking of the corner states depends on the effective coupling
between the corners; for instance, reducing $\delta$ makes the corner states
less localized and increases the critical value of $\gamma$. For larger
$|\gamma|$, the other eigenstates eventually undergo $\mathcal{PT}$ breaking,
acquiring complex eigenvalues. Fig. 2(a) shows the phase diagram for the
number of non-real eigenvalues versus $\gamma$ and $\delta$.

Aside from $\mathcal{PT}$ symmetry, the lattice is invariant under
$\mathcal{P}^{\prime}$, reflection across the horizontal axis, as indicated in
Fig. 1(a). The non-degenerate branch I is antisymmetric under
$\mathcal{P}^{\prime}$, and the two-fold degenerate branch II contains a
$\mathcal{P}^{\prime}$-symmetric and a $\mathcal{P}^{\prime}$-antisymmetric
mode. The symmetry or antisymmetry under $\mathcal{P}^{\prime}$ persists
through the $\mathcal{PT}$ breaking transition. Fig. 2 shows three exemplary
mode profiles at $\gamma=0.012$ and $\delta=0.6$, for which branch I is
$\mathcal{PT}$ symmetric and branch II is $\mathcal{PT}$ broken.

We now develop the tight-binding model into a laser model by introducing gain
saturation according to

$\gamma_{A}=\frac{\Gamma}{1+|\psi_{n}|^{2}}-\gamma_{0}.$ (2)

This describes a combination of saturable gain and background loss, where
$\Gamma$ parameterizes the pump strength (taken to be the same on all A
sites), $|\psi_{n}|^{2}$ is the intensity on site $n$, and $\gamma_{0}$ is a
background loss. Each B site is assigned a fixed loss $\gamma_{B}$. The losses
are due to some combination of outcoupling loss and material absorption, such
that the laser’s output power can be assumed to be proportional to the total
intensity $I_{\mathrm{tot}}=\sum_{n}|\psi_{n}|^{2}$.

Figure 4: (a) Plot of $\mathrm{Im}[E]$ versus pump parameter $\Gamma$ for all
eigenstates of a selectively pumped lattice, in the linear regime with no gain
saturation. The corner states are denoted by solid lines, and the areas marked
in orange and green indicate where the edge and bulk states reside (the
numerous individual edge and bulk states are not plotted for clarity). Inset:
schematic of the lattice with pumped sites marked by red circles. The other
sites all have loss $\gamma_{A}=\gamma_{B}=0.017$. In this and subsequent
subplots, we take $\gamma_{0}=0.005$ and $\delta=0.6$. (b) Time dependence of
the total intensity $I_{\mathrm{tot}}=\sum_{n}|\psi_{n}|^{2}$ for
$\Gamma=0.008$ (blue curve) and $\Gamma=0.012$ (brown curve), with gain
saturation enabled. Each site is initialized to $(\alpha+i\beta)f_{0}$, where
$\alpha,\beta$ are drawn from the standard normal distribution and
$f_{0}=0.01$. The right panels show snapshots of the intensity distributions
taken at $t=40000$.

In the linear limit $|\psi_{n}|^{2}\rightarrow 0$ (i.e., at threshold), this
describes a lattice with $\mathcal{PT}$ symmetric gain/loss
$\gamma=(\Gamma-\gamma_{0}+\gamma_{B})/2$ plus additional loss
$\gamma^{\prime}=(\gamma_{B}+\gamma_{0}-\Gamma)/2$ on all sites. Fig. 3(a)
shows the plot of $\mathrm{Im}[E]$ versus $\Gamma$ for the corner states. With
increasing $\Gamma$, the corner states undergo bifurcations that map to the
$\mathcal{PT}$ breaking transitions described above. The bulk and edge modes
map to $\mathcal{PT}$ unbroken modes and have damping rate
$\mathrm{Im}[E]=-\gamma^{\prime}$. The corner state bifurcations produce
branches with higher effective gain (and partner branches with lower effective
gain). The laser threshold occurs when the first curve reaches
$\mathrm{Im}[E]=0$.

Above threshold, the gain saturation in Eq. (2) must be accounted for. Fig.
3(b) shows time-domain simulations for the total intensity $I_{\mathrm{tot}}$
at pump strength $\Gamma=0.012$, above the lasing threshold. The initial
complex amplitude on each site $n$ is set to
$\psi_{n}(t=0)=(\alpha_{n}+i\beta_{n})f_{0}$, where $\alpha_{n}$ and
$\beta_{n}$ are drawn independently from the standard normal distribution and
$f_{0}$ is an overall scale factor. The three curves in Fig. 3(b) correspond
to different choices of $f_{0}$; they all converge to the same value of
$I_{\mathrm{tot}}$ at large times, indicating that the system operates as a
stable single-mode laser. The steady state intensity profile is shown in Fig.
3(c). Comparing this to Fig. 2(c), we conclude that the laser mode is a
superposition of the $\mathcal{P}^{\prime}$-symmetric and
$\mathcal{P}^{\prime}$-antisymmetric modes of branch II. We verified
numerically that the steady-state field distribution is an equal-amplitude
linear superposition of the two degenerate corner states shown in Fig. 2(c).

Fig. 3(d) shows the laser’s “input-output” curve of total intensity versus
pump strength $\Gamma$, obtained by time domain simulations. The state is
initialized with $f_{0}=0.01$, and the total intensity is averaged over a time
window $t\in[9000,10000]$, after transient oscillations have died away. The
lasing is single-mode up to about 2.4 times threshold ($\Gamma\approx 0.022$);
thereafter, the input-output curve develops a kink corresponding to the onset
of multimode lasing. Fig. 3(e) shows the complicated oscillations of
$I_{\mathrm{tot}}$ in this regime. Fig. 3(f) shows the intensity profile at a
typical instant, revealing that the extra laser modes arise from bulk states.
The corner state mode can still be clearly seen, but most of the output power
is emitted from the more numerous bulk sites.

To demonstrate the qualitative difference between the
$\mathcal{PT}$-stabilized HOTI laser and a conventional design with selected
pumping, Fig. 4 shows the performance of a HOTI laser with a spatially
selective pump. The passive part of the HOTI lattice is the same as before,
but all the sites near one domain corner are pumped, as illustrated in the
inset of Fig. 4(a). The plot of $\mathrm{Im}(E)$ versus pump strength
$\Gamma$, in the absence of gain saturation, is shown in Fig. 4(a). We see
that the edge states and bulk states acquire effective gain comparable to and
in some cases higher than the topological corner states. Fig. 4(b) shows the
time-dependent total intensity for two values of $\Gamma$; at $\Gamma=0.008$,
prior to the topological corner states reaching threshold, a mode evolving out
of the edge state bands has already lased, as shown in the lower right panel
of Fig. 4(b). The mode frequency, determined via Fourier transformation, is
$E=-0.57$, outside the band of topological corner states. At a slightly larger
pump, $\Gamma=0.012$, the lasing is multi-modal.

In conclusion, we have demonstrated that by applying gain and loss in a
$\mathcal{PT}$ symmetric manner within a HOTI lattice, the topological corner
modes can be made to lase preferentially in a single mode over a wide dynamic
range of the laser. The competing edge and bulk states have significantly
higher thresholds, even in the presence of spatial hole burning due to local
gain saturation, and absence of frequency selectivity in the gain medium. The
higher gain experienced by the topological laser mode arises from
$\mathcal{PT}$ symmetry breaking, which has previously been shown to stabilize
non-topological laser modes in a similar fashion [46, 47, 48, 49]. In future
work, it would be interesting to study how other non-Hermitian symmetry
breakings, and other inherently non-Hermitian phenomena, could be usefully
incorporated into topological laser devices.

This work was supported by the Singapore MOE Academic Research Fund Tier 3
Grant MOE2016-T3-1-006, Tier 2 Grant MOE2018-T2-1-176, Tier 2 Grant
MOE2016-T2-2-159, Tier 2 Grant MOE2019-T2-2-085, Tier 1 Grants RG187/18, and
National Research Foundation (NRF) Competitive Research Program (CRP) (NRF-
CRP18-2017-02 and NRF-CRP23-2019-0007).

## References

  * [1] L. Lu, J. Joannopoulos and M. Soljačić, Topological photonics, Nat. Photonics 8, 821 (2014).
  * [2] B. Qi, H. Chen, L. Ge, P. Berini and R. Ma, T. Ozawa, H. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. Rechtsman, D. Schuster, J. Simon, O. Zilberberg and I. Carusotto, Topological photonics, Rev. Mod. Phys. 91, 015006 (2019).
  * [3] F. Haldane and S. Raghu, Possible Realization of Directional Optical Waveguides in Photonic Crystals with Broken Time-Reversal Symmetry, Phys. Rev. Lett. 100, 013904 (2008).
  * [4] S. Raghu and F. Haldane, Analogs of quantum-Hall-effect edge states in photonic crystals, Phys. Rev. A 78, 033834 (2008).
  * [5] Z. Wang, Y. Chong, J. Joannopoulos and M. Soljačić, Observation of unidirectional backscattering-immune topological electromagnetic states, Nature 461, 772 (2009).
  * [6] M. Hafezi, E. A. Demler, M. D. Lukin, and J. M. Taylor, Robust optical delay lines with topological protection, Nat. Phys. 7, 907 (2011).
  * [7] M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor, Imaging topological edge states in silicon photonics, Nat. Photon. 7, 1001 (2013).
  * [8] S. Barik, A. Karasahin, C. Flower, T. Cai, H. Miyake, W. DeGottardi, M. Hafezi and E. Waks, A topological quantum optics interface, Science 359, 666 (2018).
  * [9] Y. Yang, Y. Yamagami, X. Yu, P. Pitchappa, J. Webber, B. Zhang, M. Fujita, T. Nagatsuma and R. Singh, Terahertz topological photonics for on-chip communication, Nat. Photon. 14, 446 (2020).
  * [10] V. Peano, M. Houde, F. Marquardt, and A. A. Clerk, Topological Quantum Fluctuations and Traveling Wave Amplifiers, Phys. Rev. X 6, 041026 (2016).
  * [11] D. Smirnova, S. Kruk, D. Leykam, E. Melik-Gaykazyan, D.-Y. Choi, and Y. Kivshar, Third-Harmonic Generation in Photonic Topological Metasurfaces, Phys. Rev. Lett. 123, 103901 (2019).
  * [12] P. St-Jean, V. Goblot, E. Galopin, A. Lemaître, T. Ozawa, L. Le Gratiet, I. Sagnes, J. Bloch, and A. Amo, Lasing in topological edge states of a one-dimensional lattice, Nat. Photon. 11, 651 (2017).
  * [13] S. Longhi, Non-Hermitian gauged topological laser arrays, Ann. Phys. (Berl. ) 530, 1800023 (2018).
  * [14] H. Zhao, P. Miao, M. H. Teimourpour, S. Malzard, R. E. Ganainy, H. Schomerus and L. Feng, Topological hybrid silicon microlasers, Nat. Commun. 9, 981 (2018).
  * [15] M. Parto et al., Edge-Mode Lasing in 1D Topological Active Arrays, Phys. Rev. Lett. 120, 113901 (2018).
  * [16] Y. Ota, R. Katsumi, K. Watanabe, S. Iwamoto and Y. Arakawa, Topological photonic crystal nanocavity laser, Commun. Phys. 1, 1 (2018)
  * [17] G. Harari, Miguel. A. Bandres, Y. Lumer, M. C. Rechtsman, Y. D. Chong, M. Khajavikhan, D. N. Christodoulides, M. Segev, Topological insulator laser: theory, Science 359, 6381 (2018).
  * [18] M. A. Bandres, S. Steffen, G. Harari, M. Parto, J. Ren, M. Segev, D. N. Christodoulides, M. Khajavikhan Topological insulator laser: Experiments, Science 359, 6381 (2018).
  * [19] M. Parto, Y. Liu, B. Bahari, M. Khajavikhan and D. Christodoulides, Non-Hermitian and topological photonics: optics at an exceptional point, Nanophotonics 10, 403 (2021).
  * [20] Y. Zeng, U. Chattopadhyay, B. Zhu, B. Qiang, J. Li, Y. Jin, L. Li, A. G. Davies, E. H. Linfield, B. Zhang, Y. D. Chong and Q. J. Wang, Electrically pumped topological laser with valley edge modes, Nature 578, 246 (2020).
  * [21] X. Gao, L. Yang, H. Lin, L. Zhang, J. Li, F. Bo, Z. Wang and L. Lu, Dirac-vortex topological cavities, Nat. Nanotechnol. 15, 1012 (2020).
  * [22] Z. Yang, E. Lustig, G. Harari, Y. Plotnik, Y. Lumer, M. A. Bandres and M. Segev, Mode-Locked Topological Insulator Laser Utilizing Synthetic Dimensions, Phys. Rev. X 10, 011059 (2020)
  * [23] W. Zhang, X. Xie, H. Hao, J. Dang, S. Xiaom, S. Shi, H. Ni, Z. Niu, C. Wang, K. Jin, X. Zhang and X. Xu, Low-threshold topological nanolasers based on the second-order corner state, Light. Sci. Appl. 9, 109 (2020).
  * [24] Z. Shao, H. Chen, S. Wang, X. Mao, Z. Yang, S. Wang, X. Wang, X. Hu and R. Ma, A high-performance topological bulk laser based on band-inversion-induced reflection, Nat. Nanotechnol. 15, 67, (2020)
  * [25] W. Noh, H. Nasari, H. Kim, Q. Van, Z. Jia, C. Huang and B. Kanté, Experimental demonstration of single-mode topological valley-Hall lasing at telecommunication wavelength controlled by the degree of asymmetry, Opt. Lett. 45, 4108 (2020).
  * [26] D. Smirnova, A. Tripathi, S. Kruk, M. Hwang, H. Kim, H. Park and Y. Kivshar, Room-temperature lasing from nanophotonic topological cavities, Light. Sci. Appl. 9, 127 (2020).
  * [27] D. Smirnova, D. Leykam, Y. D. Chong and Y. Kivshar, Nonlinear topological photonics, Appl. Phys. Rev. 7, 021306 (2020).
  * [28] Y. Ota, K. Takata, T. Ozawa, A. Amo, Z. Jia, B. Kante, M. Notomi, Y. Arakawa and S. Iwamoto, Active topological photonics, Nanophotonics 9, 547 (2020).
  * [29] S. Malzard and H. Schomerus, Nonlinear mode competition and symmetry-protected power oscillations in topological lasers, New J. Phys. 20, 063044 (2018).
  * [30] S. Longhi, Y. Kominis, and V. Kovanis, Presence of temporal dynamical instabilities in topological insulator lasers, Europhys. Lett. 122, 14004 (2018).
  * [31] Y. V. Kartashov and D. V. Skryabin, Two-Dimensional Topological Polariton Laser, Phys. Rev. Lett. 122, 083902 (2019).
  * [32] M. Seclì, M. Capone and I. Carusotto, Theory of chiral edge state lasing in a two-dimensional topological system, Phys. Rev. Res. 1, 033148 (2019).
  * [33] M. Seclì T. Ozawa, M. Capone and I. Carusotto, Spatial and Spectral Mode-Selection Effects in Topological Lasers with Frequency-Dependent Gain, arXiv:2012.10709.
  * [34] I. Amelio and I. Carusotto, Theory of the Coherence of Topological Lasers, Phys. Rev. X 10, 041060 (2020).
  * [35] P. Zapletal, B. Galilo and A. Nunnenkamp, Long-lived elementary excitations and light coherence in topological lasers, Optica 7, 1045 (2020).
  * [36] S. Gündoğdu, J. Thingna and D. Leykam, Edge mode bifurcations of two-dimensional topological lasers, Opt. Lett. 45, 3673 (2020).
  * [37] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Quantized electric multipole insulators, Science 357, 61 (2017).
  * [38] F. Schindler, A. M. Cook, M. G. Vergniory, Z. J. Wang, S. S. P. Parkin, B. A. Bernevig, and T. Neupert, Higher-order topological insulators, Sci. Adv. 4, eaat0346 (2018).
  * [39] M. Serra-Garcia, V. Peri, R. Süsstrunk, O. Bilal, T. Larsen, L. Villanueva and S. Huber, Observation of a phononic quadrupole topological insulator, Nature 555, 342 (2018).
  * [40] S. Mittal, V. Orre, G. Zhu, M. Gorlach, A. Poddubny and M. Hafezi, Photonic quadrupole topological phases, Nat. Photonics 13, 692 (2019).
  * [41] X. Zhou, Z. Lin, W. Lu, Y. Lai, B. Hou and J. Jiang, Twisted quadrupole topological photonic crystals, Laser Photonics Rev. 14, 2000010, 2020.
  * [42] W. A. Benalcazar, T. H. Li and T. L. Hughes, Quantization of fractional corner charge in ${C}_{n}$-symmetric higher-order topological crystalline insulators, Phys. Rev. B 99, 245151 (2019).
  * [43] X. Chen, W. Deng, F. Shi, F. Zhao, M. Chen and J. Dong, Direct observation of corner states in second-order topological photonic crystal slabs, Phys. Rev. Lett. 122, 233902 (2019).
  * [44] A. Hassan, F. Kunst, A. Moritz, G. Andler, E. Bergholtz and M. Bourennane, Corner states of light in photonic waveguides, Nat. Photonics 13, 697 (2019).
  * [45] B. Y. Xie, G. X. Su, H. F. Wang, F. Liu, L. M.. Hu, S. Y. Yu, P. Zhan, M. H. Lu, Z. L. Wang and Y. F. Chen, Higher-order quantum spin Hall effect in a photonic crystal, Nat. Commun. 11, 3768 (2020).
  * [46] H. Schomerus, Topologically protected midgap states in complex photonic lattices, Opt. Lett. 38, 1912 (2013).
  * [47] S. Longhi, Parity-time symmetry meets photonics: A new twist in non-Hermitian optics, EPL (Europhysics Letters) 6, 64001 (2017).
  * [48] L. Feng, R. El-Ganainy and L. Ge, Non-Hermitian photonics based on parity-time symmetry, Nat. Photon. 11, 752 (2017).
  * [49] R. El-Ganainy, K. G. Makris, M. Konstantinos, Z. H. Musslimani, S. Rotter and D. N. Christodoulides, Non-Hermitian physics and PT symmetry, Nat. Phys. 14, 11 (2018).
  * [50] M. Miri, P. LiKamWa and D. N. Christodoulides, Large area single-mode parity-time-symmetric laser amplifiers, Opt. Lett. 37, 764 (2012).
  * [51] L. Feng, Z. J. Wong, R. M. Ma, Y. Wang, X. Zhang, Single-mode laser by parity-time symmetry breaking, Science 346, 972 (2014).
  * [52] H. Hodaei, M. A. Miri, M. Heinrich, D. N. Christodoulides and M. Khajavikhan, Parity-time–symmetric microring lasers, Science 346, 975 (2014).
  * [53] M. H. Teimourpour, M. Khajavikhan, D. N. Christodoulides, and R. El-Ganainy, Robustness and mode selectivity in parity-time (PT) symmetric lasers, Sci. Rep. 7, 10756 (2017).
  * [54] W. Liu, M. Li, R. S. Guzzon, E. J. Norberg, J. S. Parker, M. Lu, L. A. Coldren and J. Yao, An integrated parity-time symmetric wavelength-tunable single-mode microring laser, Nat. Commun. 8, 1 (2017).
  * [55] R. Yao, C. S. Lee, S. Chi, V. Podolskiy and W. Guo, Electrically injected parity time-symmetric single transverse-mode lasers, Laser Photon. Rev. 13, 1800154 (2019).
  * [56] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Exceptional topology of non-Hermitian systems, Rev. Mod. Phys. 93, 015005 (2021).
  * [57] W. L. Hua and H. Xiao, Scheme for Achieving a Topological Photonic Crystal by Using Dielectric Material, Phys. Rev. Lett. 114, 223901 (2015).