Oscillation tomography of the Earth with solar neutrinos and future experiments

We study in details the Earth matter effects on the boron neutrinos from the Sun using recently developed 3D models of the Earth. The models have a number of new features of the density profiles, in particular, substantial deviation from spherical symmetry. In this connection we further elaborate on relevant aspects of oscillations ($\epsilon^2$ corrections, adiabaticity violation, entanglement, {\it etc.}) and the attenuation effect. The night excesses of the $\nu e-$ and $\nu N-$ events and the Day-Night asymmetries, $A_{ND}$, are presented in terms of the matter potential and the generalized energy resolution functions. The energy dependences of the cross-section and the flux improve the resolution, and consequently, sensitivity to remote structures of the profiles. The nadir angle ($\eta$) dependences of $A_{ND}$ are computed for future detectors DUNE, THEIA, Hyper-Kamiokande and MICA at the South pole. Perspectives of the oscillation tomography of the Earth with the boron neutrinos are discussed. Next generation of detectors will establish the integrated day-hight asymmetry with high confidence level. They can give some indications of the $\eta-$ dependence of the effect, but will discriminate among different models at most at the $(1 - 2)\sigma$ level. For the high level discrimination, the MICA-scale experiments are needed. MICA can detect the ice-soil borders and perform unique tomography of Antarctica.


I. INTRODUCTION
Oscillations of the solar neutrinos in the Earth [1] - [32] have the following features.
1. Due to loss of the propagation coherence, the solar neutrinos arrive at the surface of the Earth as independent fluxes of the mass eigenstates [8], [9,15,24].
2. Inside the Earth the mass states oscillate in multilayer medium with smoothly (adiabatically) changing density within layers and sharp density change at the borders between the layers.
3. The oscillations proceed in the low density regime which is quantified by small parameter where V (x) = √ 2G F n e (x) is the matter potential, n e is the electron number density of the medium. For E = 10 MeV at the surface of the Earth equals ∼ 0.03. 4. The oscillation length l m ≈ l ν ≈ 330 km E 10 MeV 7.5 · 10 −5 eV 2 ∆m 2 21 is comparable to a section of trajectory in a layer d i for trajectories with nadir angles η close to π/2: d i = r i / cos η, where r i ∼ 10 km is the width of the layer in the radial direction. The highest sensitivity is to structures of the density profile of the size ∼ l m /2.

5.
The attenuation effect is realized in the order due to the finite neutrino energy resolution (reconstruction) in the experimental setup [26,32]. It means loss of sensitivity to remote structures of the Earth density profile. Consequently, only structures sufficiently close to a detector, and therefore to the surface of the Earth (crust, upper mantle), are mostly relevant for observations. This means that with the boron neutrinos, deep structures, like the core of the Earth, are not seen at the level. The attenuation effect is absent in the order 2 . Thus, the solar neutrino tomography is essentially sensitive to the small scale structures in the crust and mantle of the Earth.
In previous computations, (see, e.g., [24], [20]) the density profile of the one-dimensional PREM model [33] was used. In this model, borders between layers have forms of ideal spheres. Recently several new three dimensional Earth models have been developed. They show several new features of the density profiles which have not been taken into account previously: (i) the borders between layers are not spherically symmetric but have irregular deviations from spheres; (ii) the profiles depend on the azimuthal angle; (iii) The profiles are non-symmetric with respect to the center of neutrino trajectory. The horizontal sizes of these structures are comparable to oscillation length which means that effectively they can smooth borders between layers as well as produce some new parametric effects in oscillations.
In the present paper, we study how these new features modify the observational effects. We compute the Earth arXiv:2001.08030v1 [hep-ph] 22 Jan 2020 matter effect using new models. This allows us to assess a possibility do distinguish the models with solar neutrino detectors. At the same time, our computations quantify errors of the computed effects due to uncertainty in the density profile.
Presently there are first (about 3σ) indication of the Earth matter effect by SuperKamiokande [29], and this situation will stay till the next generation of experiments will start to operate. Here we consider solar neutrino studies by future detectors DUNE [34], Hyper-Kamiokande (HK) [35], THEIA [36] [37] and MICA [38].
The paper is organized as follows. In Sec. II, we present oscillation formalism relevant for our computations and elaborate on some new features, such as high order corrections, entanglement, etc. We introduce the generalized energy resolution functions and study their properties. The Day-Night asymmetry is presented in terms of these resolution function and potential. In Sec. III, new models of the density distribution in the Earth are described. In Sec. IV, we present results of computations of the Earth matter effect for future detectors. Conclusions are given in Sec. V.

II. RELATIVE EXCESS OF THE NIGHT EVENTS AND ATTENUATION Coherence and entanglement
Loss of the propagation coherence is due to spatial separation of the wave packets that correspond to the mass eigenstates originated from the same flavor state. Although separated, these wave packets belong to the same wave function and therefore entangled. If one of the eigenstates is detected the parts of the wave function, which describe two other eigenstates, collapse. It can be easily shown that observational result is the same as in the case of independent fluxes of mass eigenstates once total flux of these states is normalized on the total flux of the originally produced flavor neutrinos. Coherence is not restored in a realistic detector.

Corrections to probability
Recall that the ν e survival probability during a day, as function of the neutrino energy, equals P D (E) = 1 2 c 4 13 1 + cos 2θ 12 cos 2θ m 12 (E) + s 4 13 , where c 13 ≡ cos θ 13 , s 13 ≡ sin θ 13 , and cos 2θ m 12 is the mixing parameter averaged over the boron neutrino production region in the Sun [40]: Here¯ andV is the averaged matter potential in the 8 B neutrino production region. For high energy part of the boron neutrino spectrum, where¯ 1, we have So, dependence on E is weak. At the solar neutrino energies the matter effect on the 1-3 mixing is negligible, During a night the probability equals P N = P D + ∆P , where the difference of the night and day probabilities is given to the order 2 by [26,41] is slowly changing function of E, and is a correction of the order 2 , since in (7) each integral over x is of the order . The integration in (7) proceeds along a neutrino trajectory In new models of the Earth apart from the nadir angle η the density and potential profiles depend, also on position of the detector x D and azimuthal angle φ a : V = V (x, x D , η, φ a ). Correspondingly, for a given detector and a given of moment of time ∆P = ∆P (x D , η, φ a ). In Eq. (6) is the adiabatic phase acquired from a given point of trajectory x to a detector at L. ∆ m 21 (x) is the level splitting and in our calculations we use it up to the first order in : Here ∆ 21 ≡ ∆m 2 21 /2E is the splitting in vacuum. Consequently, the oscillation phase (8) equals Introducing the average density along a neutrino trajectoryρ(η), we can rewrite Eq. (9) as where φ m 0 ≡ ∆ 21 (L − x) is the zero order phase and is the phase shift due to the − correction. For ∆m 2 = 7.5×10 −5 eV 2 and ρ = 5 g/cm 3 the relative size of the correction (second term in Eq. (10)) is about 3%. For large φ m 0 the phase shift δφ m can be observable. E.g., if φ m 0 = 5π, we find δφ m = 27 • . The correction δφ m leads to the shift of oscillatory pattern in the η scale. Since δφ m = ∆ 21 δL(η) and L = 2R cos η we obtain δη = δφ m 2R sin η∆ 21 .
For η = 70 • we obtain δη = 0.2 • , while period of oscillatory dependence in the η scale for this η equals 2.8 • , i.e., the shift is by 1/14 of the period. δη increases with decrease of η. Let us consider I 2 -the second term in (6). For constant density it can be computed explicitly Apart from¯ 2 , this term contains additional small factor 0.5 cos 2θ 12 ≈ 1/6. As a result, I 2 is about 0.015% and therefore can be neglected. Our computational relative errors are of the order of 0.1%. Thus, the largest correction to the probability follows from φ m .

Comments on adiabaticity
In the lowest order in , the sensitivity to structures of the Earth matter profile, its deviation from constant density, appears due to borders between layers which strongly (maximally) break adiabaticity. Indeed, in the adiabatic case the oscillation probability would depend on density at the surface of the Earth and on the oscillation phase. However, in the lowest (zero) order in the phase coincides with the vacuum phase. The matter correction to the phase is proportional to which then appears as 2 in the probability. So, in the adiabatic case there is no sensitivity to the profile in the order.
In general, deviations of borders between layers from spherical form may produce effective smearing of borders for neutrino trajectories with large η, and consequently, to decrease of the adiabaticity violation. That would lead to partial loss of sensitivity to the density profile.
If deviation from spherical form in radial direction, ∆h, and in horizontal direction, l f , are such that neutrino trajectory at certain η crosses the border between the same layers several (many) times, the density gradient along the trajectory will decrease. For density jump in a border ∆ρ the gradient equals ∆ρ cos η/∆h. The scale of density change should be compared with the oscillation length in the adiabaticity condition. As we will see, typical scale of deviation of, e.g., the border between the crust and mantle from spherical form is ∆h ∼ 5−10 km and the horizontal size of the structures is l f ∼ (70−150) km. This gives the slope of the structure η f ∼ ∆h/l f ∼ (2 − 7) • . Therefore double crossing can occur for the trajectories with η > 83 • . For parameters of new Earth models, however, adiabaticity is still strongly broken and multiple crossing of borders can occur only in very narrow intervals of η.
In the lowest order, the result for ∆P (E) in (6) can be reproduced as a result of interference of the "oscillation waves" emitted from borders between layers [41]. For ith wave, the phase is determined by distance from border to a detector L − x i and vacuum oscillation length, while the amplitude is proportional by the density jump ∆ρ i in the border. Then ∆P (E) is the sum of the waves over borders which neutrino trajectory crosses. This representation gives simple interpretation of results of numerical computations.

Attenuation and generalized energy resolution functions
The Earth matter effect can be quantified by the Day-Night asymmetry or the relative excess of night events (events rate) in energy range ∆E as function of the nadir angle η: Here N N (η) and N D are the numbers of night and day events (rates) correspondingly. The nadir η and azimuthal φ a angles are fixed by the detection time of an event. According to new models N N (η) depends also on the position of a detector. In experiments, the observables are the electron energy and direction. Therefore, ∆E is determined by the observed energy interval of the produced (or recoil) electrons. In practice, we will use the energy of electrons above certain threshold. Thus, information on the density profile is encoded in the nadir angle dependence of the night excess. We will not consider the direction of electron.
Sensitivity of oscillations to the Earth density profile is determined by the sensitivity of a given experimental set-up to the true energy of neutrino E. This can be described by the generalized energy resolution function G ν (E r , E) such that where E r is the observed (reconstructed) neutrino energy or certain energy characteristic which can be measured in experiment. In (16) D is the factor which includes characteristics of detection: fiducial volume, exposure time, etc. It cancels in the expression for the relative excess A N D . The resolution function is normalized as G ν (E r , E)dE = 1. Similarly, one can write expression for N D .
G ν (E r , E) includes the neutrino energy resolution function: g ν (E r , E), the energy dependence of the neutrino flux f B (E) [40] and cross-section σ(E): It should also include the energy dependent efficiency of detection.
Integration over the neutrino energy with the resolution function in Eq. (16) leads to the attenuation effect [24,32]. Pluging expression for ∆P (E) from (6) into (16) and neglecting I 2 we obtain for ∆N Here integrations over x and E are interchanged. In this form the dependence of difference of events on structures of density profile is immediate.
Let us introduce the attenuation factor F (L − x) [24] such that the integral over E in Eq. (18) equals In general, this equality can not be satisfied, but it is valid for special cases and under integral over x. Then the expression for ∆N in (18) becomes For the Gaussian form of G ν (E r , E), the attenuation factor is given by where is the attenuation length, and l ν is the oscillation length in vacuum According to (20) and (21) for d λ att the attenuation factor F (d) ≈ 0, and therefore contributions of remote structures to the integral (20) and therefore to observable oscillation effect is suppressed. For d = λ att the factor F (d) = e −2 ≈ 0.14, and the attenuation becomes significant. Consequently, the Day-night asymmetry depends mainly on the shallow structures of the Earth which are close to a detector.
The attenuation length is the distance at which oscillations integrated over the energy resolution interval σ E are averaged out, or the difference of the oscillation phases for E and E + σ E becomes larger than 2π [32].
Expression (18) factorizes different dependences: The generalized resolution function encodes external characteristics: neutrino flux, cross-section, energy resolution of a detector. V (x) gives information about the density profile, oscillation probability is reduced to sin φ m .
In what follows we will find expressions for the generalized reconstruction functions and present numbers of events in the form (18) separately for the ν−nucleon and ν − e scattering.

Neutrino-nuclei scattering
We consider the charged current neutrino-nuclei interactions and the corresponding resolution function G νN . If transitions to excited states are neglected, the energies of electron and neutrino are uniquely related (upto negligible nuclei recoil): E e = E −∆E. Here ∆E ≈ ∆M +m e is the threshold of reaction. If transitions to excited states are significant but the energy of de-excitation is not measured, an additional uncertainty in reconstruction of the neutrino energy appears which should be included into G νN .
The night-day difference of numbers of events with the observed energy of electron E r e is given by e , E e ) is the electron energy resolution function with E e and E r e being the true and the observed energies correspondingly.
Introducing also E r ≡ E r e + ∆E and changing integration in (24) to integration over the neutrino energy E we . The equation (25) can be rewritten as and z being the normalization factor. Inserting expression for ∆P (E) from (6) into (26) and permuting integrations over x and E we obtain Integration over the energy can be removed introducing of the attenuation factor, as in (19), which gives Finally, integration over the interval of observed energies of electrons gives where we again substituted integration over E e by integration over E.
For the day signal, which does not depend practically on η, we have Notice that if threshold ∆E is low enough, the second integral over the resolution function is ≈ 1, so that The factors Dz cancel in the expression for A N D . Let us consider G νN (E r , E) defined in Eq. (27) in details. We take the Gaussian form for g ν (E r , E) with central energy E c = E r and σ E = 0.1E r . Fig. 1 (upper panel) shows dependence of G νN on energy for several values of E r computed with σ ∝ Ep. We compare this dependence with Gaussian form G νN = g Gauss (E r , E). The flux of Boron neutrinos, f B (E), is taken from [40]. For convenience of comparison we normalized G νN (E r , E) (dashed lines) such that G νN (E r , E) max = g Gauss (E r , E) max , and the y-axis is in arbitrary unit.
Dependences of the neutrino energy resolution

solid lines) and
Gν (E r , E) = g Gauss (E r , E) (dashed lines) on true neutrino energy E for different values of the reconstructed neutrino energy E r . For gν (E r , E) we take the Gaussian distribution with width σE = 0.1E r . Upper panel: GνN -for experiments based on the νe− nuclei scattering. Lower panel: Gνe -for experiments based on the ν − e scattering with cut Ee > 6.5 MeV.
The product σ(E)f B (E) has the form of a wide asymmetric peak with maximum at ∼ 11 MeV. Consequently for E r = 11 MeV the generalized function G νN is close to the corresponding Gaussian form with energy of maximum E c ≈ E r , while for E r > 11 MeV the factor σ(E)f B (E) shifts G νN to lower energies E c < E r and reduces the width: According to Fig. 1 for E r = 12 MeV the energy of maximum E c = 11.7 MeV and the relative width σ E /E = 0.074 instead of 0.1 in g ν . The change becomes more profound with increase of E r . For E r = 14 MeV we find E c = 13.1 MeV and σ E /E = 0.069. Thus, the energy dependence of σf B leads to better energy resolution and therefore to increase the attenuation length which means the improvement of sensitivity to remote structures.
According to Fig. 1, for Gaussian g ν , the whole effective resolution function G νN can be well approximated by the Gaussian function with appropriately choosen energy of maximum, E c = E c (E r ) = E r , and width σ E = σ(E r ). In fact, a priory the form of g ν (E r , E) is not known and eventually will be determined in experiment. Therefore in our computations we will use the generalized reconstruction function in the Gaussian form: Under integration over the neutrino energy E the difference of results for A N D computed with the Gaussian G νN (33) and G νN with Gaussian g ν is negligible. Using the PREM model we find that the relative difference results for A N D is smaller than 0.3%.

Neutrino-electron scattering
In this case the energies of neutrino and electron are not uniquely related, but correlated via the differential cross-section dσ(E, E e )/dE e . Correspondingly, expression for the effective resolution function in (18) will differ from G νN . [49] The difference of numbers of the night and day events with a given observed energy of electron E r e equals where is the difference of the charged current, dσ e /dE e , and neutral current, dσ µ /dE e , differential cross-sections. Interchanging integrations over E e and E in Eq. (34) we obtain where and The generalized reconstruction function can be introduced similarly to (27): or explicitly, inserting g ν from (37), as . (40) The only difference from (27) is that here in g ν the electron resolution function is integrated with the differential cross-section.
Instead of E r e we can introduce the "observable" neutrino energy E r = E r (E r e ) defined as the energy of maximum of G νe for a given E r e : In terms of G νe (E r e , E) the N-D difference of numbers of events can be presented as As in the νN − case, we insert explicit expression for ∆P (E) and interchange integration over x and E. Then the integration over E can be removed introducing the attenuation factor which gives The difference of numbers of events with the observable energy of electrons in the interval ∆E r e ≡ (E r,min (44) and E r = E r (E r e ) is determined by (41). The number (rate) of events with the observed electron energy E r e during a day equals Here (46) The total cross-sections are given by Expression (45) can be simplified assuming g µ ν ≈ g e ν ≈ g ν : Let us consider G νe (E, E r e ) in detail. In the bottom panel of Fig. 1 we show G νe (E, E r e ) as function of E computed according to Eq. (40). For the ν − e scattering the product σ(E)f B (E) has wide peak with maximum at E = 10 MeV, and additional weak E -dependence comes from the integral in (40). Therefore the smallest deviation of G νe (E r , E) from the Gaussian form is at E r ∼ 10 MeV. For E r < 10 MeV the maximum of G νe is shifted to higher energies, while E r < 10 MeV -to lower energies. In both cases the width of G νe decreases. According to Fig. 1 (bottom) for E r = 8 MeV the maximum of G νe is shifted with respect to E r to higher energy by 0.3 MeV, and the width is slightly smaller. For E r = 12 MeV, inversely, the maximum is shifted to E c = 11.7 MeV, and the width becomes σ E /E = 0.085. This trend (due to fast decrease of the flux with energy above 10 -11 MeV) is even more significant for larger E r : at E r = 14 MeV, we find E c = 12.6 MeV and σ E /E = 0.08. Again, taking into account the energy dependence of σ and f B improves the energy resolution, but weaker than in the νN case.
The biggest contribution to oscillation effect comes from the energy range (10 -12) MeV, where G νe is rather close to the Gaussian form. Therefore in computations we will use the Gaussian form for G νe , and consequently, the attenuation factor in the form (21). Inclusion of the flux and cross-section energy dependences narrows the resolution function.
In expressions for ∆N the φ a dependence appears in two places: in the potential: V = V (x, η, φ a ) and in the phase φ m = φ m (φ a ). For each η and position of the detector we performed averaging of ∆N over the azimuthal angle φ a . If φ a dependence of the phase is neglected, in the first approximation, the averaging of ∆N over φ a is reduced to averaging of the potential.

III. MODELS OF THE EARTH AND DENSITY PROFILES
In computations, we used density profiles reconstructed from recently developed 3D models of the Earth.
Due to the attenuation effect, the Day-Night asymmetry mainly depends on shallow density structures: crust, upper mantle and crust-mantle border called Moho, or Mohorovicic discontinuity. There are two types of crust: the oceanic crust and the continental one. The width of oceanic crust is about (5 -10) km, while the continental crust is thicker: (20 -90) km [46,47]. The predicted depth of Moho, h M oho , significantly varies for different models. In contrast, the density change in the Moho is nearly the same for all the models. Beneath Homestake the jump is from 2.9 gr/cm 3 to 3.3 gr/cm 3 .
A brief description of relevant elements of the models is given below.
1. The Shen-Ritzwoller model (S-R) [42] is based on joint Bayesian Monte Carlo inversion of geophysical data. It gives the density profile of the crust and uppermost mantle beneath the US, in area with latitudes (20 • −50 • ) and longitudes (235 • − 295 • ). In the radial direction it provides the density change from the sea level surface down to the depth of 150 km with h M oho = 52 km beneath the Homestake (see Fig. 2).
2. FWEA18, the Full Waveform Inversion of East Asia model [43], covers the latitudes 10 • − 60 • and longitudes 90 • − 150 • . It gives the density profile from the surface down to 800 km, and h M oho = 33 km beneath Kamioka.
3. SAW642AN [44] is a global (all latitudes and longitudes) radially anisotropic mantle shear velocity model based on a global three-dimensional tomography of the Earth. The model gives the density profile of mantle starting from the depth of Moho, h M oho = 24 km, down to 2900 km. No crust structure is available.
4. CRUST1 [45] is a global 3D model, that presents data with 1×1 degree grid in latitude and longitude at the surface. It gives the density and depth of borders of eight layers of the crust: water, ice, upper sediments, middle sediments, lower sediments, upper crust, middle crust, lower crust. The model predicts the depths of Moho h M oho = 48 km and h M oho = 40 km beneath Homestake and Kamioka respectively and nearly constant density of the upper mantle down to 100 km. It provides also the density distribution above the sea level.
Using these models we reconstructed the density, and consequently V (x), profiles along neutrino trajectories determined by position of detectors, η and φ a . Maximal depths h max down to which the models provide data are h max (S − R) = 150 km, h max (CRUST1) ≈ 80 km, h max (FWEA18) = 800 km, h max (SAW642AN) = 2900 km. Therefore we reconstructed the density profiles using the following prescription: • for the S-R, CRUST1 and FWEA18 models with relatively small h max we take the SAW642AN profile in the range h = h max − 2900 km.
• Below 2900 km for all the models we use the PREM profile. Recall that PREM -the Preliminary reference Earth model is a one-dimensional model that represents the average (over solid angle) density of the Earth as a function of depth. The depth of Moho in the PREM model equals h M oho = 24.4 km.
Due to attenuation effect possible uncertainties related to these compilations of the profiles do not change results significantly even for small nadir angles.
All the models, but CRUST1, give the density below the sea level. In all simulations, except the case of MICA, we consider the surface of Earth as perfect sphere and take zero density above the sea level. Effect of these simplifications is much smaller than sensitivity of all experiments (but MICA) due to restricted statistics. In the case of MICA, we have taken into account the Earth structures above sea level.
In Fig. 2, we present the S-R and CRUST1 density profiles beneath Homestake for fixed latitude 44.35 • . Both models provide data for this place down to 80 km. Shown is the depth of layers with a given density as function of longitude (azimuthal angle). Notice that at the latitude 44 • the 1 • of longitude corresponds to 76 km at the surface. The black curves show Moho depth, where density jumps approximately from 2.9 to 3.3 g/cm 3 .
Few comments are in order. 1. The surfaces of equal density, and in particular, borders between layers deviate from spherical form.
3. There are narrow spikes of large amplitude and wide regions ∼ 10 • , where the depth increases by 30% with respect to average value. 4. Two models give rather similar density distributions: the average depths and lengths are similar. At the same time, variations of S-R and CRUST1 models are not correlated.
In the case of spherical inner structures the nadir angle η c at which neutrino starts to cross a given border between layers with the depth h equals where r E = 6371 km is the radius of the Earth. For η < η c neutrino crosses this border twice. Neutrino The noticeable difference between the S-R (CRUST1) profile and SAW642AN profile appears above the S-R Moho depth h > 52 km. Below S-R Moho all three models give similar results.
According to Fig. 2 there are deviations of Moho from of ideal sphere of two types: 1) Relatively small variations of 2 − 5 • scale which would correspond to (150 -400) km at the DUNE latitude and the size (depth) ±(2 − 5) km.
2) Long (continental) scale variations of size 50 • with depth 20 km such that the smallest depth, h min = 32 km, is close to ocean and the bigger depth h max = 52 km is in the center of continent. This means that the Moho border varies within the shell (we call it Moho shell) restricted by spherical surfaces with depth 32 − 52 km and average depth 42 km.
The length of neutrino trajectory within the Moho shell equals ≈ 2 2r E (h max − h min ) ≈ 710 km which is 2 times bigger than the oscillation length. According to (48) borders of the Moho shell are seen from a detector site at η min = 84.2 • and η max = 82.7 • . So that for η > η min there is no crossings of Moho: in the interval η = (η min − η max ) one may expect multiple crossing of Moho and since horizontal scale of variations of the border is comparable to the oscillation length, parametric effects are expected. However, averaging over azimuthal angle washes out these effects. For η < η max neutrino trajectory crosses the Moho shell twice, and within each crossing it can be more than one crossing of the Moho border. Substantial effect due to Moho crossings is expected at η ∼ 83 • .
Below 83 • neutrinos cross the Moho in all the models. For smaller η the differences in these models become small.
As an example, in Fig. 3, we show the reconstructed density profiles of three models along the neutrino trajectory which ends at Homestake with η = 75 • on September 23. The length of trajectory equals 3295 km. According to Fig. 3  Clearly the profiles are not symmetric. Moreover, the density decreases to the middle of trajectory, especially for Homestake. This is related to thicker crust in the middle of a continent. FIG. 4: The same as in Fig. 3 but for the detector located at Hida.

IV. PREDICTIONS FOR FUTURE EXPERIMENTS
We compute the oscillation probability during a day time, P D (E), according to Eq. (2). The rate of events is found using Eq. (31) for νN − scattering and Eq. (47) for the νe− scattering. The excess of night event rate was computed using expression in (32) for the νN − scattering and the one in (44) for the νe− scattering. These expressions correspond to ∆P with neglected I 2 , while the phase was computed keeping the correction.
We take the generalized resolution functions G νN (E r , E) and G νN (E r , E) in the Gaussian form with certain values of the relative widths, σ E /E. The nadir angle and A N D (η, φ a ) are computed with one minute time intervals during a year. Then we averaged A N D (η, φ a ) over the azimuthal angle φ a .
To check the stability of the results with respect to the time binning, we calculated A N D (η) for several values of η with time intervals 15 seconds. Since the starting time of calculation is arbitrary, A N D was also computed with 30 seconds shift. We find that the results are stable with respect to these variations and their relative changes are less than 0.1%.
We compute numerically the annual exposures for detectors at Homestake, Hida, and MICA as functions of nadir angle with ∆η = 0.1 • (see Fig. 5). The exposure functions for Homestake is in agreement with that in Ref. [41]. The asymmetry averaged over the year is given by integration of A DN with the exposure (weight) function W (η) over η: We used exposure functions to compute the expected experimental errors for different η− intervals. The value ∆m 2 21 = 7.5×10 −5 eV 2 is used unless specially indicated.
For this process we use a generic form of cross-section where A is a factor irrelevant for the relative excess, p e is the momentum and E e = E ν − ∆M is the energy of electron with ∆M = 5.8 MeV being the reaction threshold [41]. Only 9.7% of 8 B neutrinos have energy E ν > 11 MeV but due to strong energy dependence in (50) the corresponding fraction of detected events is 0.9. Therefore, we use the threshold 11 MeV to achieve higher energy reconstruction. For g e we use the width σ E /E = 0.1. With this arrangement σ E /E = 7%, and consequently, the attenuation length equals λ att = 1800 km for the average energy 12 MeV. The nadir angle at which the length of trajectory L > λ att is η att = 82 • . For η < η att the Earth structures on the remote part of a neutrino trajectory become invisible. Results of computations of A N D (η) with the S-R, CRUST1 and SAW642AN density profiles are presented in Fig. 6.
Generic features of the η dependence of A N D are the following: (i) Oscillations in crust: Regular oscillatory pattern for η > η M oho , i.e. η ∼ 85 • − 90 • with decreasing depth due to averaging. The third oscillatory peak can be affected by small density jumps in the crust. This quasi-regular oscillatory pattern is broken at at η M oho .
(ii) Moho interference: At η < η M oho neutrino trajectory crosses the Moho border twice leading to interference of oscillation waves from two crossings. For some models and values of ∆m 2 21 the destructive interference of the waves leads to a dip at η dip (for DUNE) which depends on η M oho . This can also be interpreted as a parametric suppression of oscillations [41].
(iii) Rise of asymmetry: For η < η dip , the asymmetry A N D increases with decrease of η. The increase is due to the fact that for small η the section of the neutrino trajectory in the crust becomes much smaller than the oscillation length, and so the effective initial and final densities (averaged over the oscillation length) become larger, being determined by the mantle density.
(iv). In the region η < η dip there are bump and another dip due to effect of density jumps in the mantle at the depths 400 and 670 km.
(v) The core of the Earth η core = 33 • is not seen practically, producing ∼ 2 effect at η < η core .
We find that about 27000 ν e events (49) can be detected annually with E ν > 11 MeV in the 40 kt fuducial volume according to the CRUST1 model. The crosses show the expected errors A N D (η) after twenty years of data taking. Statistical errors (computed using the exposure function) are taken into account only and no background was considered. As follows from Fig. 6, the largest difference between SAW642AN and S-R models as well as SAW642AN and CRUST1, is in the interval η = 60 • −75 • and it originates mainly from different depths of Moho. The difference equals ∆A N D (η) ∼ 0.006 (15%) which is slightly above 1σ C.L., after 20 years of data taking. The difference between CRUST1 and S-R models is practically negligible. Averaging of A N D (η) over η leads tō A N D = 0.040, 0.040 and 0.043, for CRUST1, S-R, and SAW642AN models, respectively, and precision of measurement ofĀ N D will be 0.002.
The dependence of A N D on η in DUNE experiment computed with SAW642AN model (red line Fig. 6) is similar to that in [41] for the PREM model. It has a dip at η dip = 82 • and then increase of A N D with decrease of η. Another dip appears at η = 44 • . In our present computations (SAW642AN) the dependence A N D (η) is smoother than in [41] below the dip.

THEIA
THEIA is a proposed 30 kT water-based liquid scintillator detector loaded with 1% 7 Li [37]. It will be placed in Homestake. Neutrinos can be detected by the chargedcurrent process ν e + 7 Li → 7 Be + e.
(51) The cross-section of this process is known with high precision [36] [37]. About 5000 events are expected annually with E ν > 5 MeV. For THEIA we assume σ E /E = 12%. Therefore much lower statistics and worser energy resolution make THEIA in the present configuration less suitable for tomography than DUNE.
Since THEIA and DUNE are in the same place the results for A N D (η) are similar (see Fig. 7). The difference between A N D in THEIA and DUNE is due to lower energy threshold in THEIA, which means that effective neutrino energy, and consequently, the oscillation as well as the attenuation lengthes are smaller. This, in turn, leads to different interference effects and lower sensitivity to remore=te structures in THEIA. The difference disappears when the same energy thresholds are taken.
For THEIA maximal difference of A N D (η) computed with S-R and SAW642AN models (and also between CRUST1 and SAW642AN) is about A N D = 0.005, which is smaller than 1σ precision after 20 years of exposure. The difference between S-R and CRUST1 profile results is much smaller. The values of A N D averaged over η with exposure taken into account equal to 0.024 (CRUST1), 0.024 (S-R) and 0.027 (SAW642AN).
The discrimination between the S-R and CRUST1 models can be improved if for each nadir angle η the range of azimuthal angle φ a is divided into two parts: in the first partρ SR >ρ CRU ST 1 , and in the second onē ρ SR <ρ CRU ST 1 . Then calculating A N D in each of these parts separately and summing up moduli of differences one can avoid averaging. The same as in Fig. 6, but for THEIA.

Hyper-Kamiokande
Hyper-Kamiokande (HK) will detect the solar neutrinos by the ν −e elastic scattering with 6.5 MeV threshold [35]. We take σ E /E = 15% as a tentative value. This gives the attenuation length λ att = 700 km for E = 10 MeV.
In Fig. 8, we show the excess of night events computed with FWEA18, SAW642AN and CRUST1 density profiles. For d M oho = 33 km (FWEA18) the nadir angle η M oho = 84.15 • , and the length of the trajectory L = 1300 km, so, remote half of this trajectory will not contribute to the oscillation effect. The dip appears at η dip = 78 • which is intermediate between CRUST1 and SAW64AN.
According to Fig. 8 maximal difference of A N D in HK computed with FWEA18 and SAW642AN: ∆A N D = 0.003, appears in the wide range of nadir angles: η = 10 • − 80 • . The difference is more than 1σ C.L. after 20 year of exposure with fiducial volume 0.56 Mton. For SAW642AN model the η dependence in HK is similar to that in THEIA detector. CRUST1 and FWEA18 have the biggest difference ∆A N D = 0.004 in narror range η = 75 • −80 • . Notice that CRUST1 does not produce the dip which is a model dependent feature. The difference between CRUST1 and FWEA18 as well as SAW642AN is 1σ C.L. The expected averaged asymmetry A N D in HK equals 0.020 (FWEA18), 0.022 (CRUST1) and 0.024 (SAW642AN). Precision of measurements ofĀ N D will be 0.002.
The absolute value of asymmetry is substantially smaller than that for DUNE for two reasons: damping due to contribution from NC scattering, which is 0.76, and difference of averaged energies E HK /E DU N E = 0.75.

MICA
The Megaton scale Ice Cherenkov Array (MICA) is a proposed detector at Amundsen-Scott South Pole station [38] in the same place as ICECUBE. The latitude and longitude of MICA are 89.99 • south and 63.45 • west correspondingly. Crustal structures under Antarctica are not well known due to a lack of seismic data [48], and therefore it is interesting to explore potential of a solar neutrino detector to determine this structure.
The detection is based on the ν − e elastic scattering. In our calculations we took the characteristics of MICA from Ref. [38]: 10 Mton fiducial mass and 10 MeV energy threshold for the kinetic energy of the recoil electron. With these parameters, we find that about 5×10 5 solar νe− scattering events are expected per year. For the energy resolution we use σ E /E = 15%. We consider the MICA detector at a depth of 2.25 km below the icecap (as the Deep Core). The height of icecap at the location of MICA is 2.7 km above the sea level.
The smallest nadir angle for MICA is 66.5 • . About 35% of the neutrinos have the nadir angle in the interval 66.5 • −70 • . These neutrinos propagate through the Earth with a maximal depth of 500 km. For η = 75 • (where the largest difference of A N D from CRUST1 and SAW642AN is expected) neutrinos propagate with a maximal depth of 200 km. Neutrinos reach this angle on May 4 for the first time in a year. According to CRUST1 for η = 75 • , the depth of Moho is 35 km, with the density jump from 2.9 to 3.4 g/cm 3 .
In Fig. 9 we show A N D (η) computed with CRUST1 and SAW642AN models. CRUST1 allows to take into account the Earth density above the sea-level. Since there is no data available for SAW642AN, for this region, we take zero density above the sea-level. After 20 years of data taking MICA will collect 10 7 solar neutrino events, and it will be sensitive to the ice-soil border.  Fig. 9) are real. In this model the surface of the Earth is not spherically symmetric and the density of the Earth above the sealevel is given. Therefore neutrinos enter the Earth at different height from sea-level, which leads to ripples due to change of the baseline with η. Such ripples are far from being detected experimentally. The ripples of A N D are absent in the SAW642AN model (the red curve).

Depandence on ∆m 2 21 ; PREM model results
There is significant difference in values of ∆m 2 21 determined by KAMLAND and from global fit of the solar neutrino data. In this connection we performed computations of A N D (η) using the "solar" value ∆m 2 21 = 5×10 −5 eV 2 (Fig. 10). The changes are twofold: the overal asymmetry increases as 1/∆m 2 21 , i.e. becomes 1.5 times larger than before. The oscillation and attenuation lengths increase by the same factor 1.5. This, in turn, leads to (i) some change of the interference picture, (ii) enhancement of sensitivity to remote structures and bigger densities. As a result, at small η enhancement factor of the asymmetry is bigger than than 1.5.
Let us compare results computed for DUNE with the S-R model for two different ∆m 2 21 (blue line in Fig. 10 and black line in Fig. 6). As expected, for large η the amplitude of oscillations of A N D and its average value are 1.6 times larger than those for large ∆m 2 21 . The dip at 77 • disappears. The peak at 50 • is higher by factor 1.8. For deeper trajectories (smaller η) the enhancement factor is 1.80 − 1.85. The reason for this additional increase of the asymmetry above factor 1.5 is that due to larger oscillation length for deep trajectories the effective intial and final densities (averaged over the oscillation length) become larger. For HK and CRUST1 model the results of ∆m 2 21 change are similar: For shallow trajectories the asymmetry increases by factor 1.5, while for deep trajectories (small η) -by factor 2.  for DUNE and THEIA using the S-R model, and for HK and MICA with CRUST1 model.
Most of the previous computations were performed with PREM model which has two layers in the crust (0 -15) km and (15 -24.4) km and density jumps from 2.6 to 2.9 g/cm 3 at 15 km, and 2.9 to 3.38 g/cm 3 at 24.4 km (Moho). The 3 km layer of water is neglected. In Fig. 12 (upper panel) we compare results of PREM (black line) and CRUST1 (blue line) models for DUNE. The difference is is mainly related to the depths of Moho: η M oho = 48 km for CRUST1, which is two times larger than in PREM. Correspondingly, in the CRUST1 model, the dip of A N D is shifted to smaller η and for η < η dip the asymmetry is smaller. The latter is due to smaller effective density (averaged over the oscillation length) near the detector in CRUST1.
The PREM result is similar to that in [41]. Less profound oscillatory modulations than in [41] are related to different treatment of the energy resolution. As we mentioned before, the PREM model result is close to that of SAW642AN model which has similar depth of Moho.
For comparison in Fig. 12 we show also result for PREM model with outer water layer. That would correspond to a detector near the ocean cost. Large difference appears for η > 88 • i.e. for trajectories in water: the depth of oscillations and averageĀ N D are smaller since they correspond to small water density 1.02 g/cm 3 . Similar situation is for HK Fig. 12 (bottom). According to CRUST1 the dip is absent, A N D is larger in the range η = 75 • − 85 • , while at η < 75 • the asymmetry is 10% smaller (by 0.002) than for PREM.
The results show that usage of PREM model causes up-to 10% relative systematic error in A N D . 30  1. We performed detailed study of the Earth matter effects on solar neutrinos using recent 3D models of the Earth. Interesting and non-trivial oscillation physics is realized which is related to complicated density profiles along neutrino trajectories. The Day-Night asymmetry as a function of the nadir angle has been computed for future experiments DUNE, THEIA and Hy-perKamiokande, as well as for possible next-after-next generation experiment MICA. This allows us to assess feasibility of tomography of the Earth with solar neutrinos.
2. We estimated corrections to A N D of the order ∼ 2 . Corrections ∼ 2 from I 2 can be neglected due to additional small coefficient, while the correction to the oscillation phase can be relevant.
3. We further elaborated on the attenuation effect.
The night excess of events and A N D (η) are expressed in terms of the matter potential and the generalized energy resolution function which, in turn, determines the attenuation factor. This form is the most appropriate for tomography. We have found that inclusion of energy dependence of the boron neutrino flux and crosssection into resolution function improves the resolution, and therefore sensitivity to remote structures. It is the generalized resolution function that determines sensitivity of oscillation results to the density profile.
Further improvement of the sensitivity can be achieved imposing high enough energy threshold for detected electrons. The gain is twofold: (i) The Earth matter effect increases as E; (ii) the attenuation becomes weaker. At the same time loss of statistics is rather moderate.
4. Using recently elaborated 3D models of the Earth we reconstructed the density, and consequently, potential profiles along neutrino trajectories characterized by coordinates of a detector, nadir and azimuthal angles. The key feature of the models is the absence of spherical symmetry. Averaging over φ a leads to dumping of oscillatory modulations.
The key feature of profiles that determines the A N D (η) is the depth of Moho (border between crust and mantle). The depth differs substantially in different models, and furthermore the border substantially deviates from spherical form. 5. Difference of results for different models of the Earth at DUNE and THEIA at Homestake is about 10%. After 20 years of DUNE exposure that would correspond to 1σ. So, the models cannot be discriminated. Similar conclusion is valid for HK.
6. MICA will be sensitive to the ice-soil border. It can discriminate between the CRUST1 and SAW642AN models at 4σ C.L. after 20 years of data taking. 7. With decrease of ∆m 2 21 the overall excess increases as 1/∆m 2 21 . Also η dependence changes which is related to increase of the oscillation length and therefore decrease of the oscillation phase: for deep trajectories the enhancement with decrease of ∆m 2 21 is stronger than 1/∆m 2 21 . 8. The difference of results obtained for Homestake with S-R and CRUST1 from those of PREM model, which was used in most of previous studies, is that the dip in the nadir angle distribution does not appear and for deep trajectories the asymmetry is 10% lower.
In conclusion, future experiment DUNE, THEIA, HK will certainly establish the integrated Earth matter effect with high significance. They may observe some generic features of the η dependence such as dip and slow increase of the excess with decrease of η. However, they will not be able to discriminate between recent models. For this megaton scale experiments like MICA are needed.
The Earth can be considered as a sphere with a very small compared to the Earth radius deviations from the sphere. So, the distance of a given point at the surface from the centre of the Earth equals r E (θ, φ) = 6371 km+ H(θ, φ), where H(θ, φ) is the height from the sea-level of the location. Here, θ and φ are the latitude and longitude of the point respectively. Let us introduce coordinates x, y in the plane perpendicular to the axis of rotation of the Earth and z being along the axis. The axis is tilted by about α = 23.4 • relative to the Earth orbital plane. In these coordinates location of a point on the Earth surface at a given moment of time t is determined by y = r E (θ, φ) cos θ sin(φ + ωt) cos α − r E (θ, φ) sin α sin θ, z = r E (θ, φ) sin θ cos α + r E (θ, φ) cos θ sin(φ + ωt) sin α, where ω is the angular frequency of the Earth rotation. Location (latitude and longitude) of DUNE, THEIA (Homestake) is 44 In all the cases except for MICA we have considered the Earth surface as a perfect sphere (H(θ, φ) = 0), and the detectors located at the surface of the Earth. In the case of MICA, we used CRUST1 model, which allow to take into account H(θ, φ), and the detector is location 2.25 km below the ice surface.
The coordinates of the Earth in the solar system are X = r a cos(Ωt + Φ 0 ), Y = r a sin(Ωt + Φ 0 ), where Ω is 2π/(365.256 days), and r a = a(1 − b cos Ωt) is the distance between Earth and Sun. Here a=1 is the astronomical unit, and b=0.0167 is the eccentricity of the Earth orbit. For the starting point, t = 0, at the 23rd of September the phase equals Φ 0 = − π 2 . Let x D and y D be the coordinates of the detector and x, y and z are the coordinates of the point at which neutrino enters the Earth. The neutrino trajectory inside the Earth is determined by solving the following quadratic equation: where m ≡ Y /X and r 2 D = x 2 D + y 2 D . Taking into account tilt α, the latitude and longitude of the entering point to the Earth and consequently, the trajectory of the neutrino inside the Earth as well as the nadir angle are determined.
To perform a precise calculation of the neutrino trajectory for MICA we use the CRUST1 model. In this case, the Earth is not a perfect sphere. Therefore we solved the quadratic equation first with r E that includes H d the depth of the detector from the sea-level. In this way, we obtained the entrance point of the neutrinos into the Earth, θ 0 and φ 0 . Then we have solved Eq. (54) once again with H(θ 0 , φ 0 ).