# Energy levels for $\mathcal{PT}$-symmetric deformation of the Mathieu
equation

E. Cavalcanti<EMAIL_ADDRESS>Centro Brasileiro de Pesquisas Físicas/MCTI,
22290-180 Rio de Janeiro, RJ, Brazil  N.M. Alvarenga<EMAIL_ADDRESS>Instituto de Física, Universidade do Estado do Rio de Janeiro, 20559-900 Rio
de Janeiro, RJ, Brazil  F. Reis  Departamento de Matemática Aplicada,
Universidade Federal do Espírito Santo, 29932-540 Campus São Mateus, ES,
Brazil  J.R. Mahon  Instituto de Física, Universidade do Estado do Rio de
Janeiro, 20559-900 Rio de Janeiro, RJ, Brazil  C.A. Linhares  Instituto de
Física, Universidade do Estado do Rio de Janeiro, 20559-900 Rio de Janeiro,
RJ, Brazil  J.A. Lourenço  Departamento de Ciências Naturais, Universidade
Federal do Espírito Santo, 29932-540 Campus São Mateus, ES, Brazil

###### Abstract

We propose a non-Hermitian deformation of the Mathieu equation that preserves
$\mathcal{PT}$ symmetry and study its spectrum and the transition from
$\mathcal{PT}$-unbroken to $\mathcal{PT}$-broken phases. We show that our
model not only reproduces behaviors expected by the literature but also
indicates the existence of a richer structure for the spectrum. We also
discuss the influence of the boundary condition and the model parameters in
the exceptional line that marks the $\mathcal{PT}$ breaking.

##  I Introduction

The study of $\mathcal{PT}$-symmetric systems was originally born of pure
formal curiosity. This category of models is defined by the invariance of the
Hamiltonian over the successive application of a parity operator
($\mathcal{P}$) and a time-reversal operator ($\mathcal{T}$). A remarkable
property of $\mathcal{PT}$-invariant systems is that they can have real energy
levels even though the Hamiltonian is non-Hermitian. Not only that,
$\mathcal{PT}$-symmetric Hamiltonians are suited to investigate the boundary
between closed and open systems. This is indicated by inspection of the
spectrum: we can have a $\mathcal{PT}$-symmetric phase (where all the energy
levels real, closed scenario) and a $\mathcal{PT}$-broken phase (where some of
the energy levels become complex, a transition to an open scenario). Such
change is controlled by some deformation parameter in the Hamiltonian, which
occurs because the $\mathcal{PT}$ operator is non-linear. Therefore, although
it commutes with the Hamiltonian, there is no guarantee that they will share
the same eigenstates. This richness motivated a plethora of physical
applications in optics, superconductivity, microwave cavity, lasers,
electronic circuits, chaos, graphene, metamaterials,…, and the subject has
experienced an overwhelming growth in the past twenty years (see references in
[1, 2, 3, 4]).

The investigation of $\mathcal{PT}$-symmetric systems started with the study
of the spectrum of a deformed anharmonic oscillator [5] and the reality of the
spectrum was later rigorously proved [6]. A common behavior of such systems is
that, depending on the “strength” of the deformation, pairs of energy levels
join at a so-called exceptional point. This behavior has been observed since
the first examples were considered and it has repeatedly appeared in the
literature over the years. Around the same time, it seems that Ref. [7],
through the use of a non-Hermitian deformation of a periodic potential
($(i\sinh x)^{\alpha}\cosh^{\beta}(x)$), was the first to indicate the
presence of what the community nowadays call $\mathcal{PT}$ breaking, although
that work did not really discuss the emerging field of $\mathcal{PT}$
symmetry.

There is also a different kind of behavior, in which the energy levels form a
“energy loop”. This simply means that there is a lower and an upper bound
determining the region where the energy level is real. The appearance of such
a loop was first exhibited in Ref. [8] for the strong coupling regime of the
anharmonic oscillator with a linear deformation. Later on, now discussing from
the perspective of the growing field of $\mathcal{PT}$ symmetry, this loop
behavior was reported for periodic potentials in Refs. [9, 10, 11], which use
a Hamiltonian with a potential term of the type $ig\cos(\theta)$. The behavior
is also found for different periodic potentials as one can see in Refs. [12,
13].

In the present work we deal with a periodic potential as given by the Mathieu
equation, which has a high number of applications in both engineering and
physics, from classical phenomena (inverted pendulum, trap of charged
particles, radio frequency quadrupole,…[14]) to quantum phenomena (quantum
particle moving in a periodic potential, quantum pendulum,…). We intend to
explore the consequences of a non-Hermitian, but $\mathcal{PT}$-symmetric,
deformation on the Mathieu equation. In the past few years, different
approaches were considered in the literature which proposed some kind of
deformation [15, 16, 9, 10, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29, 28, 27]. In the context of $\mathcal{PT}$ symmetry, the first works to
discuss a periodic potential were Refs.[15, 16], which considered a Mathieu-
like potential $\cos(\lambda\phi)+i\varepsilon\sin(\lambda\phi)$, where the
strength of the deformation is controlled by the parameter $\varepsilon$.

Instead of being a mere theoretical curiosity, a deformed Mathieu potential
seems to appear in the study of graphene [25] and is directly used in the
context of optics [17, 30, 18, 19, 20, 31, 21, 22, 23]. Ref. [25] shows that
in order to describe charge carriers moving in curved/bent graphene one can
derive a generalized Dirac equation. This equation, when considering a
monochromatic sinusoidal bending of the graphene, can be recast as a Mathieu
equation with a complex driven parameter.

In the context of optical lattices, the electric-field envelope of the beam
(paraxial equation of diffraction) satisfies a Schrödinger-like equation [17,
18, 19] in which the refractive index takes the place of the ‘potential’. The
scenario with complex refractive index $n=n_{R}+n_{I}$ – where the imaginary
component $n_{I}$ represents gain or loss of energy – provides a natural
setting to perform experiments with non-Hermitian physics [17, 31].

##  II Deformed Mathieu Equation

Let us consider the classical Mathieu equation [32],

$u^{\prime\prime}(x)+\left(a-2q\cos(2x)\right)u(x)=0.$ (1)

In the language of quantum mechanics it can be taken as a time-independent
Schrödinger equation, where $a$ is related to the energy levels. It is a
differential equation associated to the Hamiltonian operator

$\mathcal{H}_{0}=p^{2}+2q\cos(2x).$ (2)

One of the simplest deformations we can introduce is to add a sine
contribution to the potential, so that we still have a periodic potential. To
keep some degree of generality we can consider both a parameter $\delta$,
controlling the strength of the deformation, and a parameter $j\in\mathbb{N}$
which gives the frequency of the new term (where $j=1$ is the simplest
scenario):

$\mathcal{H}=p^{2}+2q\left(\cos(2x)+i\delta\sin(2jx)\right).$ (3)

This operator is non-Hermitian but preserves $\mathcal{PT}$ symmetry as can be
readily checked111The parity operator $\mathcal{P}$ acts on positions as
$x\rightarrow-x$ and the time-reversal operator $\mathcal{T}$ acts as
$i\rightarrow-i$. That this sign change amounts to time reversal was shown by
Wigner a long time ago as a condition to preserve the Heisenberg algebra under
$\mathcal{PT}$ transformations [2].. The differential equation for our
deformed version of the Mathieu equation is then

$u^{\prime\prime}(x)+\left[a-2q\left(\cos(2x)+i\delta\sin(2jx)\right)\right]u(x)=0.$
(4)

The above equation is related to some developments in the literature. As
mentioned in the Introduction, in the context of curved graphene the authors
of Ref. [25] obtained an expression similar to our model for the choice
$j=1/2$, but in their setting the boundary conditions are (in our notation)
chosen at $x=0$ and $x=2\pi$ (see in Sec. III our choice of boundary
conditions). Also, Refs. [17, 18, 19, 30] considered a similar equation but
with a fixed value of $q$. We discuss a bit more about this at the end of Sec.
IV.

Let us comment here that for $j=1$ one can choose a transformation of the
coordinate parameter into the complex plane ($x=z+i\alpha$, with a fine-tuned
choice of $\alpha$) such that the non-Hermitian Hamiltonian turns into a
Hermitian one. The existence of a possible transformation relating non-
Hermitian and Hermitian Hamiltonians is well-known222It seems the relationship
between Hermitian and non-Hermitian Hamiltonians was studied by Dyson back in
the 1950s. and has already been considered for periodic potentials [33, 7,
16]. However, as far as we know, no physical explanation was yet provided for
it and it seems to be useful only from a mathematical perspective333This was
commented back in 2015 by C. M. Bender at a conference at the Israel Institute
for Advanced Studiesit seems that no substantial improvement has been done so
far.. Notice that after the transformation the new variable $z$ does not
belong to the physical (and real) coordinate space ($x\in\mathbb{R}$).

We remark that at $\delta=1,j=1$ the potential reduces to $2qe^{i2x}$. An
intuitive picture using classical orbits is provided by Ref. [34], where it is
shown that in this case classical orbits are periodic and open. According to
the established knowledge on $\mathcal{PT}$ breaking the energy levels are
real if the classical orbits are closed. However, one must notice that the
argument was produced for a different scenario. When discussing Mathieu-like
equations we are usually restricted to a “zone” given by the boundary
conditions. It then becomes indistinguishable whether the orbit is closed or
if it is open and the period has the size of the “zone”. This point of view
means that the choice $\delta=1,j=1$ still might produce a meaningful
behavior. Indeed, this can be seen for example in Ref. [24], where a band
structure for an exponential potential is discussed. Also, Ref. [35] discusses
how such a potential depends on the choice of the contour and produces an
example with a set of real eigenvalues.

In what follows we first consider the simpler scenario $j=1$ and discuss how
the deformation parameter $\delta$ affects the energy levels and the breaking
of $\mathcal{PT}$ symmetry. After this, we discuss the influence of $j\neq 1$,
and we observe what seems to be a new kind of behavior: the mixing of more
than two energy levels.

##  III Energy levels ($j=1$)

One of the simplest ways to get a feeling of how the eigenfunctions behave is
to approach the special case $q=0$ (simple harmonic oscillator), for which
periodic solutions exist only if $a=n^{2}$ (with $n=0,1,2,\dots$). To each
eigenvalue we can associate pairs of respectively even and odd solutions
$\\{\cos(nx),\sin(nx)\\}$. These can be $\pi$\- or $2\pi$-periodic functions.
For general $q$, the pairs of solutions which reduce to these as $q\rightarrow
0$ are written in terms of the Mathieu functions $ce_{n}(x,q)$ and
$se_{n}(x,q)$. Nonetheless, if $q\neq 0$, each of these eigenfunctions has a
distinct eigenvalue. Therefore, there is at most one periodic solution for
each eigenvalue. The usual notation labels the eigenvalues associated with the
even solutions $ce_{n}(x,q)$ as $a_{n}(q)$, with $n=0,1,2,\dots$, and those
associated with the odd solutions $se_{n}(x,q)$, with $n=1,2,\dots$, are
labelled by $b_{n}(q)$.

By definition, the eigenvalues of the classical Mathieu equation, Eq. (1), are
determined by the following boundary conditions [32]:

$\displaystyle u^{\prime}(0)$
$\displaystyle=u^{\prime}\left(\pi\right)=0,\qquad$ $\displaystyle a_{n}(q);$
(5a) $\displaystyle u(0)$ $\displaystyle=u\left(\pi\right)=0,\qquad$
$\displaystyle b_{n}(q).$ (5b)

That is, while the $a_{n}$ energy levels are given through the application of
Neumann boundary conditions at the endpoints, the $b_{n}$ energy levels are
obtained by Dirichlet ones. We shall employ the same boundary conditions to
deal with the deformed version of the Mathieu equation, Eq. (4).

In the scenario without deformation the usual Mathieu spectrum is formed by
the curves $a_{n}(q)$ and $b_{n}(q)$, separating regions of stability (regions
centered in the origin) and regions of instability. As soon as the deformation
parameter $\delta$ is turned on, what we observe is that we lose a small part
of the spectrum, the eigenvalues cease to be real for sufficiently negative
values of $q$ (for $\delta=0.1$, see Fig. 1, left). As we increase the value
of the deformation the allowed region diminishes (for $\delta=0.5$, see Fig.
1, middle) and, at some point, we start to see closed energy loops (for
$\delta=2$, see Fig. 1, right). The contour of these energy loops are the real
eigenvalues of the problem and inside this region, we expect stable solutions.

We can also look directly at the eigenvalues as a function of both $q$ and
$\delta$ by a 3D plot, see Fig. 2 so that we can easily see the smooth
transition from the usual Mathieu scenario to the deformed scenario where the
allowed energy levels describe loops enveloping a region of stability.

In Fig. 1 the smallest energy loop is from the smallest energy level. This
behavior should already be expected as in the classical Mathieu those regions
of stability increase in size when we increase the eigenvalue.

Figure 1: Energy levels for the deformed Mathieu equation. From left to right:
$\delta=0.1,0.5,2.0$. The dotted curves are the $a_{n},b_{n}$ curves from
classical Mathieu, the gray regions are regions of instability. Figure 2:
First two energy levels $b_{1},b_{2}$ from the $b_{n}$ set as a function of
both the $q$ parameter and the deformation $\delta$. We can see clearly the
behavior from the classical Mathieu (at $\delta=0$) and that the curves seem
to smoothly fold as the deformation parameter increases.

Notice that we lose reflection symmetry in the parameter $q$, in comparison
with the classical Mathieu scenario, but there is an invariance under a
reflection of the deformation parameter, that is,
$a_{n}(q,\delta)=a_{n}(q,-\delta)$ and $b_{n}(q,\delta)=b_{n}(q,-\delta)$.

With respect to the loss of reflection symmetry in $q$, one must recall that
this relationship was indeed between different sets of energy levels, that is
$a_{2k}(q)=a_{2k}(-q)$, $b_{2k}(q)=b_{2k}(-q)$, $a_{2k+1}(q)=b_{2k+1}(-q)$.
However, once we turn on the deformation parameter $\delta$ there is a mixing
between neighbor energy levels. This is best illustrated in Fig. 2, as we can
see the evolution of the levels $b_{1}$ and $b_{2}$. Already in the undeformed
scenario, they seem to asymptotically merge, then when we turn on the
deformation the point where they merge (the exceptional point) approaches
smaller values of $q$. In other words, the reflection symmetry in $q$ was not
really expected in our scenario as it would mean the merging of two
independent sets of energy levels ($a_{n}$ produced by Neumann b.c. and
$b_{n}$ produced by Dirichlet b.c.).

On the other hand, the reflection symmetry on the deformation seems to
indicate an overall balance of the system (at least its spectrum) between gain
and loss. If we look back to Eq. (3) and the boundary conditions Eq. (5) we
see that $\delta>0$ means that in the first-half region ($x\in[0,\pi/2]$) the
hamiltonian has a positive imaginary component and in the second-half region
($x\in[\pi/2,\pi]$) the Hamiltonian has a negative imaginary component. That
is, for $\delta>0$ our system first experiences gain and afterwards
experiences loss. When we take $\delta<0$ we change the ordering to
loss$\rightarrow$gain. This reasoning suggests that the spectrum of the system
does not distinguish whether it first experiences gain or loss.

##  IV PT transition ($j=1$)

Regarding the reality of the spectrum it is usual to define that the
parametric region where the eigenvalues are real is $\mathcal{PT}$ symmetric
and the parametric region where at least one eigenvalue is lost is
$\mathcal{PT}$ broken. Therefore, the $\mathcal{PT}$ transition occurs exactly
when one of the eigenvalues is lost. Our scenario is a bit richer due to the
existence of the $q$ parameter. We can, therefore, choose to discuss the
$\mathcal{PT}$ transition as the point where for each value of $q$ at least
one eigenvalue is lost. Recalling that the $\mathcal{PT}$ transition means a
transition from a closed to an open system, here we are describing a
transition from a scenario where the “ground state” stability is attainable to
a scenario where it is completely unstable. This might be tracked by looking
at the first eigenvalue, which has the smallest range of stability. So, for
each possible value of the deformation parameter $\delta$ we check the value
of $q$ ($q>0$) at which the eigenvalues become complex. This can be done
numerically and is exhibited in Fig. 3. We observe that there is a constraint
in the allowed value of the deformation $\delta$, above which the system
becomes $\mathcal{PT}$ broken. This upper bound is $\delta=1$ for almost all
values of $q$, which seems to be related to the condition we found earlier
when we compared the deformed and classical Mathieu equations. Only when this
parameter is small, large values of $\delta$ are allowed. This seems to be
related to some “resistance” of the system as it comes near the $q\rightarrow
0$ limit, which is related to the harmonic oscillator limit. We refer to
“resistance” in the sense that the system resists undergoing the
$\mathcal{PT}$ breaking.

Figure 3: Numerical evaluation of the exceptional line. In the filled region
the $\mathcal{PT}$ symmetry remains unbroken and all energy levels are real.
Outside this region at least one energy level is lost and $\mathcal{PT}$
symmetry is broken. The line of phase separation is the exceptional line,
where a coalescence of energy levels occur. This image has reflection symmetry
for $\delta\rightarrow-\delta$.

In Fig. 3 we completely exhibit the phase diagram for $j=1$ according to
variations of the parameters $\delta$ and $q$. There are two phases, the
unbroken $\mathcal{PT}$ phase being represented by the filled region inside
the exceptional lines, and the broken $\mathcal{PT}$ phase being outside the
exceptional lines. For negative values of $\delta$ we just apply the
reflection symmetry $\delta\rightarrow-\delta$. It is noticeable that the
behavior for positive and negative values of $q$ is distinct. While for $q<0$
the phase is already broken for small values of $\delta$, for $q>0$ we stay in
the unbroken phase as long as $\delta<1$.

We remark that if someone is interested in a system in which the parameter $q$
varies freely in the positive real axis, it makes sense to say that $\delta=1$
is an upper bound above which $\mathcal{PT}$ symmetry is broken. We can avoid
this ‘restriction’ only for small values of the parameter $q$, where larger
values of $\delta$ are still allowed.

In the literature, Refs. [17, 30, 18] discuss a potential
$4(\cos^{2}(x)+iV_{0}\sin(2x))$, while Ref. [19] rewrites it as
$4(\cos(2x)+2iV_{0}\sin(2x))$. According to them, at $V_{0}=1/2$ the
$\mathcal{PT}$ breaking occurs, but they consider a fixed value of the
potential amplitude. If we compare these results with our scenario, we see
that it still holds for a choice of $q=2$ for our parameter (this reproduces
the fixed amplitude from Refs. [17, 30, 18, 19]). However, we also generalize
it for different values of $q$.

Based on the discussion in the literature [17, 30, 18, 19] and the findings of
the present paper, we suggest an experimental measurement. We propose to
consider the effect of the deformation parameter $\delta$ on the
$\mathcal{PT}$ breaking when the amplitude $2q$ is reduced. As indicated in
Fig. 3 we expect that the smaller the amplitude the higher the “stability” of
the system. One can understand it intuitively by noticing that the limit of
small values of $q$ is equivalent to reducing the system to a harmonic
oscillator, which is stable. Therefore, this setting should be more
“resistant” against $\mathcal{PT}$ breaking.

##  V Different scenarios ($j\neq 1$)

Figure 4: Energy levels $a_{n}$ ($u^{\prime}(0)=u^{\prime}(\pi)=0$) for $j=2$.
(Upper left) For $\delta=0.33$, pairs of energy levels merge at the left
plane. (Upper right) For $\delta=0.34$, four energy levels merge, changing the
structure. (Bottom left) For $\delta=0.43$, an “isolated loop” appears at the
left plane, two energy levels become degenerate at a point in the right plane.
(Bottom right) For $\delta=0.44$, a new energy loop appears.

In what follows we only show the behavior of the energy levels $a_{n}$
produced by the Neumann boundary conditions, Eq. (5a). The second set of
energy levels produced by Dirichlet boundary conditions, Eq. (5b), has a
similar behavior and do not add any useful information at this moment. We
return to it at the end of the section to discuss the $\mathcal{PT}$-broken
regime.

In Fig. 4 we can see that the scenario $j=2$ starts similarly to the $j=1$
case, the energy levels join for $q<0$ for small values of the deformation
parameter $\delta$. However, the increase in the deformation parameter
$\delta$ which controls the strength of the imaginary contribution to the
Hamiltonian produces a richer structure for the spectrum.

With the increase of the deformation, we notice lines merging for $q>0$, as
can be seen for $\delta=0.33$ (upper left in Fig. 4), this is the energy loop
behavior commented in the previous section and already reported in many works
about $\mathcal{PT}$-symmetric models. However, we see a new behavior, where
four originally independent energy levels merge as a single line for
$\delta=0.34$ (See the upper right image in Fig. 4). Increasing a bit further
the deformation we also see the emergence of new loops (still connected to the
origin at $q=0$) and “isolated loops” (that are disconnected from the origin
$q=0$). Notice that around $\delta=0.43$ (bottom left in Fig. 4) two energy
levels approach each other and become degenerated. This degeneracy marks the
formation of the closed loop. Also, notice that for $\delta=0.44$ (bottom
right in Fig. 4) the system is well defined (the first six energy levels are
present) for small $q$, then at some point $q$ two of them disappear (this
means a transition from a closed to an open system) and after a small increase
in $q$ we recover all energy levels.

The behavior of the energy levels $a_{n}$ for $j=2$ tells us another different
aspect in comparison with $j=1$. Recall that for $j=1$ the breaking of
$\mathcal{PT}$ symmetry could be checked by looking only to the lowest energy
level (the composition of $a_{0},a_{1}$). This occurs because the lowest
energy levels produces the smallest energy loop (smallest values of the lower
and upper boundary $q$ given by the coalescence points). However, at $j=2$
this is not the case. We see clearly in Fig. 4 that at $\delta=0.44$ the first
energy level is not even limited from the right and there is already a loop
produced by joining the original 3rd and 4th energy levels. This indicates the
need to change our methodology to check the $\mathcal{PT}$-breaking.

Figure 5: The $a_{n}$ energy levels for $j=2$. At $\delta=0.76$ (Left) the
$\mathcal{PT}$-breaking occurs in the degeneracy joining the 3rd and 4th
energy levels , at $\delta=0.77$ (Right) a jump occurs and the
$\mathcal{PT}$-breaking occurs at coalescence point of the 2nd and 3rd energy
levels.

The breaking of $\mathcal{PT}$ symmetry is controlled by the loop produced by
the 3rd and 4th energy levels until this loop “collides” with the loop
produced by the 1st and 2nd energy levels. This occurs around $\delta=0.76$
(see Fig. 5, left) where a new degeneracy occurs. After this “collision” both
loops merge into a single energy loop with six coalescence points (see Fig. 5,
right). This behavior will be responsible for a jump when we try later on to
determine the dependence between $\delta$ and $q$. In the scenario exhibited
in Fig. 5 the jump occurs from the coalescence between the 3rd and 4th energy
levels to the coalescence between the 2nd and 3rd energy levels. We can also
see it in Fig. 7.

We remark that such discontinuity does not necessarily occur at any time at
which independent energy loops merge. For example, at $j=3$, $\delta=0.8$
(Fig. 6, left) the breaking point is signalled by the 3rd and 4th energy
levels, and it is still the case even after the first six energy levels are
merged into a single loop (Fig. 6, right).

Figure 6: The $a_{n}$ energy levels for $j=3$. At $\delta=0.80$ (left) we see
three closed loops produced by the first six energy levels. At $\delta=2$
(right) these energy levels merge into a single loop. Notice that in both
figures the $\mathcal{PT}$-breaking is controlled by the coalescence of the
3rd and 4th energy levels.

Instead of looking at just one energy level, as we did for $j=1$, we see,
based on the previous discussion, that we must check all the first energy
levels and identify a $\mathcal{PT}$-breaking if anyone of them gets complex.
For all scenarios the $a_{n}$ energy levels (Neumann b.c.) break first (small
values of $q$) than the $b_{n}$ energy levels (Dirichlet b.c.). This is why we
focus our attention first on the $a_{n}$ set. In Fig. 7 we see how the
exceptional value of $q$ varies both with $\delta$ and $j$. Notice that for
$j=2$ there is a discontinuous jump around $\delta\approx 0.76$: this is
exactly the jump we discuss on Fig. 5. These jumps also occur for different
choices of $j$ as we increase the observed range of possible values of the
parameter $q$.

Figure 7: Exceptional lines at which a coalescence of energy levels $a_{n}$
occurs, for different values of the parameter $j$. (Upper) First quadrant of
the $\delta\times q$, above each exceptional line we have a
$\mathcal{PT}$-broken phase and at least one energy level becomes complex.
(Bottom) Fourth quadrant of the $\delta\times q$; under each exceptional line
$\mathcal{PT}$ symmetry is broken. We do not show the 2nd and 3rd quadrants
because there is a reflection symmetry $\delta\rightarrow-\delta$ in the
figure.

The curve $q(\delta,j)$ represents what we call an “exceptional line”, and it
marks a transition between a phase of unbroken $\mathcal{PT}$ symmetry (under
each curve, all energy levels are real) and a phase of broken
$\mathcal{PT}$-symmetry (above each curve, at least one energy level is
complex). This transition also means a transition between closed and open
systems. It is important to notice that our diagram does not contemplate the
recovery of the energy levels, although we know this behavior occurs as shown
in Figs. 4, 5 and 6. This was an active choice as we are mostly interested in
the unbroken $\mathcal{PT}$ phase that is connected with the origin at $q=0$;
in a more complete phase diagram, “islands” of unbroken $\mathcal{PT}$ phases
are expected.

Regarding the “stability” of the system, in the sense that we have a bigger
parametric region $(q,\delta)$ where a $\mathcal{PT}$-symmetric phase
connected to the origin occurs, we can see two main behaviors. For $q>0$ and
small values of $\delta$, the increase of the factor $j$ means a loss of
“stability” as indicated in Fig. 7 by the diminished area under the
exceptional line. Anyway, we cannot extract any clear hierarchy between
different choices of the parameter $j$ due to the occurrence of the jumps. On
the other hand, both for large values of the deformation $\delta$ and for
negative values of the parameter $q$ we observe a clear hierarchy where the
increase in the parameter $j$ means an increase in the parametric region with
a $\mathcal{PT}$-symmetric phase.

For large values of $\delta$ the exceptional lines exhibit a smooth behavior
for all choices of the parameter $j$. However, although the only free
parameter is the value of $j$, it looks like the decaying behavior is not
related to the value of $j$. If we make a power law fitting of the type
$q(\delta,j)=A_{j}\delta^{-\alpha_{j}}$ for each exceptional line we get

$\displaystyle q(\delta,j=1)=0.72\,\delta^{-1.37},$ (6a) $\displaystyle
q(\delta,j=2)=1.07\,\delta^{-0.88},$ (6b) $\displaystyle
q(\delta,j=3)=1.92\,\delta^{-0.79},$ (6c) $\displaystyle
q(\delta,j=4)=3.03\,\delta^{-0.80}.$ (6d)

Figure 8: Power-law fit for the exception line $q(\delta,j)$ for different
choices of $j$. The markers are numerical points (square for $j=1$, circle for
$j=2$, triangle for $j=3$ and diamond for $j=4$) and the curve is the fitting
$q(\delta,j)=A_{j}\delta^{-\alpha_{j}}$. See Eqs. (6) for the values of
$A_{j}$ and $\alpha_{j}$.

However, as commented, the exponents of this power law do not seem to the
authors to have any relationship with the value of the parameter $j$.

For completeness, let us now show, see Fig. 9, the exceptional lines for the
second set of energy levels ($b_{n}$) produced by the Dirichlet boundary
conditions, Eq. (5b). In this scenario, we do not see any discontinuous jump
for $j=2$. Also, for all values of $j$ the region of unbroken $\mathcal{PT}$
is broader than the previous scenario (Fig. 7, Neumann boundary condition).
This means that for the deformed Mathieu equation the Dirichlet b.c. turns the
system more stable than the Neumann boundary condition. Intuitively, we
usually say that the Dirichlet b.c. means a closed boundary condition while
the Neumann b.c. is a open boundary condition. In this sense, it is very
natural that the Neumann b.c., as an open b.c., breaks $\mathcal{PT}$ first.

Figure 9: Exceptional lines at which a coalescence of energy levels $b_{n}$
occurs, for different values of the parameter $j$ in the $\delta\times q$
plane. In the filled region we have an unbroken-$\mathcal{PT}$ phase. Recall
that there is a reflection symmetry $\delta\rightarrow-\delta$ in the figure.

##  VI Conclusion

We studied in this article a non-Hermitian but $\mathcal{PT}$-symmetric
deformation of the Mathieu equation (which might be reproduced, for example,
in optical experiments), that has as free parameters the deformation strength
$\delta$ and a parameter $j$ related to the frequency of the imaginary
contribution. We managed to show that the $\mathcal{PT}$-breaking point
depends on the deformation parameter $\delta$ and we discussed a phase diagram
where an exceptional line signals the transition from the unbroken
$\mathcal{PT}$ phase to the broken $\mathcal{PT}$ phase. We also show that the
change of the parameter $j$ turns the structure of the spectrum much richer
than the simple $j=1$ scenario.

Concerning the transition from a closed system (unbroken $\mathcal{PT}$ phase)
to an open system (broken $\mathcal{PT}$) we see that:

  1. 1.

The increase of $q$ reduces “stability” and favors the transition. Recall that
$q=0$ would make our model return to the classic harmonic oscillator.

  2. 2.

The increase of the deformation parameter $\delta$ also diminishes “stability”
and favors the transition. This is expected, as this is the parameter which
controls the strength of the imaginary contribution to the Hamiltonian.

  3. 3.

The increase in the frequency, controlled by the parameter $j$ of the sine
function makes “stability” grow.

  4. 4.

Dirichlet boundary conditions produce a stabler set of energy levels than
Neumann boundary conditions, which suggests that boundary conditions play a
significant role.

These conclusions suggest possible experimental investigations to be performed
for physical scenarios which can be modelled by a Mathieu equation. Also, we
discuss in our model, there is an impact of both the choice of boundary
conditions and a parameter related to the frequencies. We indicate, from a
formal perspective, the importance of investigating the different choices of
$\mathcal{PT}$-symmetric extensions one can make for a selected model.

###### Acknowledgements.

E.C. thanks the Brazilian agency Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq) for financial support through the PCID-B
scholarship (Programa de Capacitação Institucional). N.M.A. was financed by
the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil
(CAPES) – Finance Code 001.

## References

  * Bender [2020] C. M. Bender, $\mathcal{PT}$-symmetric quantum field theory, J. Phys. Conf. Ser. 1586, 012004 (2020).
  * Bender [2016] C. M. Bender, $\mathcal{PT}$ symmetry in quantum physics: From a mathematical curiosity to optical experiments, Europhysics News 47, 17 (2016).
  * Coulais  _et al._ [2020] C. Coulais, R. Fleury, and J. van Wezel, Topology and broken hermiticity, Nature Physics 17, 9 (2020).
  * El-Ganainy  _et al._ [2018] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Non-Hermitian physics and $\mathcal{PT}$ symmetry, Nature Physics 14, 11 (2018).
  * Bender and Boettcher [1998] C. M. Bender and S. Boettcher, Real spectra in non-Hermitian Hamiltonians having $\mathcal{PT}$ symmetry, Phys. Rev. Lett. 80, 5243 (1998), arXiv:physics/9712001 .
  * Dorey  _et al._ [2001] P. Dorey, C. Dunning, and R. Tateo, Spectral equivalences, Bethe Ansatz equations, and reality properties in $\mathcal{PT}$-symmetric quantum mechanics, Journal of Physics A: Mathematical and General 34, 5679 (2001).
  * Fernández  _et al._ [1999] F. M. Fernández, R. Guardiola, J. Ros, and M. Znojil, A family of complex potentials with real spectrum, Journal of Physics A: Mathematical and General 32, 3105 (1999).
  * Delabaere and Pham [1998] E. Delabaere and F. Pham, Eigenvalues of complex Hamiltonians with $\mathcal{PT}$-symmetry. II, Physics Letters A 250, 29 (1998).
  * Bender and Kalveks [2011] C. M. Bender and R. J. Kalveks, Extending $\mathcal{PT}$ symmetry from Heisenberg algebra to E2 algebra, Int. J. Theor. Phys. 50, 955 (2011), arXiv:1009.3236 [hep-th] .
  * Bender and Jones [2012] C. M. Bender and H. F. Jones, WKB analysis of $\mathcal{PT}$-symmetric Sturm-Liouville problems, J. Phys. A 45, 444004 (2012), arXiv:1201.1234 [hep-th] .
  * Jones-Smith and Kalveks [2013] K. Jones-Smith and R. Kalveks, Vector Models in PT Quantum Mechanics, Int. J. Theor. Phys. 52, 2187 (2013), arXiv:1304.5692 [quant-ph] .
  * Meisinger and Ogilvie [2013] P. N. Meisinger and M. C. Ogilvie, $\mathcal{PT}$ Symmetry in Classical and Quantum Statistical Mechanics, Phil. Trans. Roy. Soc. Lond. A 371, 20120058 (2013), arXiv:1208.5077 [math-ph] .
  * Bender [2019] C. M. Bender,  _$\mathcal{PT}$ symmetry: In quantum and classical physics_ (World Scientific, 2019).
  * McLachlan [1947] N. W. McLachlan,  _Theory of Application of Mathieu Functions_ (Oxford University Press, 1947).
  * Bender  _et al._ [2005] C. M. Bender, H. F. Jones, and R. J. Rivers, Dual $\mathcal{PT}$-symmetric quantum field theories, Phys. Lett. B 625, 333 (2005), arXiv:hep-th/0508105 .
  * Bender [2007] C. M. Bender, Making sense of non-Hermitian Hamiltonians, Rept. Prog. Phys. 70, 947 (2007), arXiv:hep-th/0703096 .
  * Makris  _et al._ [2008] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, Beam dynamics in $\mathcal{PT}$ symmetric optical lattices, Physical Review Letters 100, 103904 (2008).
  * Midya  _et al._ [2010] B. Midya, B. Roy, and R. Roychoudhury, A note on the invariant periodic potential, Physics Letters A 374, 2605 (2010).
  * Jones [2011] H. F. Jones, Use of equivalent Hermitian Hamiltonian for $\mathcal{PT}$-symmetric sinusoidal optical lattices, Journal of Physics A: Mathematical and Theoretical 44, 345302 (2011).
  * Lin  _et al._ [2011] Z. Lin, H. Ramezani, T. Eichelkraut, T. Kottos, H. Cao, and D. N. Christodoulides, Unidirectional invisibility induced by $\mathcal{PT}$-symmetric periodic structures, Physical Review Letters 106, 213901 (2011).
  * Graefe and Jones [2011] E.-M. Graefe and H. F. Jones, $\mathcal{PT}$-symmetric sinusoidal optical lattices at the symmetry-breaking threshold, Physical Review A 84, 013818 (2011).
  * Moreira  _et al._ [2013] F. C. Moreira, V. V. Konotop, and B. A. Malomed, Solitons in $\mathcal{PT}$-symmetric periodic systems with the quadratic nonlinearity, Physical Review A 87, 013832 (2013).
  * Brandão [2019] P. Brandão, Stability boundaries of a Mathieu equation having $\mathcal{PT}$ symmetry, Physics Letters A 383, 3043 (2019).
  * Longhi [2010] S. Longhi, Spectral singularities and Bragg scattering in complex crystals, Physical Review A 81, 022102 (2010).
  * Kerner and Naumis [2012] R. Kerner and G. Naumis, Interaction of flexural phonons with electrons in graphene: A generalized Dirac equation in corrugated surfaces, Physica B 407, 2002 (2012), arXiv:1102.0746 [cond-mat.mtrl-sci] .
  * Ge  _et al._ [2014] L. Ge, M. Shen, C. Ma, T. Zang, and L. Dai, Gap solitons in $\mathcal{PT}$-symmetric optical lattices with higher-order diffraction, Optics express 22, 29435 (2014).
  * Longhi [2019a] S. Longhi, Metal-insulator phase transition in a non-Hermitian Aubry-André-Harper model, Physical Review B 100, 125157 (2019a).
  * Longhi [2019b] S. Longhi, Topological phase transition in non-Hermitian quasicrystals, Physical Review Letters 122, 237601 (2019b).
  * Xu  _et al._ [2019] T. Xu, Y. Chen, M. Li, and D.-X. Meng, General stationary solutions of the nonlocal nonlinear Schrödinger equation and their relevance to the $\mathcal{PT}$-symmetric system, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 123124 (2019).
  * Musslimani  _et al._ [2008] Z. Musslimani, K. G. Makris, R. El-Ganainy, and D. N. Christodoulides, Optical solitons in $\mathcal{PT}$periodic potentials, Physical Review Letters 100, 030402 (2008).
  * Rüter  _et al._ [2010] C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, Observation of parity–time symmetry in optics, Nature physics 6, 192 (2010).
  * Hille [1997] E. Hille,  _Ordinary differential equations in the complex domain_ (Courier Corporation, 1997).
  * Ciafaloni [1978] M. Ciafaloni, Instanton contributions in reggeon quantum mechanics, Nucl. Phys. B 146, 427 (1978).
  * Bender  _et al._ [2014] C. M. Bender, D. W. Hook, N. E. Mavromatos, and S. Sarkar, Infinite class of $\mathcal{PT}$-symmetric theories from one timelike Liouville Lagrangian, Physical Review Letters 113, 231605 (2014), arXiv:1408.2432 [hep-th] .
  * Cannata  _et al._ [1998] F. Cannata, G. Junker, and J. Trost, Schrödinger operators with complex potential but real spectrum, Physics Letters A 246, 219 (1998).


# Optimized FLUKA cross sections for cosmic-ray propagation studies

P. De La Torre Luquea M. N. Mazziottab a The Oskar Klein Centre, Department of
Physics, Stockholm University, AlbaNova
SE-10691 Stockholm, Sweden
bIstituto Nazionale di Fisica Nucleare, Sezione di Bari, via Orabona 4,
I-70126 Bari, Italy

###### Abstract

The current great precision on cosmic-ray (CR) spectral data allows us to
precisely test our simple models on propagation of charged particles in the
Galaxy. However, our studies are severely limited by the uncertainties related
to cross sections for CR interactions. Therefore we have developed a new set
of cross sections derived from the FLUKA Monte Carlo code, which is optimized
for the treatment of CR interactions. We show these cross sections and the
main results on their application for CR propagation studies. Finally, we
discuss the prediction of a low-energy break in the electrons spectra inferred
from gamma-ray data.

##  1 Introduction

Galactic CRs, injected and accelerated in astrophysical sources, propagate
throughout the Galaxy for millions of years, occasionally interacting with the
gas in the interstellar medium (ISM) through spallation reactions that produce
secondary particles, such as gamma rays, neutrinos or secondary nuclei that we
call secondary CRs. The amount of secondary CRs formed provides direct
information about the mean grammage traversed by the so-called primary CRs
(i.e. the amount of gas per unit area that primary CRs cross during their
journey) and, therefore, about the time that these CRs reside in the Galaxy.
Measuring the spectra of secondary CRs such as Li, B, Be or F provides further
information on how CRs interact with the magneto-hydrodynamical (MHD)
turbulence in the ISM plasma and the energy dependence of these interactions.

The transport of CRs in the Galaxy is conventionally studied as a diffusive
process characterised by a diffusion coefficient which is, basically, a power-
law in energy ($D\propto E^{\delta}$) and whose spectral index $\delta$ is
intimately related to the MHD interactions that CRs undergo. The ratios of the
fluxes of secondary CRs to primary CRs, such as the B/C ratio, are used to
determine the diffusion coefficient, given that they are directly related to
each other, as can be seen from the approximate relation
$N_{sec}/N_{prim}(E)\propto\sigma(E)/D(E)$, where $\sigma(E)$ is the
spallation cross section of production of the secondary CR involved.
Therefore, a precise evaluation of these ratios is crucial nowadays to unveil
the process of transport of charged particles in the Galaxy and its main
features, which also requires a good knowledge of cross sections of CR
interactions in a broad energy range. However, while current CR data is
extremely precise since the past decade (at the level of a few percent),
recent studies have proved the need of improving cross sections on spallation
reactions in order to reduce the uncertainty related to the determination of
the diffusion coefficient from the secondary-to-primary flux ratios [1, 2, 3].
In fact, the accuracy of spallation cross sections measurements for the energy
range between the GeV and the TeV is very poor (20%) and frequently there is
not data at all for the energies of interest. On top of this, given that the
nuclear models describing these kind of interactions are mainly adjusted to
reproduce accelerator data, the cross sections computed from fundamental
models are not totally consistent with the available data in the GeV-TeV range
and, hence, we usually rely on parameterisations, fitted to the very scarce,
limited and uncertain experimental data, instead.

Nevertheless, these simulation codes, based on Monte Carlo simulations and
event generators experienced a positive boost in the last years, usually
driven by radiological and medical applications, which need to accurately
describe the transport of ions in different materials and their interactions.
Recently, the FLUKA Monte Carlo nuclear code 111http://www.fluka.org/fluka.php
has been optimized to be used in different kinds of astrophysical studies,
with special attention to CR interactions [4, 5, 6, 7]. FLUKA is a general
purpose tool that can be used to transport particles in arbitrarily complex
geometries and magnetic fields and study their nuclear interactions with
hadrons and nuclei from the MeV/n up to $10$ PeV.

Here, we overview the main results presented in De la Torre et al., (2022)
[8], where the full cross-section network for CR interactions up to iron
nuclei (Z=26) was computed using FLUKA. These cross sections are tested
against data and implemented in a customised version 222A similar version of
this code can be downloaded at https://doi.org/10.5281/zenodo.4461732 of the
DRAGON2 code [9, 10] in order to study the different ratios of the secondary
CRs B, Be and Li and the propagation parameters derived from these predictions
through a Markov chain Monte Carlo (MCMC) analysis [11] of the most recent
data from the AMS-02 collaboration [12]. Furthermore, we demonstrate that
gamma-ray data indicate the need of including a break at low energies ($<10$
GeV) in the injection spectrum of electrons.

##  2 The derived FLUKA cross sections

Inelastic and inclusive cross sections of all stable isotopes from protons to
iron impinging on helium and hydrogen as targets (representing the composition
of the ISM gas,) from $1$ MeV/n to $35$ TeV/n, are calculated with FLUKA using
176 bins equally spaced in a logarithmic scale.

Inelastic cross sections are those related to the probability of destruction
of a nucleus when it interacts with another particle and its effects in the
spectrum of CRs become only significant at low energies, when the time-scale
of inelastic collisions ($\tau_{inel}^{-1}\sim vn_{ISM}\sigma_{inel}(E)$,
where $v$ is the speed of the particle, $n_{ISM}$ is the target number density
and $\sigma_{inel}$ is the inelastic cross section of the interaction) is
smaller than the diffusion time-scale ($\tau\sim H^{2}/2D(E)$, where $H$ is
Galactic halo size and $D(E)$ is the diffusion coefficient). Experimental data
on inelastic cross sections with proton as target are measured with a
precision $\leq 15\%$ in the GeV range. The inelastic cross sections computed
with FLUKA show a good agreement with data and with other dedicated
parameterisations, consistent in all the energy range studied within $\sim
25\%$ discrepancies. For interactions of nuclei heavier than Ne (Z=10) with
protons, for which experimental data in the GeV range is mainly absent, the
cross sections predicted from FLUKA differ from those predicted from dedicated
parameterisations by a $\sim 20\%$ factor, roughly constant in energy [8].

In turn, inclusive cross sections are that fraction of the inelastic cross
sections which go into the creation of a secondary nucleus and regulate the
rate of production of secondary CRs ($\tau_{p\rightarrow s}^{-1}\sim
v_{p}n_{ISM}\sigma_{{p\rightarrow s}}(E)$, where $v_{p}$ is the speed of the
projectile nucleus, p, and $\sigma_{{p\rightarrow s}}$ is the cross section of
production of the secondary particle, s, from the interaction of the
projectile particle with the target nucleus). These are poorly known because
of the experimental difficulties to perform the measurements and, therefore,
the uncertainties associated to inclusive cross-section data are of $\sim
20\%-30\%$ in the GeV range. On top of this, inclusive cross sections should
also include the decay of ghost nuclei (i.e. short-lived nuclei generated in
spallation reactions, which have a negligible lifetime compared to typical CR
propagation times, and decay into the secondary nucleus that we are
considering) in the current CR propagation codes. The inclusive cross sections
including this effect are referred to as cumulative inclusive cross sections.

Figure 1: Spallation cross sections of CRs interaction with ISM gas computed
with FLUKA compared to the most widely used parameterisations, for the
production of isotopes of B, Be and Li from 12C and 16O as projectiles.

Figure 1 shows the total (i.e. cross sections of interactions with a gas with
ISM composition) cumulative cross sections of production of the some isotopes
of B, Be and Li from C and O (the main producers of these secondary CRs).
Here, we compare the computed FLUKA cross sections with well-known dedicated
cross sections parameterisations, commonly used in CR codes: the GALPROP [13]
and DRAGON2 [1] parameterisations and the cross sections calculated with the
WNEW03 [14] and YIELDX [15] (TS98 in the legend) codes. As we see, the FLUKA
cross sections are very consistent with the predictions from the other
parameterisations, which is remarkable since our predictions completely rely
on theory-based interaction models and not in fits to cross sections data. We
highlight that both the normalization and energy dependence of these cross
sections is compatible with the most updated parameterisations in the whole
energy range. In addition, the position of the predicted resonances are also
in good agreement with those parameterised from data. We observe that such
good consistency is found in most channels of production of isotopes of B, Be
and Li from interactions of C, N, O, Ne, Mg and Si with gas, which constitute
the main primary CRs producing these secondary CRs. However, we observe that
the channels of production of F (produced mainly from Ne, Mg and Si) predict
much lower cross sections (by a factor of $\sim 2$), which seem to be not
compatible with the few existent data. Finally, we remark the importance of
having reliable cross sections predictions for those channels in which
experimental data is totally absent (meaning that the parameterisations are
merely extrapolations in these channels), as happens for the channels
production of heavy secondary cosmic rays, like 26Al, for example.

##  3 CR transport predictions using the FLUKA cross sections

To study the spectra of B, Be and Li derived with the FLUKA cross sections and
the propagation parameters inferred from their ratios, we implement these
cross sections in the DRAGON2 code. The same injection and propagation set-up
as presented our previous works [1, 11, 8] is used, which includes
reacceleration and neglects convection. Solar modulation is treated with the
Force-field approximation [16]. In this work, we focus on the results obtained
with the diffusion coefficient parameterised as in Equation 1. In this
equation, the free parameters entering the fit are the normalization, $D_{0}$,
the exponent of $\beta$ (speed of particle in units of the speed of light),
$\eta$, and the spectral index, $\delta$, while $R_{0}$, the rigidity at which
the diffusion coefficient is normalized, is set to $4$ GV. Then, we fix
$\Delta\delta=0.14$, $R_{b}=312$ GV and $s=0.04$ [11]. The use of the
diffusion coefficient of Eq. 1 requires using as injection a broken power-law
with a break set at $8$ GeV.

$D=D_{0}\beta^{\eta}\frac{\left(R/R_{0}\right)^{\delta}}{\left[1+\left(R/R_{b}\right)^{\Delta\delta/s}\right]^{s}}\,\,\,,$
(1)

The starting point of this analysis consist of finding the propagation
parameters that allow us to reproduce the boron-over-carbon (B/C) spectrum
reported by AMS-02. The injection spectra of nuclei up to iron are adjusted,
recursively, along with the propagation parameters, to reproduce AMS-02 data
[12]. This analysis yielded values of the propagation parameters entering in
the fit which are reported in Table 1. This table also reports those
parameters found for the same analysis for the DRAGON2 and GALPROP
parameterisations, for comparisons. As we can see, the parameters obtained
using the FLUKA cross sections are consistent within $1\,\sigma$ with those
obtained using dedicated parameterisations, except for the normalization of
the diffusion coefficient. This is something remarkable and achieved for the
first time with a set of cross sections derived from fundamental models of
nuclear interactions.

B/C best-fit parameters | FLUKA |  GALPROP |  DRAGON2
---|---|---|---
D0/H ($10^{28}$ cm2 s-1 kpc-1) | 0.82 $\pm$ 0.03 | 0.94 $\pm$ 0.04 | 0.97 $\pm$ 0.04
$v_{A}$ (km/s) | 23.3 $\pm$ 2.3 | 24.4 $\pm$ 3. | 22\. $\pm$ 3.6
$\eta$ | -0.67 $\pm$ 0.13 | -0.66 $\pm$ 0.12 | -0.78 $\pm$ 0.13
$\delta$ | 0.45 $\pm$ 0.01 | 0.45 $\pm$ 0.01 | 0.44 $\pm$ 0.01
Table 1: Propagation parameteres found in the MCMC analysis of the B/C flux
ratio. The halo size value, H, obtained from the 10Be flux ratios is $7.54$
kpc. The errors shown here correspond only to statistical uncertainties.

To evaluate the size of the magnetised halo we analyse the ratios of 10Be to
9Be and to the total flux of Be. Full details can be found in our reference
work [8]. This fit yields a value of the halo size of $\sim
7.5^{+1.13}_{-0.95}$ kpc, similar to estimations with other cross sections
($H\sim 6.76\pm 1$ kpc for the DRAGON2 cross sections and $H\sim 6.93\pm 0.98$
kpc for the GALPROP ones).

One of the most important tools for the study of the cross sections used to
evaluate the spectra of secondary CRs is the flux ratios among them (Be/B,
Li/B and Li/Be), which are shown in Fig. 3. At high energy, these ratios
mainly depend on the ratio of the cross sections of the CRs involved and the
injection spectra of primary CRs, which are adjusted to reproduce AMS-02 data.
In this figure, we show a band of statistical uncertainty associated to the
uncertainty related to the determination of the FLUKA cross sections. As can
be seen in the figure, while the Li/Be ratio reproduces within $1\,\sigma$
uncertainty AMS-02 data above $\sim 2$ GeV, the ratios involving B are
overproduced by a roughly constant $20-25\%$ above $3$ GeV. We remark that
this discrepancy is yet within the typical uncertainty on cross sections
measurements above the GeV. However, this discrepancy means that, as happens
with the parameterisations commonly used in CR studies, they are not able to
reproduce simultaneously B, Be and Li. A simultaneous solution to this
discrepancy and the discrepancy found in the value of $D_{0}$ with respect to
the parametrisations would be to renormalise (to lower) the cross sections of
B production by around a $20\%$. The result of applying this scaling to the B
production cross sections is shown in the figure as a dashed line and leads to
a good simultaneous reproduction of the AMS-02 spectra of B, Be and Li above
$\sim 3$ GeV.

Figure 2: Be/B, Li/B and Li/Be ratios obtained with the FLUKA computations,
together with the band of statistical uncertainties due to the spallation
cross sections calculation (since this is a Monte Carlo computation).

###  3.1 Combined analysis of B, Be and Li

The final test in the study of the spectra of secondary CRs evaluated with the
FLUKA cross sections is a combined analysis where we combine the ratios of B,
Be and Li to C and O (B/C, B/O, Be/C, Be/O, Li/C and Li/O) and its respective
secondary-to-secondary flux ratios (Be/B, Li/B, Li/Be), MCMC algorithm
presented in a precedent work [11]. This analysis includes scaling factors, as
nuisance parameters, that renormalise the cross sections of production of B,
Be and Li and, associated to this scaling factor, there is a penalty factor
that penalises large variations of the original cross sections. The result of
this analysis is a set of propagation parameters and scaling factors that are
able to simultaneously reproduce all the ratios of B, Be and Li within the
$1\,\sigma$ statistical uncertainties in the determination of the propagation
parameters. This is shown in Figure 3, where we report the ratios of Be/B,
Li/B and Li/Be, B/O, Be/O, Li/O and the ${}^{10}Be$/9Be ratio, included in the
fit. We highlight that the diffusion coefficient obtained from these combined
analyses allows us to reproduce, at the same time, the 3He/4He flux ratio,
shown in the lower-right panel, which is evidence that the diffusion
coefficient predicted by the light secondary CRs is compatible for the
different nuclei within uncertainties.

These analyses predict a spectral index of the diffusion coefficient $\sim
0.36$, which is significantly lower than the typical values found from B
ratios, although compatible with the value found for the combined analysis
with the GALPROP cross sections [11]. Furthermore, this value is consistent
with the basic predictions from wave-particle interactions, for which the
standard spectrum of plasma waves leads to a spectral index of the diffusion
coefficient which is $0.33\lesssim\delta\lesssim 0.5$. Finally, another
crucial point of this analysis is the cross-section scaling factors obtained.
These are $1.18$, $0.94$ and $0.93$, for B, Be and Li, respectively, with
$1\sigma$ statistical uncertainties of $\sim\pm 0.01$. However, we highlight
that there are other systematic uncertainties related to the determination of
these scale factors, like those related to the gas distribution used, and
could make the total systematic uncertainties in the determination of these
scale factors larger than $5\%$ [11].

Figure 3: Top-left panel: Secondary-over-secondary flux ratios of B, Be and Li
computed using the propagation parameters and cross sections scaling factors
found in the combined analysis, compared to the AMS-02 experimental data.
Lower-left panel: B/O, Be/O and Li/O flux ratios predicted with the parameters
determined by the MCMC combined analysis for the DRAGON2 cross sections. Top-
right panel: 10Be/9Be flux ratio compared to available experimental data.
Lower-right panel: Predicted 3He/4He flux ratio compared to Pamela and AMS-02
data. Residuals, calculated as model-data/model are also shown. This ratio is
calculated with the diffusion coefficient obtained from the combined fit.

##  4 Evidence of a low-energy break in electron injection spectrum from
gamma-ray data

Finally, we employed the Local Interstellar spectrum (LIS) predicted from the
combined analysis described above to study the diffuse gamma-ray emission,
which is computed using the Gammasky code In particular, we tested the local
HI gamma-ray emissivity spectrum using Fermi-LAT gamma-ray data [17]. The
dominant contribution to this observable at high energies is the $\gamma$-ray
emission produced from the decay of unstable particles formed via nuclear
reactions (hadronic emission), such as pions. As it is observed in Figure 4,
the hadronic emission allows us to reproduce the Fermi-LAT emissivity, within
$1\,\sigma$ statistical uncertainties, assuming an ISM composition with
relative abundance of H : He : C : N : O : Ne : Mg : Si = $1:0.096:4.65\times
10^{-4}:8.3\times 10^{-5}:8.3\times 10^{-4}:1.3\times 10^{-4}:3.9\times
10^{-5}:3.69\times 10^{-5}$.

At low energy ($E_{\gamma}<1$ GeV) bremsstrahlung emission from leptons
becomes dominant. In order to evaluate the bremsstrahlung emission, we
reproduce the AMS-02 e${}^{+}+$e- data. We tested two different models that
allow us to reproduce electrons: (i) a broken power-law with a break at $65$
GeV and a modified Force-field approximation, able to account for the polarity
and charge sign of the different particles [18]. (ii) a doubly-broken power-
law with a break at $8$ GeV (same break position as it is commonly employed
for nuclei) and a high energy break at $65$ GeV, using the conventional Force-
field approximation [16]. In this study, we found that the first option
overproduces very significantly the local $gamma$-ray emissivity at low
energies, while the second approach allows us to reproduce the emissivity
within the uncertainties associated to the Fisk potential, assuming to be of
$\phi\pm 0.1$, motivated by neutron monitors [19].

This break was observed from synchrotron radio emission [20] but never proved
to be necessary to reproduce gamma-ray data, as we determine here. To further
test the hypothesis of a break at the few GeV range in the electron spectrum
we studied the diffuse gamma-ray emission at MeV energies around the center of
the Galaxy ($|b|<47.5^{\circ},|l|<47.5^{\circ}$), where Inverse-Compton (IC)
emission produced from leptons interacting with the interstellar radiation
fields is totally dominant, as it can be seen from Figure 5 (left), where we
show the predicted gamma-ray emission around the MeV for the model including a
break at $8$ GeV in the electron spectrum. Hadronic emission in this plot
corresponds to the gamma-ray emission from exited states of nuclei, computed
using the FLUKA cross sections. As it can be seen, our predictions reproduce
well the experimental data from COMPTEL telescope [21]. Figure 5 right shows a
comparison between the predictions with and without a break in the electron
spectrum, clearly supporting the break hypothesis, as the local emissivity
does too.

Figure 4: Local HI gamma-ray emissivity spectrum for the propagation
parameters derived from the combined analysis. An uncertainty band related to
solar modulation uncertainty is also shown.

##  5 Conclusions

In this work, we presented the FLUKA cross sections for Galactic CR
propagation studies as an alternative and an improvement with respect to
dedicated cross sections parameterisations. We have tested the spectra of B,
Be and Li and evaluated the main propagation parameters inferred from AMS-02
data, demonstrating that every CR observable can be reproduced with accuracy
and at the level of precision of current CR parameterisations.

We showed that the FLUKA cross sections allow us to reproduce simultaneously
the light secondary CRs B, Be, Li and 3He, when introducing nuisance scaling
factors to renormalise the overall cross sections of production of B, Be and
Li. The scaling factors obtained are of $\sim 1.18$, $\sim 0.94$ and $\sim
0.93$ for B, Be and Li, respectively, which are below the typical average
experimental uncertainties of cross sections measurements for the best known
channels.

Finally, from the study of the diffuse gamma-ray emission we have determined
the need of including a break in the injection spectrum of electrons at a few
GeV, inferred from the local HI gamma-ray emissivity and the diffuse IC
emission at MeV energies.

Figure 5: Predicted diffuse gamma-ray emission in the Galactic region
$|b|<47.5^{\circ},|l|<47.5^{\circ}$. The left panel shows the different
contributions to the diffuse gamma-ray emission at those energies, while the
right panel shows the predicted diffuse gamma-ray emission for the hypothesis
of a low-energy break in the electron injection spectrum compared to the
hypothesis without break.

## Acknowledgments

P. De la Torre is supported by the Swedish National Space Agency under
contract 117/19. We acknowledge the FLUKA collaboration for providing and
supporting the code.

## References

## References

  * [1] Pedro De La Torre Luque et al.  Implications of current nuclear cross sections on secondary cosmic rays with the upcoming DRAGON2 code.  JCAP, 03:099, 2021.
  * [2] N. Weinrich and others.  Combined analysis of AMS-02 (Li,Be,B)/C, N/O, 3He, and 4He data.  The Astrophysical Journal, 639:A131, July 2020.
  * [3] M. Korsmeier and A. Cuoco.  Implications of Lithium to Oxygen AMS-02 spectra on our understanding of cosmic-ray diffusion.  Phys. Rev. D, 103(10):103016, 2021.
  * [4] Giuseppe Battistoni et al.  Overview of the fluka code.  Annals of Nuclear Energy, 82:10–18, 2015.  SNA + MC Conference, 2013.
  * [5] M. N. Mazziotta and otherz.  Cosmic-ray interactions with the Sun using the FLUKA code.  Phys. Rev. D, 101(8):083011, 2020.
  * [6] M. N. Mazziotta et al.  Production of secondary particles and nuclei in cosmic rays collisions with the interstellar gas using the FLUKA code.  AstroPart. Phys., 81:21–38, 2016.
  * [7] M. Ackermann et al.  Measurement of the high-energy gamma-ray emission from the Moon with the Fermi Large Area Telescope.  Phys. Rev. D, 93(8):082001, 2016.
  * [8] Pedro de la Torre Luque et al.  FLUKA cross sections for cosmic-ray interactions with the DRAGON2 code.  2 2022.
  * [9] Carmelo Evoli et al.  Cosmic-ray propagation with $\small{DRAGON2}$: I. numerical solver and astrophysical ingredients.  JCAP, 02:015, 2017.
  * [10] Carmelo Evoli et al.  Cosmic-ray propagation with DRAGON2: II. Nuclear interactions with the interstellar gas.  JCAP, 07:006, 2018.
  * [11] P. De La Torre Luque et al.  Markov chain monte carlo analyses of the flux ratios of b, be and li with the DRAGON2 code.  JCAP, 2021(07):010, jul 2021.
  * [12] M. Aguilar et al.  Observation of new properties of secondary cosmic rays lithium, beryllium, and boron by the alpha magnetic spectrometer on the international space station.  Phys. Rev. Lett., 120:021101, Jan 2018.
  * [13] Igor V. Moskalenko and S. G. Mashnik.  Evaluation of Production Cross Sections of Li, Be, B in CR.  In ICRC, volume 4 of ICRC, page 1969, July 2003.
  * [14] WR Webber, A Soutoul, JC Kish, and JM Rockstroh.  Updated formula for calculating partial cross sections for nuclear reactions of nuclei with z$\leq 28$ and e$>150$ mev nucleon–1 in hydrogen targets.  ApJ Supplement Series, 144(1):153, 2003.
  * [15] R Silberberg, CH Tsao, and AF Barghouty.  Updated partial cross sections of proton-nucleus reactions.  ApJ, 501(2):911, 1998.
  * [16] L. J. Gleeson and I. H. Urch.  A Study of the Force-Field Equation for the Propagation of Galactic Cosmic Rays.  Astrophysics and Space Science, 25(2):387–404, December 1973.
  * [17] Jean-Marc Casandjian.  Local HI emissivity measured with Fermi-LAT and implications for cosmic-ray spectra.  Astrophys. J., 806(2):240, 2015.
  * [18] Ilias Cholis, Dan Hooper, and Tim Linden.  A predictive analytic model for the solar modulation of cosmic rays.  Phys. Rev. D, 93:043016, Feb 2016.
  * [19] A. Ghelfi et al.  Neutron monitors and muon detectors for solar modulation studies: 2. $\phi$ time series.  Adv. Space Res., 60(4):833–847, 2017.
  * [20] Giuseppe Di Bernardo et al.  Cosmic ray electrons, positrons and the synchrotron emission of the galaxy: consistent analysis and implications.  JCAP, 2013(03):036–036, mar 2013.
  * [21] Siegert, T., Berteaud, J., Calore, F., Serpico, P. D., and Weinberger, C.  Diffuse galactic emission spectrum between 0.5 and 8.0 mev.  A&A, 660:A130, 2022.