text:
\section{Introduction}
\label{intro}

Aside from its intrinsic interest, a detailed description of the inner parts of the Earth is essential for a proper understanding of basic geological phenomena such as volcanology, earthquakes, plate tectonics, and mountain building
\cite{tarbuck2009earth,fowler:2005}. According to current knowledge, the Earth's internal structure is stratified and consist of successive layers with different chemical, geological, and physical properties. Most of the information about the layers and the boundaries between them has been acquired by examining how seismic waves created by natural earthquakes are refracted and reflected as they propagate through the Earth \cite{Lay:1995,Aki:2002}. Valuable information has also been obtained from measurements of magnetic and gravitational fields, observations of the planet's moment of inertia and precessional motion, and physical, chemical, and mineralogical analyzes of meteorites and xenoliths. Based on chemical composition, three main layers have been identified within our planet: the crust, the mantle, and the core. The mantle and the core are further subdivided into two regions each.

The mantle surrounds the core and extends from beneath the thin crust to a depth of 2,900 km. It consists mainly of dense silicate rocks rich in iron and magnesium. Seismic wave velocities in the mantle show three discontinuities at depths of 410, 660 and 2,700 km \cite{Helffrich:2001}. The first two correspond to the edges of the transition zone between the upper and lower mantle, and are best explained by mineral phase transformations without compositional changes. The third discontinuity is typically a 2.5-3.0\% increase for both S- and P-waves observed at the top of the D$^{\prime\prime}$ layer, a transition shell about 200-300 km thick, at the base of the mantle, which presents a variety of seismic anomalies and is presumably the source of the large mantle plumes \cite{Loper:1995,Laywg:1998}. This layer is not yet fully understood, many of its characteristic may be attributed to the discovered MgSi${\rm O}_3$ postperovskite phase \cite{Murakami:2004}, but compositional differences may also play an important role in addition to the phase transformation \cite{Hiroshe:2021}. Regarding the core, it is well established that it is composed primarily of an iron-nickel alloy, with Ni/Fe\,$\sim 0.06$. The inner core is solid, while the lower density and absence of S-wave propagation are indications of a liquid outer core. The density deficit in the outer core cannot be simply explained by a difference in the state, but it requires about 5-10 wt\% (weight percent) of lower atomic weight elements to reduce its density and melting point \cite{Birch:1964, Poirier:1994}.

Good estimates of the abundance and distribution of ``light''  elements in the core are essential to understanding the formation and evolution of the Earth, as well as how the core and mantle interact in the region around the core-mantle boundary (CMB) \cite{Litasov:2016,Hiroshe:2021}. The various processes that occur in this interface are highly influenced by heterogeneous structures at or near it. The  great disparity in density prevents direct convection movement through the CMB, and nearby regions control the transfer of heat and material \cite{Wookey:2008}. These processes strongly affect the convection in the mantle, responsible for the plate motion and continental drift, and the more vigorous convective flow in the outer core that is believed to be at the origin of the Earth's magnetic field (i.e., the functioning of the geodynamo) \cite{Lister:1995,Roberts:2013}. Due to vigorous convection, the liquid part of the core is usually assumed to have a uniform composition. However, seismological evidence indicates the possible existence of a $\sim 300$ km thick region on top of the outer core (E$^\prime$ layer), which shows anomalous low seismic velocities. This region is likely to be less dense than the rest of the outer core, but simply increasing the concentration of light elements also produces higher speeds, in contradiction with observations. None of the mechanisms that have been considered are without complications \cite{Brodholt:2017} and further observations are needed before a satisfactory explanation of the origin and nature of the E' layer can be formulated.

The study of seismic wave propagation and normal mode oscillation is undoubtedly the most effective and reliable method to search the Earth's deep structures and process. Nevertheless, the impressive progress done has not been accompanied with a concomitant improvement in the precision of the density estimations \cite{Bolt:1995}. They are done through an underdetermined inversion problem, performed in two steps: the spatial distribution of seismic wave velocities is first inferred from seismological data and then the density distribution is inferred from the seismic velocities using some empirical relation \cite{Geller:2003}.  Such procedure allows the average density along a path to be estimated with an uncertainties of about 5\% \cite{Bolt:1995} for the mantle and presumably larger for the core. For example, the density jump at the inner-core boundary, which play an important role in the maintenance of the geodynamo, has been inferred to be of 0.82 with an error of more than 20\% \cite{Master:2003}. While the density distribution can be obtained from seismological remote sensing, the compositional structure of the Earth  \cite{McDonough:1995,allegre1995chemical} has been much more difficult to determine. Thus, the compositions of the lower mantle and the core remain quite uncertain, despite significant advances in recent years.  Since in situ sampling is impossible, estimations are done by comparing density and sound velocity data from seismological observations with those from laboratory experiments and theoretical calculations \cite{Zhang:2016}. Due to technical limitations, it has been difficult to perform reliable experiments for molten samples under high pressure and temperature conditions and the available information is still insufficient to infer the core composition. The most likely light elements in the core are oxygen, silicon, carbon, nitrogen, sulfur, and hydrogen \cite{Li:2021,Hiroshe:2021}, but there is still no consensus on the nature and proportion of the components.

In addition to the metallic core, significant amounts of hydrogen can be incorporated within the mantle, in silicate minerals, melts, and hydrous fluids. This is done in a variety of chemical species (OH, H$_2$, $\cdots$) generically referred to as ``water" \cite{Hirschmann:2006,Peslier:2017,Othani:2020}. The abundance and distribution of water has influenced not only the evolution, dynamics state and thermal structure of the deep interior, but also the evolution of the crust and hydrosphere \cite{Bodnar:2013}.
Even small amount of water can affect properties like melting  temperature, rheological properties, electrical conductivity, and seismic velocities of the mantle. How water is transported into Earth's deep interior and how it is distributed are today open questions. Almost no constraints exists on the water content of the lower mantle. The existing laboratory data cannot be used to infer the water content from geophysical observations at these high pressure and temperature conditions. Low-velocity regions have been observed near the top of the lower mantle and directly above the CMB \cite{Peslier:2017}. If the low-velocity is interpreted as caused by partial melting, then some water is likely be present in these regions \cite{Garnero:2016}. Water could instead be stored in the lowest parts of the mantle in Al-postperovskite \cite{Townsend:2016}. There is a wide consensus that the mantle transition zone (MTZ), at 410-660 km deep, is a potential reservoir of water because its main mineral constituents can store up to $\sim 3\%$ wt water \cite{Bodnar:2013}. However, the amount of water contained there is poorly constrained. Some water-rich inclusions recently found in diamonds \cite{Pearson:2014,tschauner2018ice} suggest a wet MTZ, but it is not clear if they  are representative of the typical water content of the deep mantle or reflect local conditions. Geophysical methods (electrical conductivity and seismic observations) provides constraints on the water distribution in a global scale \cite{karato2011water}. The inferences with these approaches are not direct and, despite the effort involved, a wide range of values have been reported in the literature \cite{karato:2020,Fei:2017}. As noticed in Ref. \cite{karato:2020}, geophysical estimates have large uncertainties and the water content in the MTZ can be heterogeneous, having significant lateral variability as revealed by a recent novel approach \cite{munch:2020}.

From the comments in the previous paragraph it is apparent how useful it would be to have additional experimental techniques, unaffected by the same uncertainties, which could provide complementary and independent information about the deep interior of the Earth. A promising candidate is neutrino tomography \cite{Winter:2006vg}.  The basic idea is that neutrino propagation within the Earth is affected by their interactions with the particles present in the terrestrial matter. The cross section for neutrino interactions increase with energy and, in the case of absorption tomography, the density profile can be reconstructed from the attenuation of the flux of very high energy ($\gtrsim 10$ TeV) neutrinos passing through the Earth \cite{Volkova:1974xa, Wilson:1983an, Ralston:1999fz, Jain:1999kp, GonzalezGarcia:2007gg, Reynoso:2004dt, Romero:2011zzb,Donini:2018tsg}. Another option is oscillation tomography, which takes advantage of the matter effect on flavor oscillations \cite{Wolfenstein:1977ue, Barger:1980tf, Mikheev:1986gs} of lower energy (MeV to GeV) neutrinos. In a medium, the transition probabilities between active neutrinos depend on the number density of electrons $n_e$ along the neutrino trajectory, which is proportional to the product of the matter density $\rho$ times the average ratio $Z/A$ of the atomic number $Z$ and the mass number $A$ \cite{Nicolaidis:1987fe, Nicolaidis:1990jm, Winter:2006vg, Borriello:2009ad, Winter:2015zwx, Rott:2015kwa, VanElewyck:2017dga, Bourret:2017tkw, DOlivo:2020ssf, Kelly:2021jfs, Denton:2021rgt, Upadhyay:2021kzf}.

In this paper, we re-examine the feasibility of studying the internal structure of the Earth using atmospheric neutrino oscillation tomography.  We analyze the ability of a detector such as ORCA to resolve deviations both of density and composition with respect to a standard Earth modeled in terms of 55 concentric shells integrated in five main layers, corresponding to the inner and outer core, the lower and upper mantle, and the crust. Constant densities are assigned to the shells from the mean value of the PREM densities within each shell.  In our scheme, unlike other work on the subject, the densities of the main layers can be modified in a manner consistent with the well-measured total mass and moment of inertia of the Earth. In contrast, $Z/A$ is a function of the chemical and isotopic composition of the medium and is not subject to either of these constraints. Therefore, in principle, one could constrain the allowed values of the density and composition of the Earth's deeper regions  by studying the effects that changes in these quantities have on the events produced by atmospheric neutrinos after traversing the Earth, holding the total mass and moment of inertia of the planet fixed.  We focus on the possible application to obtain information on the composition and density of the lower mantle regions above the CMB, in particular, the D$^{\prime \prime}$ region, paying special attention to the content of light elements, more specifically hydrogen.  With the exception of hydrogen ($(Z/A)_H = 1$), all other light elements have an almost equal number of protons and neutrons and hence $Z/A\cong 0.5$  for them. Thus, the presence of a significant amount of hydrogen (water) would produce appreciable $Z/A$ changes compared to those in dry regions. In this sense, information from neutrino oscillation tomography could give valuable information on the water content in the lower mantle and core. We also allow for some variation in the location of the boundary between the lower and upper mantle as an effective way to account for the transition region between these two layers.

The paper is organized as follows. In Sect. \ref{sec:structure} we present a model of the Earth's structure with 55  shells integrated within 5 main layers. In Sect. \ref{sec:osc} we briefly review the formalism of matter neutrino oscillations and describe the algorithm to calculate the transition  probabilities. In Sect. \ref{sec:simu} we determine the number of $\mu$-like neutrino events in a detector
such as ORCA and the effects that changes in the composition  and density of the outer core and lower mantle have on this observable. The results and final comments are presented in Sect.\ref{sec:results}, where we carry out a Monte Carlo simulation of the number of  $\mu$-like events and apply it to test different composition models of the outer core and mantle, changing also the density.

\section{Model of the Earth's structure}
\label{sec:structure}
The observed lateral variations in the Earth's properties are much less pronounced than the vertical variations. Therefore, the internal structure of the Earth can be well approximated by one-dimensional spherically symmetric models of seismic velocities, attenuation, and density as a function of depth \cite{dziewonski1981preliminary,Kennett:1995}.The most widely used of these models for seismic tomography has been the Preliminary Reference Earth Model (PREM) \cite{dziewonski1981preliminary}. This model represents the mean properties of the Earth as a function of the radial distance $r$ and was designed to fit different data sets and some basic data of the planet (radius, mass, and moment of inertia). In this study,  we use a spherical model consisting of 55 concentric shells, each with a constant density equal to the average value of the PREM densities in the shell. The set of shells is divided into five large layers demarcated by concentric spheres of different radii: inner core ($IC$), outer core ($OC$), lower mantle ($M_1$), upper mantle ($M_2$), and crust ($C$). In Table \ref{tab:compo} we give the values for the standard composition and the corresponding radii.
\begin{figure}[h!]
\centering
\includegraphics[width=0.8\textwidth]{densidad.pdf}
\caption{PREM density profile of the Earth with 55 shells.}
\label{fig:prem55}
\end{figure}

\begin{table}[h]
\centering
\begin{tabular}{@{}|c|c|c|c|@{}}
\toprule
\textbf{Layer} & \textbf{$n^0$ of Shells} & \textbf{$R_{inf}$ - $R_{sup}$ {[}km{]}} & \textbf{Z/A} \\ \midrule
Inner Core   & 7              & 0 - 1221.5         & 0.4691 \\
Outer Core   & 13             & 1221.5 - 3480      & 0.4691 \\
Lower Mantle & $N_{M_1}$-$21$ & $3480$ - $R_{M_1}$ & 0.4954 \\
Upper Mantle & $49$-$N_{M_1}$ & $R_{M_1}$ - $6346$ & 0.4954 \\
Crust        & 6346-6371      & 6346 - 6371        & 0.4956 \\ \bottomrule
\end{tabular}
\caption{Compositions of the main Earth layers.}
\label{tab:compo}
\end{table}

The primary information on the Earth's density as a function of $r$ comes from the total mass of the Earth
$
M_{_{\hsn\oplus}} = 5.9724 \times 10^{27} {\rm g}
$
and its mean moment of inertia about the polar axis
$
I_{_{\hsn\oplus}} = 8.025 \times 10^{44} {\rm g} \hs{0.03 cm}{\rm cm}^2
$
\cite{Williams:1994}.
From these two quantities and $R_{_{\hsn\oplus}} = 6371$ km for the Earth's radius, one gets  $I_{_{\hsn\oplus}} \simeq 0.33 M_{_{\hsn\oplus}} R^2_{_{\hsn\oplus}}$ that is noticeably smaller than the moment of inertia of a homogeneous sphere of the same radius ($0.40 M_{_{\hsn\oplus}} R^2_{_{\hsn\oplus}}$). This corroborates that there must be a concentration of mass towards the center of the planet or, in other words, that the inner regions are denser than average. We require that our model satisfy both constraints:
\begin{eqnarray}
\label{masa_inertia_st}
M_{_{\hsn\oplus}} &=& \frac43 \pi \sum_{i=1}^{55} \rho_i (R_{i+1}^3-R_i^3) = M_{\chic{IC}}+M_{\chic{OC}}+M_{M_1}+M_{M_2}+M_C
\nonumber \\
I_{_{\hsn\oplus}}&=& \frac8{15} \pi \sum_{i=1}^{55} \rho_i (R_{i+1}^5-R_i^5) = I_{\chic{IC}}+I_{\chic{OC}}+I_{M_1}+I_{M_2}+I_C
\end{eqnarray}
where  $M_i$ and $I_i$, $i = IC, OC, M_1, M_2, C$, are the masses and moments of inertia of the major layer
specified above, which are given by the sums of the contributions of the shells contained in each of these divisions,
as indicated in Table 1.

To modify the densities of the outer core and the lower and upper mantle we multiply the densities of all
shells within each of these layers by the respective rescaling factor, $f_{\chic{OC}}, f_{M_1}$, and $f_{M_2}$.
This is done in such a way that neither $M_{_{\hsn\oplus}}$ nor $I_{_{\hsn\oplus}}$ change.
Then,
\begin{eqnarray}
\label{masa_inertia_md}
M_{_{\hsn\oplus}} &=& M_{\chic{IC}}+f_{\chic{OC}} \; M_{\chic{OC}}+ f_{M_1} \; M_{M_1}+ f_{M_2} \; M_{M_2}+M_{\chic{C}}
\nonumber\\
I_{_{\hsn\oplus}}&=&I_{\chic{IC}}+ f_{\chic{OC}} \; I_{\chic{OC}}+f_{M_1} \; I_{M_1}+f_{M_2} \; I_{M_2}+I_{\chic{C}}
\end{eqnarray}
Equating Eqs. \eqref{masa_inertia_st} and \eqref{masa_inertia_md} we obtain the following homogeneous system of
linear equations:
\begin{eqnarray}
\delta_{\chic{OC}} \; M_{\chic{OC}}+\delta_{M_1} \; M_{M_1}+\delta_{M_2} \; M_{M_2} &=& 0\,,
\nonumber \\
\delta_{\chic{OC}} \; I_{\chic{OC}}+\delta_{M_1} \; I_{M_1}+\delta_{M_2} \; I_{M_2} &=& 0\,,
\end{eqnarray}
where $\delta_{\chic{OC}} = f_{\chic{OC}} -1$, $\delta_{M_1} = f_{M_1} -1$, and $\delta_{M_2} = f_{M_2} -1$ are the
relative changes of the densities in the outer core, lower mantle, and upper mantle, respectively. Solving
this system, we can express $\delta_{M_1}$ and $\delta_{M_2}$ as functions of $\delta_{\chic{OC}}$:
\begin{eqnarray}
\label{reldeltaocdeltaM12}
\delta_{M_1}=-\delta_{\chic{OC}} \;\; \frac{\Delta_{M_1}}{\Delta}\,, \qquad \delta_{M_2}=-\delta_{\chic{OC}} \;\; \frac{\Delta_{M_2}}{\Delta}\,,
\end{eqnarray}
where
\begin{eqnarray}
\Delta &=& M_{M_1} \; I_{M_2} - I_{M_1} \; M_{M_2},\nonumber\\
\Delta_{M_1} &=& M_{\chic{OC}} \; I_{M_2} - I_{\chic{OC}} \; M_{M_2},\\
\Delta_{M_2} &=& I_{\chic{OC}} \; M_{M_1} - M_{\chic{OC}}\;I_{M_1}.\nonumber
\end{eqnarray}

The value of the radius $R_{M_1}$ set the position of the boundary between the regions $M_1$ and $M_2$ and varying it we can change the number of shells within each of these layers. This, in turn, modifies the
values of $M_{M_{1,2}}$ and $I_{I_{1,2}}$ and makes the quantities $\delta_{1,2}$  dependent on $R_{M_1}$.
Fig. \ref{fig:dens} shows the relative changes in the densities of layers $M_1$ and $M_2$ as a function of the relative change in the density of the outer core, for three different positions of the boundary between $M_1$ and $M_2$.
\begin{figure}[h!]
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{densM1.pdf}
\caption{}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{densM2.pdf}
\caption{}
\end{subfigure}
\vskip\baselineskip
\hfill
\caption{Relative changes in the densities of the layers (a) $M_1$ and (b) $M_2$ as a function of the relative change in the density of the outer core, for $R_{M_1}$equal to
4800 km, 5600 km, and 5871 km.}
\label{fig:dens}
\end{figure}
\begin{table}
\centering
\begin{tabular}{|c|c|c|}
\hline
\textbf{Parameter} & \textbf{Normal Ordering} & \textbf{Inverted Ordering} \\ \hline
$\Delta m^{2}_{21}$ {[}$\text{eV}^{2}${]} & $7.42\times 10^{-5}$ & $7.42\times 10^{-5}$ \\
$\Delta m^{2}_{31}$ {[}$\text{eV}^{2}${]} & $2.533\times 10^{-3}$ & $-2.437\times 10^{-3}$ \\
$\sin^{2}{\theta_{12}}$ & 0.309 & 0.308 \\
$\sin^{2}{\theta_{13}}$ & 0.0223 & 0.0232 \\
$\sin^{2}{\theta_{23}}$ & 0.561 & 0.564 \\
$\delta/\pi$ & $1.19$ & $1.54$ \\ \hline
\end{tabular}
\caption{Three-neutrino oscillation parameters obtained by averaging the best-fit values of three recent global fits of the current neutrino oscillation data \cite{Capozzi:2020,Esteban:2020cvm,deSalas:2020pgw}.}
\label{tab:constants}
\end{table}
\section{Atmospheric neutrino oscillations}
\label{sec:osc}

Neutrino oscillations are a well-verified and widely studied phenomenon that has proven beyond any doubt that neutrinos are mixed massive particles. The current experimental and observational data set can be interpreted in terms of the minimal extension of the Standard Model, where the known flavor states $|{\nu}_{\alpha}\rangle (\alpha = e, \mu, \tau)$ are linear combinations of the states $|{\nu}_i \rangle$ with masses $m_i\,(i= 1,2,3)$:
$
\label{flav&massrel}
|{\nu}_{\alpha}\rangle= \sum_i U^*_{\alpha i} |{\nu}_i \rangle.
$
The coefficients $U_{\alpha i}$ appear in the leptonic charged current and are elements of a unitary matrix $U$. For for Dirac neutrinos, this matrix can be expressed as
$
\label{eq:Ufactorization}
U= {\cal O}_{23}\rm\Gamma{\cal O}_{13}{\cal O}_{12}\rm\Gamma^*\!,
$
\noindent
with
\bea
{\cal O}_{12} =
      \left(
           \begin{array}{ccc}
            c_{12} & s_{12} & 0 \\
           -s_{12} & c_{12} & 0 \\
               0   &    0   & 1
           \end{array}
        \right),
&
\ {\cal O}_{13} =
      \left(
           \begin{array}{ccc}
            c_{13} &  0  & s_{13} \\
	       0   &  1  &   0    \\
           -s_{13} &  0  & c_{13}
           \end{array}
        \right),
\nonumber \\
{\cal O}_{23} =
      \left(
           \begin{array}{ccc}
            1   &      0   &    0   \\
            0   &   c_{23} & s_{23} \\
            0   &  -s_{23} & c_{23}           \end{array}
        \right) ,
&
\!\!\! \Gamma =
      \left(
           \begin{array}{ccc}
            1   &      0   &    0   \\
            0   &      1   &    0   \\
            0   &      0   & {\rm e}^{i\delta}
           \end{array}
        \right)  ,
\eea

\noindent
where $c_{ij}=\cos\theta_{ij}$ and $s_{ij}=\sin\theta_{ij}$. A value of $\delta$ different form 0 or $\pi$ implies CP-violation
in the leptonic sector of the theory \!\footnote{For Majorana neutrinos there are two additional physical phases, but they are not relevant in neutrino oscillations and are therefore omitted in the analysis of the phenomenon \cite{Bilenky:1980cx}.}\!.

In addition to the mixing angles $\theta_{ij}$ and the CP-violating phase $\delta$, the oscillations between the three active neutrinos are parametrized by two squared-mass differences: $\Delta m^2_{21} \equiv m^2_2 - m^2_1$ and  $\Delta m^2_{32} \equiv m^2_3 - m^2_1$. Five of the parameters ($\theta_{12}, \theta_{13}, \theta_{23}, \Delta m^2_{21}$, and $|\Delta m^2_{31}|$), have been determined with remarkable precision ($\thicksim1-5 \%$) by global fits of the data from solar, atmospheric, reactor and long baseline experiments. The issues still pending are: the sign of $\Delta m^2_{31}$, the $\theta_{23}$ octant, and the determination orfthe phase $\delta$. The sign of $\Delta m^2_{31}$ characterize the normal ordering (NO), with $m_3 > m_{1,2}$, and the inverted ordering (IO), with $m_3 < m_{1,2}$. The 3$\nu$ oscillation parameters  shown in Table 2 are the mean of the best-fit values for the allowed ranges at $1\sigma$ of the global analyses performed by three groups \cite{Capozzi:2020,Esteban:2020cvm,deSalas:2020pgw}.

Neutrinos are produced and detected as flavor eigenstates. Consider a neutrino $\nu_\alpha$ produced at time $t_0$
that propagates in vacuum. Due to slight mass differences, the phases of the mass eigenstate components of the original flavor state change at different rates and, due to this, the flavor content of the neutrino beam oscillates along the trajectory. When neutrinos propagate in a medium, the coherent forward scattering of neutrinos with electrons is different for  $\nu_e$ and $\nu_{\mu,\tau}$, resulting in different refraction indexes for the electron neutrino and the other flavors. As a consequence, neutrino oscillations can be significantly modified  in matter compared to oscillations in vacuum and new resonance enhancement effects appear. These effects are sensitive to the density and composition of the medium and we will take advantage of this in order to examine the inner parts of our planet by means of the oscillations of atmospheric neutrinos in the Earth.
Atmospheric neutrinos have played a very important role in the study and characterization of the phenomena of neutrino oscillations. They are generated around the Earth as decay products in hadronic showers that result from collisions of cosmic rays with nuclei in the upper atmosphere. This provide a continuous source of neutrinos spanning a very wide range of energies and travelled distances before detection. On this work, we concentrate on those with energies in the range of 1-10 GeV and different nadir angles.

Let a neutrino $\nu_\alpha$ that enters the solid terrestrial matter at time $t_0$. At any time $t > t_0$ the state of the system $|\psi(t)\rangle$ can be expressed as $|\psi(t)\rangle = \hat {\mathcal U}(t, t_0)|\psi(t_0)\rangle$, where $|\psi(t_0)\rangle = |\nu_{\alpha}\rangle$ and $\hat{\mathcal U}(t, t_0)$ is the evolution operator. The probability of having a neutrino of flavor $\beta$ inside the Earth, at a distance $\ell \simeq t - t_0 \hs{0.02 cm} (\hslash = c = 1)$ from the entry point, is
\begin{equation}
\label{oscprob}
P_{{\nu_\alpha} \to {\nu_\beta}}(\ell) = |\,{\mathcal U}_{\beta \alpha}(\ell)|^2\,,
\end{equation}
where the probability amplitude
$
\label{matrixelem}
{\mathcal U}_{\beta \alpha}(\ell)= \langle \nu_{\beta}|\,\hat {\mathcal U}(\ell)|\nu_{\alpha}\rangle
$
is an element of the $3 \times 3$ unitary matrix $\mathcal{U}(\ell)$ representing the evolution operator in the flavor basis.
\begin{figure}[h!]
\centering
\includegraphics[width=0.45\textwidth]{Earth_path.pdf}
\hfill
\caption{Neutrino path through the Earth.}
\label{fig:earth}
\end{figure}
Rather than solving for $\mathcal{U}(\ell)$ directly, it proves to be more convenient to determine the evolution operator in the basis of the mass eigenstates and then transform it to the flavor basis using the relation
\begin{equation}
\mathcal{U}(\ell) = \tilde U \tilde{\mathcal{U}}(\ell) \tilde{U}^{\dagger} ,
\end{equation}
where $\tilde U =  {\cal O}_{23}\rm\Gamma{\cal O}_{13}{\cal O}_{12}$.

The matrix $\tilde{\mathcal{U}}(\ell)$ obeys the equation
\!\footnote{
The operator $\tilde{\mathcal{U}}(\ell)$ evolves the wave function $\tilde{\Phi}(\ell) = {\rm\Gamma^*}\Phi(\ell)$, with $\Phi^{\mathtt T}(\ell)  =
\left(\phi_1(\ell), \phi_2,(\ell), \phi_3(\ell)\right)$, where
 $\phi_i(\ell), i =1,2,3$ are the amplitudes of the mass eigenstates.}
\begin{equation}
\label{Scheq}
i \frac{d}{d \ell}\tilde{\mathcal{U}}(\ell) = \tilde H (\ell) \tilde{\mathcal{U}}(\ell), \qquad \tilde{\mathcal{U}}(0) = I,
\end{equation}
with
\begin{equation}  \
\label{eq:Hmatter}
\tilde H (\ell) =  H_{\rm 0}  + {\mathtt V}(\ell) {\mathcal O}^{\mathtt T} Y {\mathcal O}\,.
\end{equation}
Here, $H_0 = {\rm diag}\{0, \Delta_{21}, \Delta_{31}\}$  and $Y = {\rm diag}(1,0,0)$ are diagonal matrices
and ${\mathcal O}^{\mathtt T}$ is the transpose of the orthogonal matrix ${\mathcal O} = {\cal O}_{13}{\cal O}_{12}$.
Note that the real and symmetric matrix ${\mathcal O}^{\mathtt T} Y {\mathcal O} $ does not depend on $\theta_{23}$
nor $\delta$. Explicitly,
\begin{equation}
\label{eq:matrixW}
 {\mathcal O}^{\mathtt T} Y {\mathcal O}= \left(\!
\begin{array}{ccc}
c_{\rm{\chic{13}}}^2c_{\rm{\chic{12}}}^2 & c_{\rm{\chic{12}}}s_{\rm{\chic{12}}}c_{\rm{\chic{13}}}^2 & c_{\rm{\chic{12}}}c_{\rm{\chic{13}}}s_{\rm{\chic{13}}} \\
c_{\rm{\chic{12}}}s_{\rm{\chic{12}}}c_{\rm{\chic{13}}}^2 & s_{\rm{\chic{12}}}^2c_{\rm{\chic{13}}}^2 & s_{\rm{\chic{12}}}c_{\rm{\chic{13}}}s_{\rm{\chic{13}}} \\
c_{\rm{\chic{12}}}s_{\rm{\chic{13}}}c_{\rm{\chic{13}}} & s_{\rm{\chic{12}}}c_{\rm{\chic{13}}}s_{\rm{\chic{13}}} & s_{\rm{\chic{13}}}^2
\end{array}
\!\right)\!.
\end{equation}
The first term in  Eq. \eqref{eq:Hmatter} is the Hamiltonian that drives the flavor evolution in vacuum, while
the second term accounts for the matter effects. For antineutrinos, the sign of the second term is reversed
and the matrix $\rm\Gamma$ is replaced by its complex conjugate.

In normal matter (n, p, e)
\begin{equation}
{\mathtt V}(\ell) =  \sqrt 2 G_F n_e(\ell)\,,
\end{equation}
where $n_e(\ell)$ is the electron number density and $G_F$ is the Fermi constant.
The electron number density depends on both the matter density $\rho(\ell)$ and the chemical and isotopic composition
of the medium:
\begin{equation}
\label{elecdensity}
n_e (\ell) =\dfrac{{\rho}(\ell)}{m_{\rm u}} \dfrac{Z}{A}( \ell)\,,
\end{equation}
where $m_{\rm u} = 931.494$ MeV is the atomic mass unit and
$
Z/A =  \sum_{\lambda} \digamma_{\!\!\lambda} (Z/A)_{\lambda} .
$
The summation runs over all the chemical elements present in the medium and $(Z/A)_{\lambda}$ denotes the ratio between the atomic number $Z_{\lambda}$ and the atomic mass $A_{\lambda}$ of the element that contributes a fraction $\digamma_{\!\!\lambda} $ to the mass at a given position.

The relevant quantities for us are the oscillation probabilities for ${\nu_\mu} (\bar{\nu}_\mu) \to {\nu_{\mu.\tau}} (\bar{\nu}_{\mu,\tau})$
and ${\nu_e} (\bar{\nu}_e) \to {\nu_{\mu,\tau}} (\bar{\nu}_{\mu,\tau})$, which are required to compute the $\mu$-like events produced in a detector by
``upward'' atmospheric neutrinos after traveling a distance $L = 2 R_{_{\hsn\oplus}} \cos \eta $ through the Earth (see Eq. \eqref{eq:numeroeventos}).
According to the previous considerations, because of the dependence of these probabilities on the potential ${\mathtt V}(\ell)$, atmospheric neutrinos are sensitive to changes in the density and composition of the traversed layers. In what follows, we examine the feasibility of applying such an effect to obtain  meaningful information about the deepest part of the Earth.  With this in mind, we calculate $\tilde{\mathcal{U}}(L)$ for the Earth modeled as a sphere made up of 55 concentric spherical shells with different (constant) densities, as described in Sec. \eqref{sec:structure}.

The complete evolution operator can be expressed as the product of the evolution operators for the consecutive shells through which the neutrinos pass on their way to the detector:
\begin{eqnarray}
\label{fullevolop}
\tilde{\cal U}(L) = \prod_{j=1}^{2 j_m-1}\tilde{\cal U}_j(L_j)\,,
\end{eqnarray}
where $L_J$ are the distances in each shell, such that
$L = \sum_j L_j$
and
\begin{equation}
\tilde{\cal U}_j(L_j) =
\exp(-i{\tilde H}_j L_j)\,.
\end{equation}
The effective potential in ${\tilde H}_{j}$ takes the fixed value $ {\mathtt V}_{j} = \sqrt{2}\hspace{.03cm} G_{\hs{-.03cm}F} n^{j}_e,$ where  $n^{j}_e$ denotes the constant number density of electrons in shell $j$.  Since the Hamiltonians for the different layers do not generally commute between them, the exponential factors  in Eq. \eqref{fullevolop} must be in the prescribed order. Each of these factors can be evaluated by applying the Cayley-Hamilton theorem which allows us to convert the infinity series into a polynomial:
$
\tilde{\mathcal U}_{j}(L_{j})  = a^j_0 I+a^j_1 {\tilde H}_{\ell} + a^j_2 {\tilde H}_{\ell}^2 \,,
$
where the coefficients are functions of the eigenvalues of ${\tilde H}_{j}$ (which coincide with those of $H_{j}$)
\cite{Ohlsson:1999xb, DOlivo:2020ssf}.

\begin{table}[h]
\centering
\begin{tabular}{@{}|c|c|c|c|@{}}
\toprule
\textbf{Set} & \textbf{$f_{oc}$} & \textbf{$f_{(Z/A)_{oc}}$} & \textbf{$f_{(Z/A)_{M_1}}$} \\ \midrule
I   & 1    & 1    & 1    \\
II  & 1    & 1.01 & 1    \\
III & 1.01 & 1    & 1    \\
IV  & 1    & 1    & 1.01 \\ \bottomrule
\end{tabular}
\caption{Different sets of density and composition.}
\label{tab:sets}
\end{table}

The distances $L_j$ are functions of the nadir angle $\eta$. According to Fig. \ref{fig:earth}, they
can be determined as
$ L_j=\vert \tilde{z}_{j+1} - \tilde{z}_j \vert$,
where
\begin{equation}
 \tilde{z}_j = \begin{cases}
z^+_j\,,  & 1 \le j \le j_m\,,\\
z^-_{2j_m +1 -j}\,, & j_m + 1 \le j \le 2j_m\,,
\end{cases}
\end{equation}
\begin{equation}
z^{\pm}_ j = R_{_{\hsn\oplus}} \cos\eta\pm \sqrt{r_j^2-(R_{_{\hsn\oplus}} \sin\eta)^2}\,.
\end{equation}
\begin{figure}[h!]
\centering
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{probmumu_a20_NO.pdf}
\caption{NO and $\eta = 20^{\circ}$}
\end{subfigure}\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{probmumu_a50_NO.pdf}
\caption{NO and $\eta = 50^{\circ}$}
\end{subfigure}
\caption{Survival probability $P_{{\nu_\mu} \to {\nu_\mu}}$ as a function of the neutrino energy, for normal ordering (NO) and nadir angle $\eta$ equal to (a) $20^{\circ}$
and (b) $50^{\circ}$\!.}
\label{fig:mu--mu}
\end{figure}
\begin{figure}[h!]
\centering
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{probemu_a20_NO.pdf}
\caption{NO and $\eta = 20^{\circ}$}
\end{subfigure}\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{probemu_a20_IO.pdf}
\caption{IO and $\eta = 20^{\circ}$}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{probemu_a50_NO.pdf}
\caption{NO and $\eta = 50^{\circ}$}
\end{subfigure}\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{probemu_a50_IO.pdf}
\caption{IO and $\eta = 50^{\circ}$}
\end{subfigure}
\caption{Transition probability $P_{{\nu_e} \to {\nu_\mu}}$ as a function of the neutrino energy, for normal ordering (NO) and inverted ordering (IO) and nadir angle $\eta$ equal to (a,b) $20^{\circ}$ and (c,d) $50^{\circ}$\!.}
\label{fig:mu--e}
\end{figure}
To exemplify the effect that changes in composition and density have on the probabilities of flavor oscillations, in Figs.~\ref{fig:mu--mu} and \ref{fig:mu--e} we show the probabilities for $\nu_{\mu} \to \nu_{\mu}$ and $\nu_{e}\to \nu_{\mu}$ as functions of the neutrino energy, for two different angles of incidence, both for NO and IO.  In the next section, we use these probabilities in the calculation of  the $\mu$-like events produced by atmospheric neutrinos arriving at the detector after passing through the internal regions of the Earth.

\section{Neutrino events and test of Earth's composition}
\label{sec:simu}

Let us consider a generic detector with a mass of water or ice containing a number $n_{\chic N}$ of nucleons as a target. Such a detector, like the planned ORCA, will efficiently detect the Cherenkov radiation emitted along the path of $\mu^{\pm}$ produced by the interactions of $\nu _ {\mu}$ and $\bar\nu_{\mu} $ with the nucleons in the instrumented volume or near it. In addition, it will be posible to detect showers due to the $e$ or $\tau$ produced by the associated neutrinos \cite{Stavropoulos:2021hir,Kalaczynski:2021ytv,Bourret:2017tkw}.

To generate a data set that simulates the actual observations, with the inevitable statistical fluctuations, we perform a numerical calculation of the events due to neutrinos.  Our observable is the number of $\mu$-like events ${\mathcal N}_{\mu}$ and in what follows we use it to test different hypothesis about the composition and density of the outer core and the lower part of the mantle.
To identify the most sensitive regions, we first do a scan dividing the cone under the detector into a series of angular and energy bins and calculate the number of $\mu$-like events within given angular and energy intervals:
\begin{eqnarray}
\label{eq:numeroeventos}
 {\mathcal N}_{\mu}= n_{\chic N} T \int_{E_{\text{min}}}^{E_{\text{max}}} \! \dd{E} \int_{\eta_{\text{min}}}^{\eta_{\text{max}}} \dd{\Omega}
\left[\sigma^{cc}_{\nu_{\mu}}(E)\left(P_{\nu_{\mu}\rightarrow \nu_{\mu}} {\dfrac{d\Phi}{dE}}^{\!\!\nu_{\mu}}
\!+ P_{\nu_{e}\rightarrow \nu_{\mu}}{\dfrac{d\Phi}{dE}}^{\!\!\nu_{e}}\right) +\right.
\nonumber \\
\left. \sigma^{cc}_{\bar \nu_{\mu}}(E)\left(P_{\bar \nu_{\mu}\rightarrow \bar \nu_{\mu}} {\dfrac{d\Phi}{dE}}^{\!\!\bar \nu_{\mu}}
\!+ P_{\bar \nu_{e}\rightarrow \bar \nu_{\mu}} {\dfrac{d\Phi}{dE}}^{\!\!\bar \nu_{e}} \right) + \right.
\nonumber \\
\left. \sigma^{cc}_{\nu_{\tau}}(E) Br_{\tau \rightarrow \mu} \left(P_{\nu_{\mu}\rightarrow \nu_{\tau}} {\dfrac{d\Phi}{dE}}^{\!\! \nu_{\mu}}
\!+ P_{\nu_{e}\rightarrow \nu_{\tau}} {\dfrac{d\Phi}{dE}}^{\!\! \nu_{e}} \right) + \right.
\nonumber \\
\left. \sigma^{cc}_{\bar \nu_{\tau}}(E) Br_{\bar\tau \rightarrow \bar\mu} \left(P_{\bar \nu_{\mu}\rightarrow \bar \nu_{\tau}} {\dfrac{d\Phi}{dE}}^{\!\!\bar \nu_{\mu}}
\!+ P_{\bar \nu_{e}\rightarrow \bar \nu_{\tau}} {\dfrac{d\Phi}{dE}}^{\!\!\bar \nu_{e}} \right)
 \right],
\end{eqnarray}
where $T$ is the detection time. The fluxes of atmospheric neutrinos and antineutrinos were taken from Ref. \cite{Atmnu:1996},
while the charged-current cross sections for the ${\nu}_{\mu}(\bar{\nu}_{\mu})$- nucleon and ${\nu}_{\tau}(\bar{\nu}_{\tau})$-nucleon scatterings, in the considered range of neutrino energies ($ 1-10$ GeV), are given approximately by \cite{Formaggio:2013kya}:
\begin{equation}
\begin{aligned}
\sigma^{cc}_{\nu_{\mu}}(E) &\simeq 0.75 \times 10^{-38} \left({E}/{\rm GeV}\right) {\rm cm}^2\,, \\
\sigma^{cc}_{\bar \nu_{\mu}}(E) &\simeq 0.35 \times 10^{-38} \left({E}/{\rm GeV}\right) {\rm cm}^2\,,\\
\sigma^{cc}_{\nu_{\tau}}(E) &\simeq 0.13 \times 10^{-38} \left({E}/{\rm GeV}\right) {\rm cm}^2\,, \\
\sigma^{cc}_{\bar \nu_{\tau}}(E) &\simeq 0.05 \times 10^{-38} \left({E}/{\rm GeV}\right) {\rm cm}^2 \,.
\end{aligned}
\end{equation}
In Eq. \eqref{eq:numeroeventos}, $Br_{\tau \rightarrow \mu} = Br_{\bar\tau \rightarrow \bar\mu} \simeq 0.17$ are the branching ratios of the decays $\tau \rightarrow \mu {\bar \nu}_\mu \nu_\tau$ and $\bar\tau \rightarrow {\bar\mu}{\nu_\mu} {\bar\nu}_\tau$. The dependence of ${\mathcal N}_{\mu}$ on the density and composition of the medium is incorporated through the oscillation probabilities, which are calculated from the expressions given in Sec. \ref{sec:osc}. It is understood that these probabilities are evaluated at $L$.
The values of the oscillation parameters are those given in Table \ref{tab:constants}. To estimate the impact of current uncertainties on the oscillation parameters, in addition to the best fits, we consider their values at the extremes of the 1$\sigma$ ranges \cite{Capozzi:2020,Esteban:2020cvm,deSalas:2020pgw}. Thus, the muon numbers were also calculated by evaluating the oscillation probabilities at those values of a given parameter and keeping the best fits for the rest. We find that errors in the oscillation parameters introduce only few percent variations in the number of muon events and have little effect on our analysis.

To determine the angular and energy intervals where ${\mathcal N}_{\mu}$ is more sensitive to changes in the compositions of either the outer core or the lower mantle we introduce the quantity
\begin{equation}
\label{eq:upsilon}
{\Upsilon}(E,\eta) \equiv \abs{1 - \dfrac{\mathcal{N}_{\mu}^{(Z/A)}}{\mathcal{N}_{\mu}^0}} \times 100\,,
\end{equation}
which gives the percentage difference between the number of events for the standard composition $\mathcal{N}_{\mu}^{0}$ and the number of events for a different composition $\mathcal{N}_{\mu}^{\rm (Z/A)}$. The radii of the layers and the densities of the shells in the layers are those given from the PREM.

\begin{figure}[h!]
\centering
\begin{subfigure}[b]{.50\textwidth}
\includegraphics[width=\textwidth]{SNO_101_1.png}
\caption{Normal Ordering}
\label{fig:SNO_101_1}
\end{subfigure}\begin{subfigure}[b]{.50\textwidth}
\includegraphics[width=\textwidth]{SIO_101_1.png}
\caption{Inverse Ordering}
\label{fig:SIO_101_1}
\end{subfigure}
\caption{Level surfaces of the function $\Upsilon(E,\eta)$ (Eq. \eqref{eq:upsilon}) in the ($E,\eta$) plane, for a $1\%$ change in the composition of the outer core.}
\label{fig:sens_OC}
\end{figure}

Figs.~\ref{fig:sens_OC} and \ref{fig:sens_D2} show the level surfaces of $\Upsilon$ in the ($E,\eta$) plane
for a change of $1\%$ in the composition of the outer-core and the composition of the $M_1$ region, respectively, for both NO and IO.  From the figures, it is apparent that in the case of the outer core the most sensitive regions correspond to energies around
5 GeV and angles compatible with the shadow of the outer core and $M_1$.  On the other hand, for the lower mantle the angular region has to be increased and the energy shifts to slightly lower values. The deviations in the number of events are considerably more pronounced for NO, mainly because for IO the resonance effects in matter occur for antineutrinos, whose charged cross-section is about two times smaller than that of neutrinos.

To simulate the number of events observed by a km$^3$ detector as ORCA, we follow the same procedure implemented in \cite{DOlivo:2020ssf}, which we reproduce here for completeness.
\begin{figure}[h!]
\centering
\begin{subfigure}[b]{.50\textwidth}
\includegraphics[width=\textwidth]{SNO_1_101.png}
\caption{Normal Ordering}
\label{fig:SNO_1_101}
\end{subfigure}\begin{subfigure}[b]{.50\textwidth}
\includegraphics[width=\textwidth]{SNO_1_101.png}
\caption{Inverse Ordering}
\label{fig:SIO_1_101}
\end{subfigure}
\caption{Level surfaces of the function $\Upsilon(E,\eta)$ (Eq. \eqref{eq:upsilon}) in the ($E,\eta$) plane, for a $1\%$ change in the composition of layer $M_1$.}
\label{fig:sens_D2}
\end{figure}
In order to perform the Monte Carlo simulation, based on the previous results, we consider events with $4 \,\text{GeV} < E < 11\, \text{GeV}$ and $10\degree<\,\eta\,<60\degree$. Both of these intervals are divided into $200$ subintervals. Every pseudo experiment is made up by tossing, in each square bin of the grid, a number of Poisson distributed events with the mean value equal to ${\mathcal N}_{\mu}$ as given by Eq.~(\ref{eq:numeroeventos}). Thus, each of the $n_{exp}$ pseudo experiments consists of 200 $\times$ 200 numbers corresponding to events, one for each bin. This sample of events are then distributed in angle and energy. We call them the  {\it true events} and suppose that they are distributed according to the probability distribution function (pdf)  $f^{\,i_{exp}}_t(E,\eta)$, ${i_{exp}} = 1, \cdots, n_{exp}$. For this function we take the normalized histograms constructed by means of the Monte Carlo simulation. The fractional number of true events in bin $(i,j)$ for the $i_{exp}$-th experiment is given by the integral of $f^{\,i_{exp}}_t(E,\eta)$ over the energy and angle intervals of the bin $(i,j)$.

To obtain a realistic distribution of events we must allow for a limited resolution of the detector, both in energy and angle. The net effect is the redistribution of the true events ("smearing") in the grid bins, which is implemented by folding the true distribution with a resolution function $\mathcal{S}(E_{\mathtt o}, \eta_{\mathtt o} \vert E, \eta)$. We assume a Gaussian smearing and write
\begin{eqnarray}
\label{smearing}
\mathcal{S}(E_{\mathtt o}, \eta_{\mathtt o} \vert E, \eta)=\dfrac1{2\pi \Delta \eta(E) \Delta E(E)}\exp\left(-\dfrac{(\eta - \eta_{\mathtt o})^2}{2\Delta\eta(E)^2}\right) \exp\left(-\dfrac{(E-E_{\mathtt o})^2}{2\Delta E(E)^2}\right),
\end{eqnarray}
with the detector characterized by the angular and energy resolutions $\Delta \eta(E) = \alpha_{\eta}/\sqrt{E/\text{GeV}}$ and $\Delta E(E) = \alpha_{E} \, E$, respectively. For the values of the parameters $\alpha_{\eta}$ and $\alpha_{E}$ we consider different situations according to the discussion in  Ref. \cite{Rott:2015kwa}. To keep our analysis as simple as possible we assume a detection efficiency of $100\%$.
When convoluted with the true events the kernel in Eq. \ref{smearing} gives us what we call the {\it observed events}. That is, $\mathcal{S}(E_{\mathtt o}, \eta_{\mathtt o} \vert E, \eta)$ represents the conditional pdf for the measured values to be $(E_{\mathtt o}, \eta_{\mathtt o})$ if the true values were $(E, \eta)$. Since the event is observed somewhere, this function is normalized such that
\begin{equation}
\iint \! {\mathcal S}(E_{\mathtt{o}}, \eta_{\mathtt{o}} \vert E, \eta) \dd{E}_{\mathtt{o}} \dd{\eta}_{\mathtt{o}} = 1\,.
\end{equation}

In terms of the resolution function,  the number of observed events $\mathcal{O}^{\,i_{exp}}_{m,n}$ in the bin $(m,n)$ for the $i_{exp}$-th experiment is given by
\begin{equation}
\label{observed}
\mathcal{O}^{\,i_{exp}}_{m,n} = {\mathcal N}_{tot}\sum\limits_{i,j} \iint_{\text{bin} \; m,n} \! \dd{E}_{\mathtt{o}} \, \dd{\eta}_{\mathtt{o}} \; \iint_{\text{bin} \; i,j} \! \dd{E} \, \dd{\eta}\;
\mathcal{S}(E_{\mathtt o}, \eta_{\mathtt{o}} \vert E, \eta) f^{\,i_{exp}}_t(E, \eta)\,.
\end{equation}
As an illustration, in Fig. \ref{fig:trueobsNO} and Fig. \ref{fig:trueobsIO} we show the pdf of the true and observed events as functions of the energy and nadir angle, for NO and IO, respectively, in the case of a {\it standard Earth}. The figures were made using the same (large) number of bins for both kind of events, but these numbers generally differ. In what follows, to test how well the hypothesis of different Earth compositions is in agreement with the standard Earth, we take the observed events to be distributed into five angular bins and nine energy bins.
\begin{figure*}[h!]
\centering
\begin{subfigure}[b]{0.50\textwidth}
\includegraphics[width=\textwidth]{eve_ver_NO.png}
\caption{}
\label{fig:verNO}
\end{subfigure}\begin{subfigure}[b]{0.50\textwidth}
\includegraphics[width=\textwidth]{eve_obs_NO.png}
\caption{}
\label{fig:obsNO}
\end{subfigure}
\caption{Probability distribution functions of (a) true events and (b) observed events, as functions of energy and nadir angle, for the {\it standard Earth} and  normal ordering of neutrino masses. The same number of bins was used to make both figures, but a much small number of bins was used to calculate the observed events.}
\label{fig:trueobsNO}
\end{figure*}
\begin{figure*}[h!]
\centering
\begin{subfigure}[b]{0.50\textwidth}
\includegraphics[width=\textwidth]{eve_ver_IO.png}
\caption{}
\label{fig:verIO}
\end{subfigure}\begin{subfigure}[b]{0.50\textwidth}
\includegraphics[width=\textwidth]{eve_obs_IO.png}
\caption{}
\label{fig:obsIO}
\end{subfigure}
\caption{Probability distribution functions of (a) true events and (b) observed events, as functions of energy and nadir angle, for the {\it standard Earth}  and  inverted ordering of neutrino masses. The same number of bins was used to make both figures, but a much small number of bins was used to calculate the observed events.}
\label{fig:trueobsIO}
\end{figure*}
From Eq. \ref{observed}, for each bin $(m,n)$, with  $m = 1, \cdots, 9$ and $n = 1, \cdots, 5$, we determine the number of events $\mathcal{O}^{\hspace{0.01 cm}{\mathtt s}i_{exp}}_{\alpha}$ for a standard Earth and  $\mathcal{O}^{{\hspace{0.02 cm}i_{exp}}}_{\alpha}\!({\rm H}, \varkappa)$ for an alternative Earth with different composition and density. Here $\alpha = 1, \cdots , 45$ label the two-dimensional bins ($m,n$). Thus, for each of these bins we have a sample of the observed events and determinate the corresponding mean values
\begin{eqnarray}
\mathcal {\bar O}^{\mathtt s}_{\alpha}&=&\frac1{n_{exp}}\sum_{i_{exp}}\mathcal{O}^{\hspace{0.01 cm}{\mathtt s}i_{exp}}_{\alpha}\,, \nonumber\\
\mathcal{\bar O}_{\alpha}({\rm H}, \varkappa)&=&\frac1{n_{exp}}\sum_{i_{exp}}\mathcal{O}^{{\hspace{0.02 cm}i_{exp}}}_{\alpha}\!({\rm H}, \varkappa)\,.
\end{eqnarray}

For Poisson distributed events, in terms of the likelihood function $L$ we construct the negative log-likelihood ratio function as
\begin{eqnarray}
\chi^2_{\lambda} = -\ln\left[\frac{L\left(\mathcal{\bar O}^{\mathtt s}_{\alpha};\mathcal{\bar O}_{\alpha}({\rm H}, \varkappa)\right)}{L\left(\mathcal{\bar O}^{\mathtt s}_{\alpha},\mathcal{\bar O}^{\mathtt s}_{\alpha}\right)}\right]=
2\sum_{\alpha} \left[\mathcal{\bar O}_{\alpha}({\rm H} , \varkappa)-\mathcal{\bar O}^{\mathtt s}_{\alpha} \ln\left(\frac{\mathcal{\bar O}^{\mathtt s}_{\alpha}}{\mathcal{\bar O}_{\alpha}({\rm H} , \varkappa)}\right)\right]\,.
\end{eqnarray}
According to Wilks' theorem \cite{wilks1938} the $\chi^2_{\lambda}$ distribution can be approximated  by the $\chi^2$ distribution and, from it, the goodness of the fit can be established. The statistical significance of the $\chi^{2}$ test is
given, as usual, by the $\mathbf p$-value:
\begin{equation}
{\mathbf p}=\int_{\chi^2}^\infty \! f_{\chi}(w,n_{_{\!\rm dof}}) \dd{w}\,,
\end{equation}
where $n_{_{\!\rm dof}}$ is the number of degrees of freedom and
$ f_{\chi}(w,n_{_{\!\rm dof}})$
is the chi-square distribution. In our case, $n_{_{\!\rm dof}} = 9\times5 -2 = 43$.
In this way, we can examine the levels of discrepancy between the standard Earth and different hypothesis about the composition and density of the outer core and the lower mantle. This is done in the next section, where the results and final comments are presented.

\section{Results and Final Comments}
\label{sec:results}

Our goal is to determine to what extent the presence of light elements, in particular hydrogen, in the outer core and mantle can produce measurable effects due to the modifications it introduces in the flavor transformations of atmospheric neutrinos. We also pay attention to how uncertainties in the densities of the deepest regions of the planet can complicate obtaining compositional information. As discussed in Section 2, we rescale the mantle and outer core densities so that the constraints imposed by the Earth's total mass and moment of inertia are satisfied. At the same time, we allow for some changes in the compositions of the outer core and/or the lower mantle $M_1$.

As the quantities to be fitted, we take the compositions of the outer-core and the lower mantle $M_1$ and the densities of the outer core, $M_1$ and the rest of the mantle $M_2$. To evaluate the effects that changes of two of these quantities have on our observable, we construct regions in several planes:
($\delta_{\chic{OC}},(Z/A)_{\chic{OC}}$), ($\delta_{M_1},(Z/A)_{M_1}$), and ($(Z/A)_{M_1},(Z/A)_{\chic{OC}}$),  where the statistical significance, given by the ${\mathbf p}$-value, for the discrepancy between the standard and the modified Earth is less than $1\sigma$. In Fig.~(\ref{fig:compo_dens_NO}) and (\ref{fig:compo_dens_IO}) we show these regions for the oscillation parameters given in Table \ref{tab:constants}, for 10 years operation of an 8 Mton detector as ORCA, with resolution parameters $\alpha_{\eta} =0.25$ and $\alpha_{E} = 0.2$. As can be seen, the regions for NO are more restricted than for IO. We consider two different values for the radius of the $M_1$ layer: 5600 km and 5871 km, which correspond to the lower and upper border of the transition zone between the lower and upper mantle.

\begin{figure}[h!]
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{cOCdOCNO_del.pdf}
\caption{Composition and density of the outer-core. }
\label{fig:c_oc_d_oc_NO}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{cM1dM1NO_del.pdf}
\caption{Composition and density of the $M_1$ layer. }
\label{fig:c_M1_d_M1_NO}
\end{subfigure}
\vskip\baselineskip
\centering
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{ccNO.pdf}
\caption{Outer-core and $M_1$ compositions.}
\end{subfigure}
\caption{Expected 1$\sigma$ regions for combined composition and density measures in the outer-core and  the $M_1$ region, for normal ordering and three values of the $M_1$ radii.}
\label{fig:compo_dens_NO}
\end{figure}

In the graphs of composition versus density, non-zero correlations are observed, indicating that increases in composition are compensated by decreases in density, and vice versa, such that the electron density in the shell does not change. This compensation is not complete since, for the Earth's moment of inertia and mass to remain fixed, variations in density must also occur in other layers traversed by neutrinos and, as a consequence, the allowed regions become closed. The 1$\sigma$ confidence regions for combined measures of Z/A in the outer core and $M_1$  agree with those obtained by \cite{VanElewyck:2017dga}.
Suppose that the variation of  Z/A in the outer core or mantle is associated only with the abundance of hydrogen, then, as show in Fig.\ref{fig:ZA}, $Z/A=(1-\digamma_{\!H})(Z/A)_0+ \digamma_{\!H}(Z/A)_H$, where $(Z/A)_0$ is the value with no hydrogen and $\digamma_{\!H}$ is the fractional contribution that hydrogen makes to the total mass of the layer. In this case, the 1$\sigma$ regions are compatible with 5-8 wt\% hydrogen in the lower mantle and outer core, respectively, values too high in light of geophysical estimates.

\begin{figure}[h!]
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{cOCdOCIO_del.pdf}
\caption{Composition and density of the outer-core. }
\label{fig:c_oc_d_oc_IO}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{cM1dM1IO_del.pdf}
\caption{Composition and density of the $M_1$ layer. }
\label{fig:c_M1_d_M1_IO}
\end{subfigure}
\vskip\baselineskip
\centering
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{ccIO.pdf}
\caption{Outer-core and $M_1$ composition.}
\end{subfigure}
\hfill
\caption{Expected 1$\sigma$ regions for combined composition and density measures in the outer-core and the $M_1$ region, for inverted ordering and three values of the $M_1$ radii.}
\label{fig:compo_dens_IO}
\end{figure}

As a general rule, for the non fitted parameters we have kept their PREM values. However, for the sake of completeness, we paid some attention to effects from uncertainties in the densities when fitting the compositions. Thus, we have jointly fitted $(Z/A)_{oc}, (Z/A)_{M_1}$ and $\delta_{oc}$ for the outer core density. (Notice, that the densities of the layers $M_1$ and $M_2$ vary in order to maintain the values of the Earth's mass and moment of inertia, as discussed in Section \ref{sec:structure}). The 1$\sigma$ volume obtained in this way is projected onto the composition-composition plane, giving rise to the shaded area in Fig. \ref{fig:merge}. The observed negative correlation is caused by the compensation effect in matter neutrino oscillations, with an increase in composition counterbalanced by a decrease in density, and vice versa, to keep the electron number density unchanged. In the case of the composition-density plots (Fig.\ref{fig:c_oc_d_oc_NO}, Fig.\ref{fig:c_M1_d_M1_NO}, Fig.\ref{fig:c_oc_d_oc_IO}, and Fig.\ref{fig:c_M1_d_M1_IO}), the uncertainty in the outer core density is incorporated in the lower mantle density, since they are related by Eq. \eqref{reldeltaocdeltaM12}.

\begin{figure}[h!]
\centering
\includegraphics[width=0.75\textwidth]{ccNOmerge.pdf}
\hfill
\caption{The same region of Fig.\ref{fig:compo_dens_NO}\hs{0.03 cm}c, but including the uncertainties in the outer-core densities.
}
\label{fig:merge}
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[width=0.75\textwidth]{ZA.pdf}
\hfill
\caption{The ratio Z/A as a function of the weight fraction of hydrogen.}
\label{fig:ZA}
\end{figure}
As is evident from Figs. \ref{fig:compo_dens_NO} and \ref{fig:compo_dens_IO}, an experiment like ORCA has a limited potential to reveal a non-standard composition and/or density of the Earth. This worsens at 2 and 3$\sigma$ and for the inverted ordering. A question then arises: how much exposure and how much resolution are required to constraint a non-standard composition? As a partial answer, in Fig.\ref{fig:resolution} we shown the1$\sigma$ regions for different exposures (in unity of ten years of ORCA operation) and angular and energy resolutions equal to $\alpha$. We see that with an exposition of thirty years and resolutions of $0.1$ it is possible to constraint the hydrogen content to 1\%.
\begin{figure}[h!]
\centering
\includegraphics[width=0.75\textwidth]{ccNO_H.pdf}
\caption{Fraction of hidrogen in the outer-core and the lower mantle for extra exposition $\mathtt e$ and resolution $\alpha$. }
\label{fig:resolution}
\end{figure}

Finally, we pay special attention to the $D^{\prime \prime}$ region. Since the thickness of this remote interface between the rocky mantle and the iron core is relatively thin, changes in density and composition have little effect on our observable. Fig.\ref{fig:d2comp} shows the $1\sigma$ regions for the detector size versus the composition $D^{\prime \prime}$, for a thickness of 200 km and various values of the detector resolution. As can be seen, the allowed regions are not very restrictive for reasonable values of resolution and size of the detector. The ability to constrain the composition is significantly improved in the case of a slightly thicker layer (500 km).

\begin{figure}[h!]
\centering
\begin{subfigure}[b]{0.75\textwidth}
\includegraphics[width=\textwidth]{D2_25.pdf}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.75\textwidth}
\includegraphics[width=\textwidth]{D2_26.pdf}
\end{subfigure}
\caption{$D^{"}$ composition for different resolutions and size detector. $D^{"}$ thickness $500$ km.}
\label{fig:d2comp}
\end{figure}

In summary, we have studied the possibility of conducting an oscillation tomography of the Earth based on the matter effects on the flavor oscillations of atmospheric neutrinos  propagating through the depths of the planet. Using the $\mu$-like events in a generic large Cherenkov detector as physical observables and making a Monte Carlo simulation of the energy and azimuthal angle distribution of these event, we tested possible variants with respect to a reference geophysical model with the densities as given by PREM and different composition in the outer-core and lower mantle. Unlike previous studies, the procedure we followed in this work allowed us to simultaneously vary the composition and density. When one of these quantities was fixed in the value of the reference geophysical model, our results, shown in Figs. \ref{fig:compo_dens_NO}, are compatible with those obtained by others authors in an uncorrelated way \cite{VanElewyck:2017dga,Borriello:2009ad}. As shown, the study of questions of the type examined here would benefit greatly from the application of phenomena normally studied in the realm of particle physics and the advent of new and improved neutrino telescopes, such as KM3NeT, ORCA, PINGO and HyperK \cite{Adrian-Martinez:2016fdl,Fermani:2019vac,Hyper-Kamiokande:2018ofw,Aarsten:2017}.

\begin{acknowledgements}
This work was supported in part by DGAPA-UNAM Grant No. PAPIIT IN106322 and by the Consejo Nacional de Investigaciones Cientificas y Tecnicas (CONICET) and the UNMDP Mar del Plata, Argentina.
\end{acknowledgements}

\newpage
\bibliographystyle{spphys}


text:
\section{Introduction}

Supernarrow dibaryons (SNDs) are 6-quark states, a decay of which into
two nucleons is forbidden by the Pauli exclusion principle
\cite{fil1,fil2,alek1}.
Such states satisfy the following condition:
\begin{equation}
(-1)^{T+S}P=+1
\end{equation}
where $T$ is the isospin, $S$ is the internal spin, and $P$ is the
dibaryon parity. In the $NN$ channel, these six-quark states would
correspond to the following forbidden states:
even singlets and odd triplets with the isotopic spin $T=0$ as well as
odd singlets and even triplets with \mbox{$T=1$}.
These six-quark states with the masses \mbox{$M < 2m_{N}+m_{\pi}$}
($m_N (m_{\pi})$ is the nucleon (pion) mass) can mainly decay
by emitting a photon. This is a new class of metastable six-quark states
with the decay widths \mbox{$< 1$keV}.

The experimental discovery of the SNDs would have important consequences
for particle and nuclear physics and astrophysics. This would lead to
a deeper understanding of the evolution of compact stars and the new
possibility of quark-gluon plasma observation.
In nuclear physics there would be a new concept: SND-nuclei.

In the framework of the MIT bag model, Mulders et al. \cite{muld} calculated
masses of different dibaryons, in particular, masses of $NN$-decoupled
dibaryons. They predicted dibaryons $D(T=0;J^P =0^{-},1^{-},2^{-};M=2.11$ GeV)
and $D(1;1^{-};2.2$ GeV) corresponding to the forbidden states $^{13}P_J$
and $^{31}P_1$ in the $NN$ channel.
However, the dibaryon masses obtained exceed the pion production threshold.
Therefore, these dibaryons preferentially decay into the $\pi NN$ channel
and their decay widths are larger than 1 MeV.

Using the chiral soliton model, Kopeliovich \cite{kop} predicted that
the masses of \mbox{$D(1,1^+)$} and $D(0,2^+)$ SNDs could
exceeded the two nucleon mass by 60 and 90 MeV, respectively.
These values are lower than the pion production threshold.

In the framework of the canonically quantized biskyrmion model,
Krupnovnickas {\em et al.} \cite{riska} obtained an indication on
possibility
of the existence of one dibaryon with $J=T=0$ and two dibaryons with $J=T=1$
with masses smaller than $2m_N+m_{\pi}$.

In the present paper we analyze the experimental search for SNDs, construct
and analyze the mass formula for SNDs, and suggest a possible interpretation
of exotic narrow baryons with low masses.

\section{Supernarrow dibaryons}
\label{sec:2}

We will consider the following SNDs:
$D(T=1,J^P=1^+,S=1)$ and $D(1,1^-,0)$.

It is worth noting, that the  state $(T=1, J^P=1^-)$ corresponds
to the states $^{31}P_1$ and $^{33}P_1$ in the NN channel.
The former is forbidden and
the latter is allowed for a two-nucleon state. In our work we will study
the dibaryon $D(1,1^-,0)$, a decay of which into two nucleons is forbidden
by the Pauli principle (i.e. $^{31}P_1$ state).

SNDs can be formed in processes of interaction with the deuteron
only if the  nucleons in the deuteron overlap sufficiently, such that a
6-quark state with deuteron quantum numbers can be formed. In this case, an
interaction of a photon or a meson with this state can change its
quantum numbers so that a metastable state can form.
Therefore, the probability of the production of such dibaryons is
proportional to the probability $\eta$ of the 6-quark state existing
in the deuteron.

The magnitude of $\eta$ can be estimated from the
deuteron form factor at large $Q^2$ (see for example \cite{bur,grach}). However,
the values obtained depend strongly on the model of the form factor of
the 6-quark state over a broad region of $Q^2$. Another way to estimate
this parameter is to use the discrepancy between the theoretical and
experimental values of the deuteron magnetic moment \cite{kim,bhad,kon2}.
This method is free from the restrictions quoted above and gives
$\eta\le 0.03$ \cite{kon2}.

\begin{wrapfigure}{h}{0.3\textwidth}
\epsfxsize=0.3\textwidth
\epsfysize=2cm
\epsffile{fig1.eps}
\caption{The diagram of the SND decay into $\g NN$}
\label{dec}
\end{wrapfigure}
Since the energy of nucleons, produced in the decay of the SNDs under
study with $M<2m_N +m_{\pi}$, is small, it is expected that the main
contribution to a two nucleon system should come from the
$^{31}S_0$ (virtual singlet) state (Fig. \ref{dec}).

\begin{table}
\label{table1}
\centering
\caption{ Decay widths of the dibaryons $D(1,1^+,1)$ and $D(1,1^-,0)$
at various dibaryon masses $M$.}
\begin{tabular}{|c|c|c|c|c|c|c|}\hline
$M$(GeV)     &1.904  &1.926 &1.942 &1.965 &1.985  &2.006  \\ \hline
$\a_t(1,1^-)$&0.0514 &0.327 & 0.771&1.909 & 3.495 &5.073  \\
(eV)         &       &      &      &      &       &       \\ \hline
$\a_t(1,1^+)$&0.206  &1.307 &3.083 &7.635 &13.98  &23.49  \\
(eV)         &       &      &      &      &       &        \\ \hline
\end{tabular}
\end{table}
The results of calculations of the
decay widths of the dibaryons into $\g NN$ on the
basis of such assumptions at $\eta=0.01$ are listed in Table I.

As a result of the SND decay through $^{31}S_0$ in
the intermediate state, the probability distribution of such a decay
over the emitted photon energy $\omega$ should be characterized by a
narrow peak at the photon energy close to the maximum value
$\omega_m=(M^2-4m^2_N)/2M$ (Fig.~2).
Note that the interval of the photon energy from $\omega_m$ to
$\omega_m-1$ MeV contains about 75\% of the contribution to the width
of the decay $D(1,1^{\pm})\to \gamma NN$. This leads to a very small
relative energy of the nucleons from the SND decay
and these nucleons are emitted into a narrow angle cone with respect
to the direction of the SND motion.
\begin{figure}[h]
\epsfxsize=6cm
\epsfysize=6cm
\centering{
\epsffile{fig2.eps}}
\caption{The distribution of the decay probability of the isovector
SND over $\omega_m-\omega$ at $M=1904$ MeV.}
\label{width}
\end{figure}

\begin{figure}[ht]
\epsfxsize=6cm
\epsfysize=6cm
\centering{
\epsffile{fig3.eps}}
\caption{The distribution of the decay probability of SNDs over
$\theta_{1,2}$}
\label{d-angl}
\end{figure}
Moreover, the distribution of the SND decay probability over the angle
between the final nucleons should be characterized by a narrow angular
cone.

The dependence of SND decay probability on the angle
between the final nucleons $\theta_{12}$ is shown in Fig.~3.
This figure  demonstrates that the nucleons from
the decay of SNDs should be mainly emitted in a very
narrow angular cone.
These dependencies should be taken into account when looking for SNDs.

\section{Previous works}

For the first time, SNDs have been searched for in the reactions
$pd\to p+pX_1$ and $pd\to p+dX_2$
\cite{konob,izv,yad,prc,conf1,conf2,epj,hadr}.
The experiment was
carried out at the Proton Linear Accelerator of INR with 305 MeV
proton beam using the two-arm mass
spectrometer TAMS.

Several software cuts have been applied to the mass spectra in these works.
In particular, the authors limited themselves by the consideration of
intervals of the proton energy and angles from the decay of the $pX_1$
states, and very narrow angular cone between final nucleons, which
were determined by the kinematics of the SND decay into $\gamma NN$ channel.
Such cuts are very important as it provides a possibility to suppress the
contribution from the background reactions and random coincidences
essentially.

In the works \cite{conf1,conf2,epj,hadr}, CD$_2$ and $^{12}$C were
used as targets.
The scattered proton was detected in the left arm
of the spectrometer TAMS at the angle $\theta_L=70^{\circ}$. The second
charged particle (either $p$ or $d$) was detected in
the right arm by three telescopes located at $\theta_R=34^{\circ}$,
$36^{\circ}$, and $38^{\circ}$.

\begin{figure}[ht]
\epsfxsize=7cm
\epsfysize=9cm
\centering{
\epsffile{fig4.eps}}
\caption{The missing mass $M_{pX_1}$ (a) and
$M_{X_1}$ (b) spectra of the reaction $pd\to p+pX_1$.}
\label{inr}
\end{figure}
As a result,
three narrow peaks in missing mass spectra have been observed
(Fig.\ref{inr}) at $M_{pX_1}=1904\pm 2$, $1926\pm 2$, and $1942\pm 2$ MeV
with widths equal to the experimental resolution ($\sim 5$ MeV) and
with numbers of standard deviations (SD) of 6.0, 7.0, and 6.3, respectively.
It should be noted that the dibaryon peaks at $M=1904$ and 1926 MeV had
been observed earlier by same authors in ref. \cite{prc,konob,izv,yad}
at somewhat different kinematical conditions.

On the other hand, no noticeable
signal of the dibaryons has been observed in the missing mass $M_{dX_2}$
spectra of the reaction $pd\to p+dX_2$.
The analysis of the angular distributions of the protons from the decay of
$pX_1$ states and the suppression observed of the SND decay into $\gamma d$
showed that the peaks found can be explained as a manifestation of the
isovector SNDs, the decay of which into two nucleons is forbidden
by the Pauli exclusion principle.

An additional information about the nature of the observed states has been
obtained by studying the missing mass $M_{X_1}$ spectra of the
reaction $pd\to p+pX_1$.
If the state found is a dibaryon decaying mainly into two nucleons then
$X_1$ is the neutron and the mass $M_{X_1}$ is equal to the neutron mass
$m_n$. If the value of $M_{X_1}$, obtained from the experiment, differs
essentially from $m_n$, then $X_1=\gamma+n$ and we have the additional
indication that the observed dibaryon is SND.

The simulation of the missing mass $M_{X_1}$ spectra of the reaction
$pd\to ppX_1$ has been performed \cite{conf1,conf2,epj,hadr} assuming that
the SND decays as $D\to\g+\, ^{31}S_0\to \g pn$ through two nucleon singlet
state $^{31}S_0$ \cite{fil2,prc,epj}. As a result, three narrow peaks at
$M_{X_1}=965$, 987, and 1003 MeV have been predicted. These peaks
correspond to the decay of the isovector SNDs with masses 1904, 1926, and
1942 MeV, respectively.

In the experimental missing mass $M_{X_1}$ spectrum besides the peak
at the neutron mass due to the process $pd\to p+pn$,
a resonance-like behavior of the spectrum has been observed at $966\pm 2$,
$986\pm 2$, and $1003\pm 2$ MeV \cite{conf1,conf2,epj,hadr}.
These values of $M_{X_1}$ coincide with
the ones obtained by the simulation and essentially differ from
the value of the neutron mass (939.6 MeV). Hence, for all states under
study, we have $X_1=\gamma+n$ in support of the statement that the
dibaryons found are SNDs.

On the other hand, the peak at $M_{X_1}=1003\pm 2$ MeV corresponds to
the peak found at SPES3 (Saturne) ref. \cite{tat2} and was attributed to an
exotic baryon state $N^*$ below the $\pi N$ threshold. In that work the
authors investigated the reaction $pp\to\pi^+pX$ and have found
altogether three such states with masses 1004, 1044, and 1094 MeV
with SD$\gtrsim 10$. In additional, states with masses close to 966 and
986 MeV were also extracted at SPES3, but  from a small  number of
data \cite{tat3}.

Therefore, if the exotic baryons with anomalously small masses really
exist, the observed peaks at 966, 986, and 1003 MeV might be a manifestation
of such states. The existence of such exotic states, if
proved to be true, will fundamentally change our understanding of the
quark structure of hadrons \cite{bald,walch,azim}.

However, in experiments on a single nucleon, no any significant structure
was observed \cite{lvov,jiang,kohl}. Therefore,
the question about a nature of the peaks observed in \cite{epj,tat2}
remains open at present.

In ref. \cite{khr} dibaryons with exotic quantum numbers were searched for
in the process $pp\to pp\gamma\gamma$. The experiment was performed with
a proton beam from the JINR Phasotron at an energy of about 216 MeV. The
energy spectrum of the photons emitted at $90^{\circ}$ was measured.
As a result, two peaks have been observed in this spectrum. This behavior
of the photon energy spectrum was interpreted as a signature of the exotic
dibaryon resonance $d_1^*$.
with a mass of about $1956\pm 2 stat\pm 6 syst$ MeV and possible isospin
$T=2$ or $T=1$.

On the other hand, an analysis \cite{cal} of the Uppsala proton-proton
bremsstrahlung data looking for the presence of a dibaryon in the mass
range from 1900 to 1960 MeV gave only the upper limits of 10 and 3 nb for
the dibaryon production cross section at proton beam energies of 200 and
310 MeV, respectively.

However, it is asserted in ref. \cite{khr2} that a more detailed analysis
of the data of ref \cite{cal} gives a confirmation of the existence of
$d_1^*(1956)$.

\begin{figure}[ht]
\epsfxsize=6cm
\epsfysize=4cm
\centering{
\epsffile{fig5.eps}}
\caption{The missing mass $(MM(\g,\pi^0)-m_d)$ (MeV)
spectra of the reaction $\go$ }
\label{mami}
\end{figure}
The reactions $pd\to ppX$ and $pd\to pdX$ have been studied also
in the Research Center for Nu\-clear Physics at the pro\-ton energy 295 MeV
over a mass range of 1898 to 1953 MeV \cite{tamii}. They did not observe
any narrow structure in the missing mass spectra of $pX$, $dX$, and $X$.

So, these results are at variance both with the SND observation in INR and
with the results of Tatischeff {\em et al.} \cite{tat2,tat3} on the search for
exotic baryons. However, exotic baryons were observed in Ref. \cite{tat2}
with a sufficiently high accuracy, which leads to a doubt about correctness
of the result of RCNP \cite{tamii}.

It is worth noting that the reaction $pd\to NX$ was investigated in
other works, too. However, in contrast to
the ref. \cite{konob,izv,yad,prc,conf1,conf2,epj,hadr},
the authors of these works did not study
either the correlation between the parameters of the scattered proton and
the second detected particle
or the emission of the photon from the dibaryon decay.
Therefore, in these works the relative contributions of the dibaryons under
consideration were small, which hampered their observation.

On the other hand, the preliminary analysis of the missing mass distributions
of the available data on the process $\go$, obtained at MAMI \cite{kash},
demonstrates three peaks (Fig. \ref{mami}), which good enough correspond to
the values of the SND masses found in INR \cite{prc,epj} (red lines).
Unfortunately, the statistic is very poor in this case.

\begin{figure}[ht]
\epsfxsize=6cm
\epsfysize=5cm
\centering{
\epsffile{fig6.eps}}
\caption{The re\-la\-tio of $\sigma (\alpha=90^o)/\sigma (\alpha=0^o)$
as the function of $P$.}
\label{polar}
\end{figure}
In order to argue more convincingly that the states found are
really SNDs, an additional experimental investigation of the dibaryon
production is needed.

In ref. \cite{alek1} a search for SNDs in the processes of pions
photoproduction by the linearly polarized photon was proposed.
The cross section of this process can be written as the following:
\beq
 \frac{d\sigma}{d\Omega}=A+\frac{q^2}{2}\sin^2\theta_{\pi} \, B
(1-P \,\cos 2\alpha)
\eeq
where $\alpha$ is the angle of the photon polarization relative to
the reaction plane.

The result of calculations of $\sigma (\alpha=90^o)/\sigma (\alpha=0^o)$
for the process $\overrightarrow{\g} d\to \pi^+ D$ at $M=1904$ MeV is shown
on Fig. \ref{polar} as the function of the polarization degree $P$.

So, we have a large variation between the differential cross sections of
pions propagating parallel to the initial photon polarization and pions
propagating perpendicular to the polarization. It is expected that the
cross section for vector SNDs should be substantially larger for mesons
emitted parallel to the photon polarization than for mesons emitted
perpendicularly. Thus, in the first case, the sensitivity to the contribution
of $D(1,1^-)$ significantly increased.

\begin{figure}[ht]
\epsfxsize=6cm
\epsfysize=10cm
\centering{
\epsffile{fig7.eps}}
\caption{The missing mass spectra of the reactions $pd\to p+ pX$,
and $\vec{\gamma} d\to \pi^+nn\gamma$ when $\pi^+//\vec{\gamma}$
and $\pi^+\perp\vec{\gamma}$.}
\label{norum}
\end{figure}

The calculations showed that $D(1,1^-)$ gives the main contribution
to the amplitude $B$, and its contribution to $A$ is less than 1\%.
The contributions of $D(1,1^+)$ in $A$ and $B$ are almost equal.
If we consider the SND masses up to 1960 MeV, then it is
expected that when a pion is emitted parallel to the photon
polarization, 3 peaks should be observed and only 1, when the pion
is emitted perpendiculary.

Such an experiment has been performed at LEGS \cite{legs,norum}. They
analyzed the reaction $d(\overrightarrow{\g},\pp n\g')n$ in the photon
energy range 210 -- 340 MeV. The linear polarization of initial photon
was approximately 99\%. The results of this experiment and their comparison
with the data obtained at INR \cite{prc,epj} are shown in Fig. 7.

As a result they have observed three peaks in missing mass spectrum
when the $\pp$ was emitted parallel to the polarization of the incident
photon $\overrightarrow{\g}$.

The mass values found are very close
to the values obtained in INR \cite{prc,epj}

The peak at $M=1926$ MeV in the bottom panel corresponds to the expected
value. However, this peak is too wide, possibly due to insufficient
accuracy in the determination of experimental data in the region of this
resonance.

The results obtained in this experiment,
support the possibility of the SND existence. However, these data were
limited by the resolution of the pion detection. So, they did not
produce a conclusive proof of the SND existence.

It should be noted that SNDs could be produced in the processes under
consideration, if a pion is only emitted from the 6-quark state of the
deuteron. Therefore the vertexes of $d\to \pi D$ can be written as
\beq
&&\a_{d\to\pi D(1,1^{-},0)}=\frac{g_1}{M} \sqrt{\eta}
\Phi_{\mu \nu} G^{\mu \nu} ,\\
&&\a_{d\to\pi D(1,1^+,1)}=\frac{g_2}{M}\sqrt{\eta}\varepsilon_{\mu\nu
\lambda\sigma}\Phi^{\mu\nu}G^{\lambda\sigma},\\ \nonumber
\eeq
where $\Phi_{\mu \nu}=r_{\mu}w_{\nu}-w_{\mu}r_{\nu}$,
$G_{\mu \nu}=p_{1\mu}v_{\nu}-v_{\mu}p_{1\nu}$, $w$ and $v$ are 4-vectors
of the dibaryon and deuteron polarizations, respectively; and
$r$ and $p_1$ are the dibaryon and deuteron 4-momenta.

The constants $g^2_1/4\pi$, $g^2_2/4\pi$, and $\eta$ are unknown.
However, the products of these coupling constants and $\eta$ can be
estimated from the results of work \cite{legs} where the SNDs,
were searched for in the process $\overrightarrow{\g} d \to \pp D\to\g'\pp nn$.

 As a result, we have
\begin{equation}
\eta\frac{g_1^2}{4\pi}=1.4\times 10^{-4}, \qquad
\eta\frac{g_2^2}{4\pi}=3\times 10^{-4}.
\end{equation}

\section{Mass formula for the SNDs}

Using the complete Green function of the dibaryons
$$
\Delta(p^2)=\frac{F(p)}{p^2-m^2-\delta_D(p^2)},
$$
we determine the SND mass as
\begin{equation}\label{mdb}
M^2=m^2+Re\,\delta_D(M^2),
\end{equation}
where $\delta_D(M^2)$ is the self energy of the SND under study and
$m$ is the mass of the dibaryon in the intermediate state.

The self energy of the lightest SND will be determined in one loop
approximation through the interaction of the pion with the  6-quark
state of the deuteron.
The self energy of the next SND will be obtain through the interaction of
the pion with the lightest SND and so on.

We calculate the SND self energy with the help of
the dispersion relations with two subtractions at $M^2=m^2$.
Then taking into account (\ref{mdb}), we obtain the mass formula for the
SNDs \cite{mass}
\begin{eqnarray}
\label{dsd}
M^2&=&m^2+Re\,\delta_D(m^2)+\left.(M^2-m^2)\frac{
d\,Re\,\delta_D(M^2)}{d\,M^2}\right|_{M^2=m^2}+ \nonumber \\
&&\frac{(M^2-m^2)^2}{\pi}P\!\!\!\int\limits_{(m+\mu)^2}^{\infty}
\frac{Im\,\delta_D(x)\,dx}{(x-M^2)(x-m^2)^2}.
\end{eqnarray}
Since the subtraction is carried out on the mass shell of the
dibaryon in the intermediate state, the subtraction function
$Re\,\delta_D(m^2)$ is equal to zero. Assuming that this dibaryon is
in the ground state, we have
$\left.d\,Re\,\delta_D(M^2)/d\,M^2\right|_{M^2=m^2}=0$.

Finally, the mass formula for SND can be represented as
\begin{equation}
\label{dsd1}
\frac{(M^2-m^2)}{\pi}P\int\limits_{(m+\mu)^2}^{\infty}
\frac{Im\,\delta_D(x)\,dx}{(x-M^2)(x-m^2)^2}=1   \nonumber
\end{equation}

Two subtractions in the dispersion relations provide a very good
convergence of the integrand in (\ref{dsd1}). Therefore we restrict
ourselves to consideration of one loop approximation only.

We assume that the SND under study and the dibaryon in the intermediate
state have opposite parities.
Then the vertex $D^{\prime}(1^{\mp})\to\pi+D(1^{\pm})$
can be written as
\begin{equation}
\Gamma^{(-)}=\frac{\bar{g}_1}{M}G_{\mu\nu}\Phi^{\mu\nu}. \qquad
\end{equation}

As a result of  calculations we have got the following expression for
the imaginary part of $\delta_D(x)$:
\begin{equation}\label{oppos}
Im\,\delta_D(x)=\left(\frac{\bar{g}^2_1}{4\pi}\right)
\frac{q[(x+m^2-\mu^2)^2+2m^2x]}{x^{\frac32}},
\end{equation}
where $q$ is the pion momentum equal to
$q=[(x-(m+\mu)^2)(x-(m-\mu)^2]^{1/2}/2x^{1/2}$.

The coupling constant $\bar{g}^2_1/4\pi$ in the vertex for transition of
the 6-quark state of the deuteron $(D(0,1^+))$ plus the pion to the
SND $D(1,1^-)$ has been fixed by requiring a reproduction of the mass
$M=1904$ MeV. It results in
\begin{equation}\label{const1}
\frac{\bar{g}^2_1}{4\pi}=26.5888 .
\end{equation}

Calculations within the framework of the present model yielded very close
values of the SND masses found in channels with $\pi^0$ and
$\pi^{\pm}$ mesons . Therefore, we take them equal one to another.

In order to calculate the mass of the next SND $D(1,1^+)$, we take in the
intermediate state the SND $D(1,1^-)$ with $m=1904$ MeV and the pion.
To calculate the mass of the next $D(1,1^-)$, we consider SND
$D(1,1^+)$ with $m=1925$ MeV and pion in the intermediate state and so on.

The coupling constant $\bar{g}_2$ in the vertex
$D(1,1^{\pm})\to D(1,1^{\mp})+\pi$ differs from $\bar{g}_1$.
Due to the isotopic invariance, we have
$\bar{g}^2_2/4\pi=3/4 (\bar{g}^2_1/4\pi)$.

The results of the calculations of the masses
and $J^P$ of the SNDs are listed in table \ref{dibm}.
\begin{table}
\centering
\caption{The masses and $J^P$ of the SNDs.}
\label{dibm}
\begin{tabular}{|c|c|c|c|c|}\hline
     &     &  model  &  experiment & experimental \\
$No$& $J^P$&$M$ (MeV)& $M$ (MeV)   &  works       \\ \hline
 1   &$1^-$& 1904    &$1904\pm 2$  & \cite{epj}   \\ \hline
 2   &$1^+$& 1925    &$1926\pm 2$  & \cite{epj}   \\ \hline
 3   &$1^-$& 1945    &$1942\pm 2$  & \cite{epj}   \\ \hline
 4   &$1^+$& 1965    &$1956\pm 6$  & \cite{khr}   \\ \hline
 5   &$1^-$& 1985    & $1982$ &predicted \cite{epj,tat2}\\ \hline
 6   &$1^+$& 2006    &            &               \\ \hline
\end{tabular}
\end{table}

As can be seen from the table, the values of the SND masses obtained are in
good agreement with available experimental data. The existence of the SND
with the mass $M=1982$ MeV was predicted in \cite{epj,tat2}
as a consequence of the observation of the peak in the missing mass
spectrum of the reaction $pp\to\pi^+pX$ \cite{tat} at $M_X=1044$ MeV.

An analysis of the probability of SND decay shows that due to the
smallness of the angle $\theta_{12}$ and difference in the energies
of the final nucleons, resonance-like states with  masses
$(p_1+k)^2\simeq (p_2+k)^2=M_{N^*}^2$ appear.
The table \ref{nstar} presents the values of $M_{N^*}$ for
various SND masses.
\begin{table}
\centering
\caption{The masses of the SNDs and $N^*$.}
\label{nstar}
\begin{tabular}{|c|c|c|c|c|c|c|}\hline
$J^P$&$1^-$ &$1^+$ & $1^-$ & $1^+$ & $1^-$ & $1^+$ \\ \hline
 M(D)&$1904\pm 1$&$1925\pm 1$&$1945\pm 1$&$1965\pm 1$&$1985\pm 1$
 &$2006\pm 1$\\ \hline
$M(N^*_{SR})$&$965\pm 2$&$986\pm 2$&$1005\pm 2$&$1025\pm 1$&$1044\pm 1$
 &$1063\pm 1$\\ \hline
$M(N^*_{exp})$&$966\pm 2$&$986\pm 2$&$1003\pm 2$& &$1044\pm 2$& \\
      &INR\cite{epj}&INR\cite{epj}&INR\cite{epj}& &SPES3\cite{tat2}& \\
              &  &SPES3\cite{tat3} &SPES3\cite{tat2}  & &           & \\
\hline
\end{tabular}
\end{table}

As the SNDs observed in \cite{prc,epj} decay into $NN^*$ with the small
relative momentum between $N$ and $N^*$, the SND parity has to be
determined by the parity of the $N^*$. As seen from tables II and III,
the SND parities found here agree with the results obtained for
the $N^*$ states.

Thus, the formation of $N^*$ is uniquely determined by the existence of
SNDs. $N^*$ is not a new particle. This is a kinematic effect
associated with the SND decay form. Consequently, the experimental
observation of $N^*$s can be an additional indication of the possibility
of the existence of SNDs.

\section{Summary}

The main properties of SNDs have been considered. A number of experiments
\cite{prc,epj,khr,kash,legs}, where evidence was obtained for the
existence of SNDs, has been analyzed.

Negative results obtained in RCNP \cite{tamii} are at variance both with
the observation of SNDs in INR \cite{prc,epj} and with the result of
B. Tatischeff {\em et al.} \cite{tat2} on search for exotic baryons. On the
other hand, the latter was observed in \cite{tat2} with sufficiently
high accuracy, which leads to doubt about correctness of the result of
RCNP \cite{tamii}.

The sum rules for the SND masses have been constructed. The values  of
the SND masses obtained by means of the sum rules agree very well with
the experimental data \cite{epj,tat2}.

The decay of SNDs leads to the formation of resonance-like states $N^*$.
The predicted values of the masses of $N^*$ are in good agreement with
available experimental data.

The experimental observation of $N^*$s is an additional conformation
of the possibility of the SND existence.

The author thanks V.L. Kashevarov and B. Norum for helpful discussion.

redpajama_set_name: RedPajamaArXiv


text:
\section{}

\maketitle

\section{Introduction}
In Parrondo's game~\cite{harmer99b}, in its usually form as a gambling game,
a single agent plays against a bank, with the choice of two games $A$ and $B$,
whose results are determined by the toss of biased coins.
Each of these games is losing when played in isolation
but when played alternately or in some other deterministic or random sequence
(such as $ABAB$ etc.)
can become a winning game.
Because of its anti-intuitive nature, Parrondo's games are also referred to as Parrondo's paradox.
Game $B$ is constructed in such a way that there is a form of feedback from the game state that,
when the game is played by itself,
causes the game to be a losing game.
Interspersing plays of game $A$ disrupts the feedback and the net result can be a winning game.
The feedback can be in the form of dependence on the total capital of the player~\cite{harmer99a},
on the history of previous wins and losses~\cite{parrondo00,kay02}
or, in a multiplayer setting, on the performance of neighbouring players~\cite{toral02}.

Interest in quantum models of Parrondo's games was initiated by Meyer~\cite{meyer02a,meyer02b},
who presented a model analogous to the capital-dependent Parrondo's game,
and was later extended to the history-dependent model~\cite{flitney02b,flitney04b}.
A close quantum analogue of the capital dependent Parrondo's game
was also developed by Ko\u{s}\'{i}k~\cite{kosik07}.

A quantum analogue of a classical discrete random walk
can be obtained by replacing
the classical coin by an SU(2) operator.
If the walker position on a one-dimensional line represents the player's total capital,
a quantum walk becomes a natural setting for a quantum Parrondo's game.
As distinct from the classical case,
in a quantum setting the change from winning to losing
may be achieved by interference effects between the two games.

In the quantum Parrondo's game model of Ref.~\cite{chand11}
two players A and B play alternately.
Winning for player A is defined as motion to the right
while motion to the left is defined as winning for player B.
It is then obvious that when
A and B use their coins alternately that one or the other of them will `win',
or that the expectation of the walker position will be zero,
in which case the players are joint winners, according to the authors' criteria.
This does not have the anti-intuitive aspect that is characteristic of Parrondo's
games nor is it a single player with a choice of two games.
A two-dimensional analogue of this idea has also appeared as an e-print~\cite{ampadu11}.
Despite its failure to meet the criterion as a Parrondo's game the idea
used in Refs.~\cite{chand11,ampadu11}
of biasing a quantum walk using phase factors
can be used to construct a version of a quantum Parrondo's game that relies on quantum interference
to produce the counter-intuitive motion of the walker in a one-dimensional quantum walk.

In this paper
we first introduce quantum walks
and discuss how to control them using the coin parameters and the initial state.
This is well covered in the literature
and is included here briefly only for completeness.
In Sec.~\ref{sec:our_model} the new quantum Parrondo's model is presented
and Sec.~\ref{sec:results} presents results for all combinations of games $A$ and $B$ of length four or less,
giving the regions in parameter space for which the effect can be achieved.
It will be noticed that although the expectation of the walker position
is reversed for some finite number of time steps,
eventually the original bias reasserts itself.
Thus we term this phenomena {\em transient} quantum Parrondo's games.

\section{Basics of quantum walks}
\label{sec:qwalks}
The usefulness of classical random walks in mathematics, theoretical computer science
and computational physics~\cite{papa94,motwani95}
has inspired significant interest in quantum walks,
both in continuous-time~\cite{farhi98,childs02}
and in discrete-time~\cite{aharonov93,meyer96,aharonov01,ambainis01}.
Quantum walks reveal a number of startling differences to their
classical counterparts.
In particular, diffusion on a line is quadratically faster \cite{nayak00,travaglione02},
while propagation across some graphs is exponentially faster \cite{childs02,childs03}.
Quantum walks also show promise as a means of implementing some quantum algorithms~\cite{shenvi03,ambainis04}.
An overview of quantum walks is given in Ref.~\cite{kempe03a}
and a recent comprehensive review appears in Ref.~\cite{va12}.

If a quantum particle moving in discrete time steps along a line
is updated at each step, in superposition, to the left and right,
the global process is necessarily non-unitary~\cite{meyer96}.
However, an additional degree of freedom,
the chirality, taking values $L$ and $R$,
allows unitary quantum walks to be constructed.
Let ${\cal H}_{\rm P}$ be the Hilbert space of discrete particle positions in one-dimension,
spanned by the basis $\{|x\rangle : x \in \mathbb{Z}\}$.
Let ${\cal H}_{\rm C}$ be the Hilbert space of chirality, or ``coin'' states,
spanned by the orthonormal basis $\{|L\rangle, |R\rangle\}$.
A discrete-time, discrete-space quantum walk
on the Hilbert space ${\cal H}_{\rm P} \otimes {\cal H}_{\rm C}$
consists of a unitary operation $\hat{U}$ on the coin state,
followed by the updating of the position to the left or right
according to the result, that is,
\begin{equation}
\hat{E} = (\hat{S} \otimes \hat{\cal P}_{R} \;+\; \hat{S}^{-1} \otimes \hat{\cal P}_{L})
        (\hat{I}_{\rm P} \otimes \hat{U}),
\end{equation}
where $\hat{S}$ is the shift operator in position space,
$\hat{S} |x\rangle = |x+1\rangle$,
$\hat{I}_{\rm P}$ is the identity operator in position space,
and $\hat{\cal P}_R$ and $\hat{\cal P}_L$ are projection operators on the coin space
with $\hat{\cal P}_R + \hat{\cal P}_L = \hat{I}_{\rm C}$, the coin identity operator.
For example, a walk controlled by an unbiased quantum coin
can be carried out by the transformations
\begin{equation}
\label{eq:step1}
\begin{split}
|x, L \rangle &\rightarrow \, \frac{1}{\sqrt{2}}
        \left( |x-1, L \rangle \,+\, i |x+1, R \rangle \right), \\
|x, R \rangle &\rightarrow \, \frac{1}{\sqrt{2}}
        \left(i |x-1, L \rangle \,+\, |x+1, R \rangle \right).
\end{split}
\end{equation}
To achieve an unbiased quantum walk with Eq.~(\ref{eq:step1})
we require an initial state
$|\psi_0\rangle = (|0, L \rangle - |0, R \rangle)/\sqrt{2}~$.
The scheme is then equivalent to the Hadamard quantum walk
with initial state $(|0, L\rangle + i |0, R\rangle)/\sqrt{2}$.

Tregenna {\em et al.}~\cite{tregenna03} has shown that with a general initial state
\begin{equation}
\ket{\psi(x,0)} =
	\left( \sqrt{\eta} \ket{R} + e^{i \mu} \sqrt{1-\eta} \ket{L} \right) \otimes \ket{0},
\end{equation}
and with a coin operator of maximum generality (up to an overall phase)
\begin{equation}
\label{eq:general_flip}
\hat{U}(\rho, \theta, \phi) = \begin{pmatrix}
											\sqrt{\rho} 				 & \;e^{i \theta} \sqrt{1-\rho} \\
											e^{i \phi} \sqrt{1-\rho} & \;-e^{i (\theta + \phi)} \sqrt{\rho}
			 						   \end{pmatrix},
\end{equation}
$\phi$ disappears in the final dynamics,
while $\theta$ only appears in the combination $\theta + \mu$.
Both $\rho$ and $\eta$ have non-trivial effects on the dynamics
but an unbiased walk requires $\rho = \frac{1}{2}$.

\section{Parrondo's game constructed from alternating coin phases}
\label{sec:our_model}
In the following we will chose the initial state
$|\psi_0\rangle = (|0, L \rangle - |0, R \rangle)/\sqrt{2}~$,
and a more symmetrical form of the coin operator,
\begin{equation}
\label{eq:flip}
\hat{U}(\rho, \alpha) = \begin{pmatrix} e^{i \alpha} \sqrt{\rho} & i \sqrt{1-\rho} \\
                         i \sqrt{1-\rho} & e^{-i \alpha} \sqrt{\rho}
                \end{pmatrix},
\end{equation}
obtained by multiplying the operator in Eq.~(\ref{eq:general_flip})
by an overall phase $e^{i \alpha}$,
and setting $\theta = \phi = \frac{\pi}{2} - \alpha$.
With $\rho = \frac{1}{2}$ and $\alpha=0$ an unbiased quantum walk is obtained.
However, with $\rho = \frac{1}{2}$ a bias can be introduced,
to the left or right,
by the use of a non-zero phase factor $\alpha$.
If we restrict $\alpha \in [-\frac{\pi}{2}, \frac{\pi}{2}]$,
then $\alpha > 0$ causes average motion to the left while $\alpha < 0$ causes motion to the right.
Introduce the operators
\begin{equation}
\hat{A}(\alpha) = \hat{U}(\frac{1}{2}, \alpha), \;
\hat{B}(\alpha) = \hat{U}(\frac{1}{2}, \beta),
\end{equation}
for some $\alpha, \beta \in [\frac{-\pi}{2}, \frac{\pi}{2}]$.
Connection with Parrondo's games can now be obtained by choosing to use $\hat{A}$ and $\hat{B}$
for games $A$ and $B$, respectively,
with both $\alpha, \beta > 0$.
Thus, quantum walks generated by either $\hat{A}$ or $\hat{B}$ alone
are biased towards the negative (`losing').
The loss is an approximately linear function of the number of time steps $t$.
For example,
$\langle x \rangle |_{t=50} \, \simeq -15.0 \, \sin \alpha$
and $\langle x \rangle |_{t=100} \, \simeq -29.6 \, \sin \alpha$.
However, for some choices of $\alpha$ and $\beta$,
and some combinations of $A$ and $B$,
interference effects can reverse the dynamics
giving a positive expectation value of the walker position (`winning'),
at least up to some finite number of time steps.
Alternately playing $A$ then $B$ $n$-times from the initial state $|\psi_i \rangle$
produces a final state $|\psi_f \rangle$ given by
\begin{equation}
|\psi_f \rangle = ( \hat{B} \otimes \hat{A} )^n \, |\psi_i \rangle.
\end{equation}
The final states for other sequences are produced analogously.

\section{Results}
\label{sec:results}
All repeated sequences of four games were tested.
Of these, only $AB$, $ABB$ and $ABBB$ were found to give some positive expectation for the walker
position after 100 time steps.
In each case there is only a small region of values of the phase parameters $\alpha$ and $\beta$
that give rise to the positive positions,
as indicated in Figs~\ref{f:AB100}--\ref{f:ABBB100}.
Note that interchanging games $A$ and $B$ only has the effect of reversing the role of $\alpha$ and $\beta$
and gives no new dynamics.
The results for game $A$ alone can be obtained by setting $\alpha = \beta$.
In each case this results in a more negative value of $\langle x \rangle$
for a given number of time steps $t$
than the alternating sequence for near-by values of $\alpha \ne \beta$.

For a sufficiently large number of time steps
all sequences eventually show a return to a negative expectation value for the walker position.
Some examples of the expectation value of the walker for the first 1000 time steps are shown in
Figs.~\ref{f:AB_change_ab}--\ref{f:ABBB_change_ab},
for the sequences $AB$, $ABB$ and $ABBB$, respectively.
The solid (near) straight lines in the right hand figures indicate the result for playing game $A$ only,
for $\alpha = 0.005$.
These graphs indicate that $\langle x \rangle$ has an oscillatory pattern superimposed upon a downward trend.
Decreasing $\alpha$ leads to a lessening of the downward slope,
while decreasing $\beta$ leads to longer periods of oscillation.
A short period, low amplitude oscillation is superimposed upon the longer one.
The amplitude of this oscillation becomes more pronounced as we go from the sequence $AB$ to $ABB$ to $ABBB$,
as close inspection of the figures show.

\section{Conclusion}
We have shown that a quantum walk on a line can give rise to a transient quantum Parrondo's game.
Two different biased coin operators $\hat{A}$ and $\hat{B}$
play the role of the two games in the Parrondo effect.
Bias is achieved by adding a phase factor to the coin operator.
Although each operator,
when used alone,
gives negative expectation values for the walker position as the number of steps increases,
when used in the repeated sequences $AB$, $ABB$ or $ABBB$
can give rise to movement in the positive direction,
for some values of the phase factors.
We have detailed the range of values of these phase factors that give rise
to the Parrondo effect after 100 time steps
for each of these sequences.
The effect is transitory.
For large enough times all sequences eventually result in a negative
expectation value for the  walker position.