{'fraction_non_alphanumeric': 0.045484999627782324, 'fraction_numerical': 0.013790664780763792, 'mean_word_length': 4.525812422871247, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analytical treatment of long-term observations of the day-night asymmetry for solar neutrinos

The Earth's density distribution can be approximately considered piecewise continuous at the scale of two-flavor oscillations of typical solar neutrinos, such as the beryllium-7 and boron-8 neutrinos. This quite general assumption appears to be enough to analytically calculate the day-night asymmetry factor for such neutrinos. Using the explicit time averaging procedure, we show that, within the leading-order approximation, this factor is determined by the electron density within about one oscillation length under the detector, namely, in the Earth's crust (and upper mantle for high-energy neutrinos). We also evaluate the effect of the inner Earth's structure on the observed asymmetry and show that it is suppressed and mainly comes from the neutrinos observed near the winter and summer solstices. As a result, we arrive at the strict interval constraint on the asymmetry, which is valid within quite a wide class of Earth models.


I. INTRODUCTION
The effect of neutrino oscillations in vacuum lies beyond the Standard Model and is thus interesting both from the theoretical and experimental points of view. The oscillations in medium have also been studied since Wolfenstein, who showed that neutrinos acquire a specific flavor-dependent potential due to coherent forward scattering on the matter [1]. As a result, the neutrino propagation in medium should demonstrate the conversion from one flavor into another (i.e. flavor oscillations), even if the vacuum mixing is negligible. This spectacular result is known as the Mikheev-Smirnov-Wolfenstein effect [2] and suffices to explain the deficit of observed solar electron neutrinos [3]. According to Mikheev and Smirnov, the leading-order estimation for the electron neutrino flux depends only on the points where the neutrino was born (the core of the Sun) and absorbed (the detector). However, the properties of the medium between these two points can also slightly affect the flavor composition of the observed neutrino flux, leading, in particular, to the day-night (solar neutrino flavor) asymmetry [4,5]. The latter effect being crucial for the entire flavor oscillations framework, a number of experiments were set up to catch this slight flavor composition variation resulting from the nighttime neutrino propagation through the Earth [6][7][8]. Some experiments, which should be sensitive enough to distinguish the effect, are also planned in the next decade and are now under construction [9].
However, one should hold in mind that, in order to make a conclusion on the expected presence of the daynight asymmetry, one needs not only the experimental data (in terms of event rates and energy spectra) but also a convenient theoretical prediction. First experiments in neutrino oscillations were only aimed at outlining the domains in the neutrino mass and mixing parameter space which do not contradict the observed rates. For such needs, the theoretical estimations obtained using numerical simulations are quite acceptable (see, e.g., [10]). Indeed, although numerical experiments can lack perfect accuracy in some regimes and do not give general results, the calculation of the desired effect can be performed for the entire parameter space without drastic algorithm changes. Therefore, using numerical predictions, one could more or less interpret the experimental data in terms of domains in this parameter space which are excluded. But now that various neutrino experiments have yielded a considerable amount of data on the vacuum neutrino mixing [11][12][13] and its parameters are fixed quite firmly, the question is whether other neutrino oscillation effects to be observed are consistent with the framework accepted so far. Namely, one may ask: how should one compare the present and future experimental data on the day-night asymmetry with the theory of neutrino oscillations in order to make an ultimate conclusion about their agreement within the allowed class of Earth and solar models? The latter issue requires a formalism able to give strict constraints (in a form of an interval with fixed boundaries) on the predicted asymmetry within certain approximations, but valid within quite a wide class of the Earth and solar models. Developing such a formalism constitutes the principal goal of the present paper. In particular, using it, we readily find the constraints on the day-night asymmetry for beryllium and boron neutrinos observed in such experiments as SNO, Super-Kamiokande, and Borexino [6][7][8].
Within the neutrino oscillations framework, it is common to use the Schroedinger-like equation to describe the spatial variations of the neutrino flavor [1,2]. Within the two-flavor approximation, the Schroedinger problem is posed in terms of the 2 × 2 flavor evolution matrix (operator) R(x, x 0 ), whose elements are the neutrino flavor transition amplitudes after traveling from point x 0 to x. The evolution equation and the initial condition read, respectively, where the matrix Hamiltonian is a point-dependent linear combination of the Pauli matrices, a(x) = − cos 2θ 0 + 2EV (x) Here, θ 0 is the vacuum mixing angle, E is the neutrino energy, ∆m 2 is the difference between neutrino masses squared, and V (x) = √ 2G F N e (x) is the Wolfenstein potential, which is proportional to the electron density in the medium N e and the Fermi constant G F = 1.17 × 10 −11 MeV −2 . The constant coefficient λ = ∆m 2 /4E is the reciprocal vacuum neutrino oscillation length, up to the factor of π, ℓ osc = 4πE Equation (1) defines a one-parametric subgroup of SU (2) and hence of SO(3), the translation along x being the group operation. In this sense, equation (1) is analogous to the Bargmann-Michel-Telegdi equation [14] in the spinor representation [15]. It is well known that the solution of the matrix linear ordinary differential equation, such as (1), can be represented as a time-ordered exponential (Dyson expansion) [16]. However, in the case of general N e (x) profile, the solution in terms of an analytical function of the matrix argument (i.e., without a symbolic operation such as time ordering) appears to be too challenging to find. In the recent investigations, considerable progress was made in finding exact solutions of Eq. (1) in certain special cases [17,18]; nevertheless, the general approach to this kind of equation still remains to be approximate.
It was first naturally accepted that numerical simulations could give exhaustive results here and, in a sense, resolve the analytical difficulties that arise in the context of such equations. For instance, a numerical approach was employed in [10,19] to find the domains in the neutrino mixing parameter space where the day-night asymmetry should be potentially observable. However, as mentioned above, estimations obtained with numerical techniques are not always reliable enough. In particular, in the large-∆m 2 regime, which has now proved to be realistic [12,13], the oscillation length (5) becomes small, and numerical evaluation of rapidly oscillating solutions introduces large uncertainties. The numerical time averaging of the flavor observation rates also becomes inaccurate in this regime. This effect is especially strong for low-energy neutrinos, including beryllium neutrinos; probably, due to this fact, the large-∆m 2 and low-energy area was not covered by original numerical simulations [10]. Therefore, obtaining stringent constraints on the day-night asymmetry factor for the continuous measurement of small-length neutrino oscillations favors the analytical approach.
Quite a large number of publications are devoted to the analysis of such approximate analytical solutions. Probably the most effective technique for finding the approximate solutions of matrix linear differential equations, such as (1), is the so-called Magnus expansion [20], which is a generalization of the Baker-Campbell-Hausdorff formula [21]. This approach provides the solution up to any order of approximation, as well as the constraints on the remainder terms [22]. Unfortunately, this technique [23][24][25][26], as well as other general methods (see, e.g., [27,28]), does not readily provide a way to find the solution in its explicit and practically usable form, not firmly fixing the Earth model, namely, N e (x) density distribution. It is thus desirable to find an approximate analytical solution of equation (1), which is valid under quite general assumptions about the electron density profile N e (x). This idea was developed in [29][30][31].
As it was mentioned earlier in this section, in our paper, we are not only aiming at finding the relevant approximate expressions, but also at analyzing their accuracy and applicability domain. Namely, in section II, we find the approximate solutions for the flavor evolution matrix inside the Earth, and then, in section III, we arrive at the observation probabilities for the neutrinos of different flavors. These probabilities are finally subjected to the averaging procedure due to the continuous observation of the solar neutrinos throughout the year (Sec. III B), and the results of this averaging are discussed in section IV. The magnitude of the day-night asymmetry appears to be sensitive to the non-trivial structure of the Earth's crust under the neutrino detector, so the effect of the crust is paid special attention in section III C. The uncertainties of our estimation are also discussed in section IV, in the interesting cases of beryllium and boron neutrinos. We analyze the sources of such uncertainties and finally compare our results with the numerical simulation. The detailed derivation of the approximate solutions for the evolution operator is moved to the Appendix to make the flow of the paper less complicated.

II. THE DENSITY PROFILE AND THE EVOLUTION MATRIX
In our investigation, we use the model electron density profile N e (x) with n − 1 narrow segments, where it changes steeply, separated by n wide though sloping segments. Let us call these segments cliffs and valleys, respectively. Let the cliffs be localized near points x j , j = 1, n − 1, namely, occupy segments [ and have the widths L j ≈ x j −x j−1 . In fact, this kind of model is a good approximation for the Earth's density profile known in geophysics, where it corresponds to the so-called Preliminary Reference Earth Model (PREM) [32,33]. In the present section, we will assume that the above segments can be chosen in such a way that the inequality takes place, i.e. the cliffs are narrow and the valleys are wide compared with the oscillation length. In the following, we will briefly call such density distribution piecewise continuous. Let us note in advance (see Sec. IV for details) that for the energies of beryllium neutrinos (E = 0.862MeV, ℓ osc ≈ 30 km) and, to some extent, boron neutrinos (E ∼ 5 − 10 MeV, ℓ osc ≈ 150 − 300 km), the Earth's layers can be divided into cliffs and valleys satisfying (6), as well as into a number of layers whose widths are of the order of the oscillation length, while the density variations are small. One can show that the contributions of the latter segments can be easily considered in a manner very similar to those of the cliffs, hence we do not consider them until Sec. III C where they will be paid special attention. In the present section, it is enough to note that the valley-cliff classification depends on the neutrino energy. It is also worth saying here that during the night, the neutrino ray traverses different paths through the Earth, thus, the lengths ǫ j and L j vary. However, assumption (6) holds for the most part of the night.
The total flavor evolution operator for our piecewise continuous density profile equals the matrix product of the evolution operators for all cliffs and valleys. Within each of these segments, the two small parameters arise: the first of them, due to the relatively small density of the Earth [29,34], while the second parameter due to the piecewise continuous structure of the density profile, δ = ℓ osc /L j for jth valley, ǫ j /ℓ osc for jth cliff.
The smallness of parameter δ will be discussed in Sec. IV. The evolution matrix for each segment, as well as the total evolution matrix, can be subjected to the following unitary transformations [25]: where ω(x) = a 2 (x) + b 2 and ψ(x) = λ ω(x)dx is the corresponding phase incursion. For the calculations which follow, it is also useful to introduce the effective mixing angle in the medium θ(x) ∈ [0, π/2] [2], which is defined by the expressions In terms of this angle, The transformation with matrices Z ± (x) locally diagonalizes the Hamiltonian H(x) in the point x. It thus makes the complete diagonalization in the homogeneous case N e (x) = const [2]. The transformation with the operator Y (ψ(x)) isolates the effect of the medium inhomogeneity, the transformed evolution matrix R 0 (x, x 0 ) satisfying the equation where the dot over θ denotes the gradientθ Due to the fact that the neutrino detector is homogeneous (θ(x) = 0), Eq. (14) has a well-defined and physically relevant x → +∞ limit for any fixed x 0 . Moreover, asymptotically convergent behavior of such systems of differential equations is stated by the Levinson's theorem [16].
In the homogeneous case, the equation above is trivial, , the slow change of the density N e (x) enables us to use the so-called adiabatic approximation leading to the same result [2] where the remainder term is a (generally speaking, non-diagonal) matrix and ∆η is the variation of the density parameter η over the valley. On the other hand, if the Wolfenstein potential undergoes a considerable change within the narrow cliff [x − j , x + j ] with the phase incursion ∆ψ ≪ 2π, then we get is the jump of the effective mixing angle on the jth cliff. Moreover, one can show that within the more accurate O(ηδ) approximation, the above expressions take the form (see the Appendix) where Now let us write the evolution operator for the whole neutrino path. The neutrinos observed during the day are created in the point x 0 inside the solar core, then travel to the Earth, enter the detector in the point x 1 and are finally absorbed in the point x * inside it. In the nighttime, however, after reaching the Earth in the point x 1 , the neutrinos pass through a number of Earth's layers (valleys) discussed above, and only after that they enter the detector in the point x n and are absorbed in x * . Crossing the Sun-to-vacuum interface, as well as traveling inside the Sun, does not involve steep electron density changes, moreover, the neutrino propagation is highly adiabatic there (see Sec. IV), thus we can treat the whole segment [x 0 , x 1 ] as a single valley and use the adiabatic approximation (16). As it was mentioned earlier, the flavor evolution operator for the whole neutrino path is a matrix product of the evolution operators for each segment (each valley and cliff). By the substitution of the approximate solutions (18) and (19) into representation (9), after some transformations we find the total evolution operator in the form (compare with Eq. (A20)) Here, the subscript 'Sun' refers to the point x 0 inside the solar core, where the neutrino is created, and the evolution operator inside the neutrino detector is denoted R det . The factor n in the remainder term indicates that it contains the sum over all cliffs and valleys. Moreover, we use the following notation: It is also convenient to append definition (24) with If the boundary between the Earth's crust and the detector is abrupt, x + n − x − n ≪ ℓ osc , then µ n vanishes. Quite analogously, µ 1 vanishes for the abrupt vacuum-to-Earth boundary.
By projecting the neutrino state onto the flavor eigenstates, we arrive at the observation probabilities for the electron/muon neutrino where, for nighttime neutrinos, while for daytime neutrinos we have The mixing angle immediately before the detector obviously takes the vacuum value θ 0 in this case. Due to the homogeneity of the detector substance and its smallness compared with the oscillation length, we easily find where the small parameter δ det is the ratio of the detector width L det to the oscillation length ℓ osc . The quadratic remainder terms can obviously be neglected.

A. Finding the probabilities
In order to evaluate the probabilities obtained above, let us first make the averaging over the phase ∆ψ 1 , which corresponds to the neutrino path between the creation point x 0 and the Earth. The region of the neutrino creation is extremely large compared with the oscillation length, thus, after this averaging, cos 2∆ψ 1 = sin 2∆ψ 1 = 0 with a high accuracy. Using this fact together with Eq. (31), after averaging (29) we find e iσ3∆ψ1 e −2iσ2θSun e −iσ3∆ψ1 = cos 2θ Sun , (32) T day = cos 2θ 0 cos 2θ Sun .
The latter expression constitutes the famous result of Mikheev and Smirnov [2]. On the other hand, after averaging over ∆ψ 1 for the nighttime neutrinos, we obtain where we have made use of the fact that matrices e iσ1μ1 and e −iσ2∆θ1 commute up to a negligible term of the order O (∆θ 1μ1 ).
In the above expressions, where L n,j ≡ x − n −x − j is the distance between the boundary of the jth crossed Earth's shell and the detector, measured along the neutrino ray. Finally, by substituting the derivatives (36) and (37) into Eq. (35) and using the fact that sin 2∆ψ n,n = 0, we arrive at the final result which is valid up to the terms of the order O (ηδ 2 ) and quadratic terms O (η 2 ). The term including the product of the detector width ∆ψ det and the oscillating factorμ j cos 2∆ψ n,j is of the order O (ηδδ det ) and is thus omitted.
The above expression provides a generalization of the main result of paper [30] for the case of nonzero-thickness cliffs and detector. It should be stressed, however, that Eq. (39) gives poor information on the effects to be measured. Indeed, the neutrino experiments last for years, and thus, Eq. (39) may only acquire a predictive power after some kind of averaging. The averaging procedure should take into account the axial rotation of the Earth (involving the integration over the nights), as well as its orbital motion around the Sun.

B. Averaging the probabilities
The averaging procedure can be performed analytically, if the oscillation phase incursions ∆ψ n,j vary by much more than 2π during the night. Namely, in this case, one can employ the stationary phase technique (see, e.g., [35]). In the case of the beryllium neutrinos (E = 0.862 MeV) traveling through the Earth, the oscillation length is about 30km, while the depths of the valleys L n,j vary by many hundreds of kilometers, so the oscillation phase variations are indeed large enough to use the stationary phase approximation. To some extent, the same holds for boron neutrinos. However, for both neutrino types, there are layers, to which the stationary phase approximation may not apply. These are the Earth's crust immediately under the detector and (for boron neutrinos) the upper mantle. These layers are discussed in detail in section III C and Appendix A 4 and do not interfere with the picture described in this paragraph.
Let us consider a neutrino traveling through the Earth, which, according to the PREM model [32], consists of a number of concentric spherical shells. The boundary between the valleys x j corresponds to the point where the neutrino crosses one of the interfaces between the Earth's shells; let r j be the radius of this interface (see Fig. 1). Further, the distances L ± n,j between the detector and the points where the neutrino enters/leaves the interface with radius r j are functions of the 'nadir angle' Θ N ∈ [0, π] defined as the angle between the direction to the Sun and the nadir in the point of the detector. In terms of the solar elevation angle Θ s [36], the nadir angle is Θ N = Θ s + π/2. The nadir angle, in turn, is a function of the Earth's axial rotation angle τ ∈ [0, 2π) ('time of day') and the orbital motion angle ς ∈ [0, 2π) ('season'). The dependence of the distances on the nadir angle is easily found to be where the upper/lower signs in (40) correspond to the neutrino entering/leaving the interface r j (see Fig. 1). The inequality (41) ensures that the intersection of the neutrino ray with this interface exists.
In order to find the night average of the electron/muon neutrino observation probabilities, let us note some properties of expressions (39) and (40). First, the number of interfaces crossed by the neutrino is defined via the inequality (41), so the number of the valleys and, thus, the number of terms entering the sums in (39) are changing during the night. Therefore, the night average of Eq. (39) contains the sum over all interfaces j, with the average of the jth oscillating term involving ∆ψ n,j defined as follows: Here, F (Θ N ) is some slowly changing function of the nadir angle and ∆τ night is the total duration of the night in terms of the Earth's axial rotation angle τ , namely, the length of the segment where Θ N (τ ) < π/2 (the Sun is below the horizon). Second, the duration of the night is, in turn, a function of the season ς. On the equinox, e.g., ∆τ night = π, while on the winter solstice, ∆τ night → max. However, the summer nights are just as long as the opposite winter days, so that ∆τ night (ς + π) = 2π − ∆τ night (ς), (43) and the total duration of the nights over all the year is exactly half the year. Therefore, the averaging over the year of N ς days should involve the division by the total duration of the nights, i.e., πN ς . For N ς ≫ 1, the summation over the nights can be replaced by the integration, and the averaging formula for the terms containing the phase incursion ∆ψ n,j takes the form Now we are able to apply the stationary phase technique to such an integral containing the rapidly oscillating exponential. Indeed, let us use the expression, which is valid for smooth functions f (x) and S(x) defined on segment The two leading terms come from the stationary point and the boundary, respectively. However, in the application to the integral (42), the boundary term vanishes. Indeed, the boundary of the integration domain corresponds to the neutrino ray being tangent to the interface with radius r j , hence, ∂ τ ∆ψ n,j ∝ ∂ τ L n,j (Θ N (τ )) → ∞, and the boundary term is absent. On the other hand, the stationary point is obviously achieved at midnight, when the nadir angle Θ N → min (the Sun is in its lowest position), so the integration over the night yields where the two possible signs before iπ/4 correspond to L n,j = L ± n,j (see Eq. (40)). The principal point here is that the second derivative ∂ 2 τ L n,j at midnight is suppressed for the inner Earth's shells, Now let us use the expression of the nadir angle Θ N via the Earth's axial tilt ε = 23.5 • , the latitude of the detector χ ∈ [−π/2, π/2], and the season ς ∈ [0, 2π] [36], cos Θ N (τ, ς) = cos χ sin ς sin τ + cos ε cos χ cos ς cos τ + sin ε sin χ cos ς, where ς = 0 corresponds to the winter solstice in the northern hemisphere. The minimum values of Θ N are achieved at midnights corresponding to τ = τ midnight , Using these expressions, we find the derivative of cos Θ N at midnight, Finally, the integral over the night (47) takes the form On the other hand, the midnight stationary (minimum) values of the nadir angle Θ N (τ midnight ) vary throughout the year (see Eq. (51)), being the smallest on the winter solstice (the darkest midnight) and the largest on the opposite summer solstice (the lightest midnight). Therefore, the right side of Eq. (53) is still containing a rapidly oscillating function of the season ς, and we can perform another isolation of the stationary points, namely, of the two solstices ς = 0, π: ∂ 2 ς cos Θ N (τ midnight (ς)) = sin(ε − χ) tan ε, ς = 0 (winter solstice), sin(ε + χ) tan ε, ς = π (summer solstice), In the above expressions, the Heaviside theta function ϑ r j − r n sin Θ N ensures that the stationary point for the jth interface is really reached, while the additional sign s ′ is explicitly introduced to avoid ±-expressions originating from Eq. (40). The signs specified in Eq. (54) are valid in the northern non-tropical and non-polar latitudes χ ∈ (ε, π/2−ε). Indeed, the typical neutrino detectors (SNO, Borexino, Super-Kamiokande) are situated in temperate latitudes. Here, the midnight solstice nadir angle is Θ N (τ midnight ) = χ + sε and the prefactor N (ς) = cos ε cos χ. With the use of this fact together with the above general averaging formula, one finally arrives at the year average of the night transition probability (39), 2π sin ε cos χ sin(χ + sε) In the above expression, we have omitted the terms resulting from averaging the O (ηδ) terms in Eq. (39), since they are extremely small. The phase incursions ∆ψ n,j should obviously be taken at Θ N = χ + sε, i.e. at solstice midnights.
For the Super-Kamiokande detector, θ = +36.2 • , and the prefactor equals 0.60 and 0.30 for the two stationary points, respectively. Let us briefly note that, in the tropical latitudes |χ| < ε, additional stationary points appear; in particular, on the Equator χ = 0, they correspond to the equinoxes. Two additional points meet on the winter solstice, when χ → +ε (on the Tropic), and one encounters a degenerate stationary point [35]. Although we have found the analytical expressions for the year averages in the tropical and equatorial zones (|χ| ≤ ε), we do not present them here due to their mathematical complexity and to the fact that the actual neutrino detectors are situated in the temperate latitudes. We confine ourselves to saying that, from (61), one can infer the amplification of the winter solstice contribution, as one approaches the Tropic.
One should also mention that, for low-energy neutrinos (E 1 MeV), in certain seasons, the phase incursions ∆ψ n,j for a certain layer j may differ by approximately a multiple of 2π on the successive nights. As a result, the contributions of these nights will not cancel each other, quite similarly to those of the nights near the solstices. This may be called a parametric resonance and, in principle, will lead to local anomalies in the observed neutrino flux, but, as one may see from Sec. IV, the day-night asymmetry for low-energy neutrinos is quite small and its anomalies are even more challenging to observe. Thus, we do not pay these additional effects much attention here.
Finally, we are left with the following conclusion. The terms entering Eq. (39), which contain the oscillating functions of the phase incursions 2∆ψ n,j , are suppressed as O rj rnλLn,j = O rj rn δ after averaging over the year; within the leading approximation, the resulting averages (59) come from the stationary phase points achieved on the winter and the summer solstices. The suppression becomes stronger for the inner Earth's shells.
It is spectacular that all the terms of the order O (ηδ) in the expression (39), including the one corresponding to the detector, become O (ηδ 2 r j /r n ) after the averaging. The terms of the order O (η), which are proportional to ∆θ j , become O (ηδr j /r n ), respectively. Finally, we are left with the average value By substituting this result together with the daytime average value (33) into expression (27) for the neutrino observation probabilities, we arrive at the day-night asymmetry factor A dn ≡ 2( P e,night − P e,day ) P e,night + P e,day = − where T MSW = cos 2θ 0 cos 2θ Sun = a vac a Sun /ω Sun defines the observation probabilities for the solar neutrinos due to the Mikheev-Smirnov-Wolfenstein effect [2] and V (x − n ) is the Wolfenstein potential in the Earth under the detector.
One should hold in mind that the asymmetry factor defined as (63) may not be directly measurable in the neutrino experiments, depending on the detection mechanism. Indeed, definition (63) differs from the one preferred by the experimentalists, the latter being where N day,night is the number of neutrino events observed during the day/night. These numbers may not correspond to the electron neutrino fluxes. For example, scattering experiments, which cannot separate the charged and the neutral current interactions, are unable to directly measure the electron neutrino flux. To compare prediction (63) with such experiments, one should reinterpret the observed event rates N day,night in terms of the electron neutrino fluxes using some theoretical assumptions. Let us provide an example of the connection between the day-night asymmetry factors defined as (63) and (64). Namely, in the case of the neutrino-electron scattering experiment, such as Borexino, the incident neutrinos produce recoil electrons inside the detector, and the scattering cross sections for such processes are well-known [37,38], including one-loop corrections [39]. The resulting ratio of the total electron/muon neutrino detection cross sections is a function of the neutrino energy E (as well as of the threshold T min of the recoil electron detection). For monochromatic beryllium neutrinos and the actual Borexino's threshold, this ratio is [11] σ(ν e ; E)/σ(ν µ ; E) ≈ 4.5, thus, the ratio of the event rates is N day N night = P e,day σ(ν e ) + (1 − P e,day )σ(ν µ ) P e,night σ(ν e ) + (1 − P e,night )σ(ν µ ) .
The 'experimental' day-night asymmetry factor is then easily found to be and, for beryllium neutrinos (T MSW ≈ 0.09), we obtain For boron and other types of neutrinos, which have continuous energy spectrum, the expression for the 'experimental' day-night asymmetry factor involves the integration over the neutrino energy, Here, the cross sections σ(ν e,µ ; E), the 'physical' asymmetry factor A dn , and the Mikheev-Smirnov-Wolfenstein factor T MSW depend on the neutrino energy, and ρ(E) is the normalized energy distribution of incident neutrinos. As one can see, the integration can be easily undertaken numerically using our asymmetry prediction (63) and the expressions for the effective total cross sections, if the recoil electron detection threshold T min is known (see, e.g., [38]).
For those experiments which directly observe the charged-current electron neutrino events, one formally sets σ(ν µ ; E) → 0. Then one has A (exp) dn = A dn for monochromatic neutrinos and for continuous-spectrum neutrinos.

C. The effect of the crust
The estimation (63) shown above is substantially based on the piecewise continuous structure of the density profile and shows that the asymmetry should depend only on the density of rock immediately under the detector, i.e. in the Earth's crust. At the same time, for beryllium neutrinos, the actual width of the crust is comparable with the oscillation length, and neither the valley nor the cliff approximation is valid for this layer. For boron neutrinos, the crust can be considered a cliff, however, the oscillation length becomes comparable with the thickness of the upper mantle, as well as with thicknesses of transition zones (see Fig. 1 and [33]). Moreover, the phase incursions are small within the crust (and the upper mantle), and so are their time variations, hence, we are unable to use the stationary phase approximation. However, we are still able to account for the effect of these near-surface layers on the observed day-night asymmetry factor (63), relying upon relatively small density variation ∆η within them. This feature, together with the bounded layers' thickness, makes it possible to find the closed form of the approximate flavor evolution operator for the crust (the upper mantle) [x + n−1 , x − n ] (see Appendix A 4), This result formally repeats the cliff approximation (19) up to the substitution ∆θ j → β, µ j → α. Using such a substitution, the transition zones, whose widths are comparable with the boron neutrino oscillation lengths, can be safely replaced by the cliffs, preserving the form of expression (39). Finally, to account for the effect of the crust (and the upper mantle), as well as the effects of the transition zones, we should make the following modification in (39): T night = cos 2θ Sun cos 2θ − n + 2 sin 2θ − n n−1 j=1 (∆θ j cos 2∆ψ n,j +μ j sin 2∆ψ n,j ) −4∆ψ det sin 2θ det n−1 j=1 ∆θ j sin 2∆ψ n,j + ∆T night , ∆T night = cos 2θ Sun {2 sin 2θ − n (β cos 2∆ψ n,n−1 + α sin 2∆ψ n,n−1 ) −4∆ψ det sin 2θ det (β sin 2∆ψ n,n−1 − α cos 2∆ψ n,n−1 )} = 2 cos 2θ Sun sin 2θ − Again, this result should be subjected to the averaging procedure to acquire the predictive power. However, due to the boundedness of the cosine, the correction is easily estimated (both before and after averaging), . (76)

IV. DISCUSSION
Let us begin with the review of the approximations we were using in the above calculations of the day-night asymmetry. First, the Sun was considered one big valley, i.e. the adiabatic approximation was used for it. The non-adiabatic corrections are of the order of the adiabaticity parameter γ = |θ(x)|/λω [2]. Using the fact that the typical spatial scale of the solar density variation is associated with the solar core radius R 0 ∼ 0.1R Sun ≈ 7 × 10 4 km [40], we find Thus, these corrections can be safely neglected.
Neglecting the finite width of the detector introduces a relative error of the order O L det ℓosc δ , compared with the leading term (63), as one can see from (39). This correction is minuscule even for beryllium neutrinos (ℓ osc ≈ 30 km), since the detector sizes are now L det ≤ 1 km. Moreover, this type of correction is subjected to additional suppression due to time averaging.
The valley-crust approximation considered in Sec. II relies upon the smallness of parameters η and δ. The relative error of the linear approximation in η is of the order of O (η), hence, this approximation works fine for beryllium neutrinos and quite well for boron neutrinos (see Eq. (7); note that the maximum values specified for η correspond to the inner core, while for typical neutrino detector latitudes, the Sun never descends low enough to shine through it). Thus, for boron neutrinos, the error of the linear approximation is within 10%. Further, the error due to the finite width of the jth valley is estimated as O (∆η j δ j ) = O (∆η j ℓ osc /L j ), where ∆η j is the density variation on the jth valley and L j is the width of the valley. On the other hand, if the width of the kth cliff L k ≪ ℓ osc , the resulting error is O (ηδ k ) = O (η L k /ℓ osc ). Within the PREM model, the cliffs are abrupt (δ k = 0), so we will not discuss the numerical values of the corresponding errors.
If there are deep layers in the PREM density distribution, whose widths are of the order of the oscillation length, they can be considered neither valleys nor cliffs. However, if the density variations ∆η over these layers are small, we can use the approximation considered in Sec. III C and Appendix A 4, which formally reduces these layers to the cliffs (see Eq. (71)). The error of such an approach is measured by the magnitudes of parameters α, β in (71), i.e. the error is of the order of ∆η. The contributions of such layers are additionally suppressed after averaging, due to their depth (see below). In the case of the Earth, there is a layer of this type. Namely, the core-mantle transition zone is about 200 km wide [33], which is of the order of the typical boron neutrino oscillation length. However, the corresponding density variation lies within only ∆η/η ∼ 5 − 7%, and the time-averaged contribution is obviously negligible. For beryllium neutrinos, the core-mantle transition zone can be considered a valley.
The contributions from the Earth's layers, which lie much more than the oscillation length under the detector, are suppressed due to time averaging, the principal nonvanishing contributions being provided by the winter/summer solstice stationary points. These contributions have the relative magnitude O ∆θj ∆θn rj rn ℓosc Ln,j compared with the leading term (63) (see Eqs. (39) and (59)). This suppression is quite strong for dense though deep inner Earth's shells, including its core.
Finally, for the layers which lie within O (ℓ osc ) under the detector and have the widths comparable with the oscillation length, the averaging procedure cannot be performed so as to result in the expression (59). Here, we are only able to follow the approach of Sec. (III C), which results in an unaveraged correction (75). The latter correction, in principle, can be explicitly averaged using numerical integration, which in this case involves only the near-surface Earth's structure. The strict constraint on the uncertainty introduced by neglecting the effect of this structure is given in Eq. (76). This expression predicts the 20% uncertainty for beryllium and ∼ 30% for boron neutrinos. However, it is quite obvious that the time average of the cosine in (75) is sensitive to the Wolfenstein potential in the layers within O (ℓ osc ) under the detector; moreover, (59) shows that this sensitivity asymptotically falls off as the inverse depth of the layer L n,j . Thus, if the resulting expression (63) is used for boron neutrinos, one should substitute for V (x − n ) the Wolfenstein potential averaged over the crust and the uppermost layers of the mantle, using some decreasing weighting function. The resulting potential will be 15 − 20% larger than that immediately under the Earth's surface. For beryllium neutrinos, it suffices to substitute the Wolfenstein potential in the crust.
Thus, we may conclude that, according to (63) and (75), the asymmetry has the order O (η) and is determined by the rock density in the layer of width about ℓ osc under the detector, i.e. in the Earth's crust (as well as in the uppermost part of the mantle for boron neutrinos). The principal uncertainty of expression (63) comes from the fact that the stationary point approximation we used for the analytical time averaging of Eq. (39) is inapplicable to the layers which lie within several neutrino oscillation lengths under the detector. Another correction comes from the time averaging of oscillating contributions of the deep layers; using (59), we can estimate its magnitude as ≤ 3% for beryllium neutrinos and ≤ 10 − 15% for boron neutrinos with E = 10 MeV. All other corrections, taken together, do not exceed 10%.
Using the recent data from SNO, KamLAND, and Borexino collaborations [11][12][13], namely, tan 2 θ 0 ≈ 0.46 and ∆m 2 ≈ 7.6 × 10 −5 eV 2 , and the typical electron densities in the Earth's crust N e(crust) = 1.3 mol/cm 3 [33] and in the solar core N e(Sun) ∼ 100 mol/cm 3 [40], we arrive at the numerical estimation for the day-night asymmetry factor for solar beryllium-7 neutrinos (E = 0.862 MeV) A dn ( 7 Be) = (−4.0 ± 0.9) × 10 −4 , A The uncertainty corresponds to the effect of the Earth's crust, which, as mentioned, can be explicitly evaluated by numerical averaging of Eq. (75). For boron-8 neutrinos with E = 10 MeV, substitution of the crust density into (63) leads to the asymmetry estimation The 'experimental' asymmetry factor for the electron scattering experiment with the recoil kinetic energy threshold T min = 4.5 MeV (such as Super-Kamiokande-III [41]), averaged over boron-8 solar neutrino spectrum [42], is then given by Eq. (69), the fully-numerical calculation yielding A (exp) dn ( 8 B) = 1.9 × 10 −2 . We do not present here the asymmetry predictions for 13 N, 15 O, and pep neutrinos, since their typical energies are around 1 MeV, while the fluxes at least one order smaller than that for beryllium neutrinos [40].
If, for boron neutrinos, one substitutes into (63) the mean density in the near-surface layers, instead of the crust density, the above predictions (79), (80) will be about 10% larger (depending on the weighting function chosen for the mean density evaluation) but never larger than those obtained by substituting the density of the upper mantle. Indeed, Fig. 2b shows that the two analytical curves corresponding to the two 'surface densities' discussed here enclose the fits obtained using two techniques of numerical simulation. One of these simulations involves the numerical time averaging of the leading O (η) terms in expression (39) (the valley-cliff approximation), while the other one includes both the numerical solution of the evolution equation (14) and the subsequent time averaging.
The numerical curves shown in Fig. 2 were calculated for the Gran Sasso laboratory, where the Borexino detector is operating; the curves for Kamioka (Super-Kamiokande) and Sudbury (SNO) lie very close to that for Gran Sasso, so we do not show them in this figure. Instead, a comparison of numerical results for selected detector latitudes is presented in Fig. 3 (again, we do not include the SNO latitude, since the results for it lie very close to those for Gran Sasso). We have also included in Fig. 3 the numerical results for the Northern tropic, since, according to our analytical estimations, near it, the subleading contributions to the asymmetry may become considerable, which come from the solstices (see Eq. (59)). We will address this issue further in this section. From Fig. 3, one may infer that the leading-order analytical estimation (63) is in good agreement with the numerical results even for boron neutrinos. Contrary to these numerical results, however, our analytical estimation is modelindependent, in particular, it does not contain the latitude of the neutrino detector. The possible errors of our predictions can also be easily estimated. In the case of beryllium neutrinos, the estimation of the errors becomes strict enough to result in a fixed-boundary interval (not a confidence interval) for the day-night asymmetry shown in Eq. (78), which is useful for experimental purposes. Our results are also in agreement with numerical day-night asymmetry predictions presented by other authors (see, e.g., [10]).
As predicted, the agreement of the numerical results with the analytical one becomes better for low-energy neutrinos. We also observe that the asymmetry vanishes for low-energy neutrinos, which is a result of the applicability of the adiabatic approximation for such neutrinos. Indeed, within this approximation, the neutrino flavor observation probabilities depend only on the creation and absorption points.
It is worth emphasizing here that such a directly measurable quantity, as the day-night asymmetry factor, does depend both on the neutrino regeneration effect inside the Earth and on the Mikheev-Smirnov-Wolfenstein effect inside the solar core, where the neutrino is created. Another useful physical quantity, namely, the regeneration factor f reg , describes solely the Earth effect in the day-night asymmetry [30], Indeed, one observes that, unlike the regeneration factor (82), the day-night asymmetry (63) depends on the solar effect manifested in the energy-dependent quantity T MSW . As a result, at the energy E ∼ 2.0 MeV, which corresponds to the Mikheev-Smirnov resonance in the solar core, we have cos 2θ Sun = 0 and T MSW = 0, and, as a consequence, the day-night asymmetry factor (63) changes sign (see Fig. 2a); on the other hand, the regeneration factor is always positive. Thus, it is the resonance inside the Sun which makes the asymmetry, observed on the Earth, vanish. From Fig. 2a, 3a, one may observe that beryllium neutrinos (E = 0.862 MeV) would be indeed quite useful for the study of the matter effects in the neutrino oscillations, since they correspond to almost maximum asymmetry magnitude in the unusual domain of its negativity. Moreover, solar beryllium neutrinos are highly monochromatic (in contrast, e.g., to the boron neutrinos) and their flux is considerably larger [40]. However, we are able to conclude that the day-night effect needs at least a 10-20 times improvement of current experimental resolution to be distinguished at a considerable confidence level for such neutrinos. In particular, the result of the Borexino experiment A (exp) dn ( 7 Be) = 1 ± 12(stat.) ± 7(syst.) × 10 −3 was reported in April, 2011 [8] and demonstrates the strong dominance of the uncertainties over the expected effect. One may still hope that the 50 kton LENA detector which should come into operation around 2020 [9] and is expected to observe as much as 10 4 solar beryllium neutrino events per day, could distinguish the day-night effect for such neutrinos. Strictly speaking, the Poisson statistics results in the relative errors ∆N/N ∼ 1/ √ N , and, for a year-long experiment, one reaches the statistical error of the order of 1/ √ 365 · 10 4 ∼ 0.5 × 10 −3 . However, one could employ the adaptive processing of the experimental data, taking into account the expected form of the curve T night (Θ N ) (see (39) and (40)), i.e. the dependence of the asymmetry on the nadir angle. This processing technique may be efficient in extracting the day-night effect from under the noise even for small event rates. Quite a similar technique was recently suggested in [43] as a search tool for periodic time variations in the 7 Be solar neutrino flux observed at LENA.
In Fig. 3b, one can also observe the role of the stationary phase points in the time average of the day-night asymmetry factor. Namely, the numerical curve corresponding to the Northern tropic demonstrates specific oscillations which come from the amplified contribution of the winter solstice stationary point to the year average of the asymmetry factor (see Eq. (59)). The curve for the Kamioka mine, which is situated about 1.5 times closer to the Tropic than Gran Sasso, also demonstrates oscillations, compared to the Gran Sasso curve. In view of this effect, it would be quite prospective to build a detector close the Tropic, which could be able to observe high-energy solar neutrinos with an energy resolution about 0.5 MeV. A favorable place for such a high-technology project could be, e.g., near São Paulo, Brazil (latitude χ = −23 • 33 ′ , i.e. exactly on the Southern tropic!), especially under the potential support of the local university.
One should hold in mind here that the positions of the interference peaks on Fig. 3 substantially depend on the radii of the Earth's shells. Nevertheless, the approximate smoothness of the three numerical curves in Fig. 3 indicates the (approximate) stability of the day-night asymmetry factor with respect to slight variations of the parameters of the PREM model, namely, the radii of the Earth's shells and the density jumps. As seen from Fig. 3, this stability becomes stronger for low-energy neutrinos, as well as for the detectors operating far from the Tropic.
It is also interesting to study the effect of the local Earth's crust inhomogeneities under the detector on the observed day-night asymmetry. Such inhomogeneities could be associated, for instance, with oil-bearing horizons. Let the inhomogeneity be described by the variation δN where the inhomogeneity size δx i ≪ ℓ osc . Then the contribution of this inhomogeneity to the asymmetry factor is given by Eq. (75), where L (i) = x − n −x i is the depth of the inhomogeneity under the detector. One can see that, although the regeneration effect is stronger for higher-energy neutrinos, this effect is insensitive to the near-surface local inhomogeneities for such neutrinos due to the smallness of the sine and ratio 2πδx i /ℓ osc in the above expression. Therefore, exploration of the Earth's crust based on the neutrino oscillations would require an extreme improvement of current measurement techniques [44], and its prospects seem obscure in the nearest future.

V. CONCLUSION
Let us make a brief summary of our initial goals concerning the analytical approach to the day-night asymmetry and the results of our investigation. We have attempted to develop a framework able to give interval constraints on the day-night flavor asymmetry, which are independent of the details of the density distribution inside the Earth. Of course, some approximations should be made to analytically obtain such general results; in our case, the principal assumptions were the relatively small density of the Earth (η ≪ 1) and its spherically-symmetric layered structure (manifested in the small parameter δ). Although the actual approximation parameters may be not quite small, using the framework developed, we can readily estimate the corresponding inaccuracies in our predictions (as done, e.g., in Sec. IV).
Our analysis shows that the day-night asymmetry is insensitive to the structure of the deep Earth's layers, including its core; we found that this sensitivity falls off as the inverse layer's depth 1/L n,j (see (59)). The day-night effect averaged over time (i.e. the effect directly measured in neutrino experiments) is principally determined by the mean electron density of the Earth within 1 − 2 oscillation lengths under the neutrino detector, i.e. in the Earth's crust for beryllium neutrinos, as well as in the upper mantle for boron neutrinos. The corresponding leading-order analytical curves plotted for these two electron densities are shown in Fig. 2 and 3 (the solid and the dashed line, respectively), together with the results of the numerical simulations. One may see that the leading approximation works quite well within the energy range E ∼ 0.5 − 12 MeV under major interest in the field. Moreover, in Fig. 2b, one may notice that high-energy neutrinos, as it was expected, 'sense' the deeper Earth's layers. Indeed, it is indicated by the fact that the high-energy segment of the numerical curve approaches the dashed theoretical curve plotted for the (higher) upper mantle density. If, in view of this fact, one substitutes into the leading-order expression (63) the mean Earth's density within 1 − 2 oscillation lengths under the detector, the resulting estimation will agree with the numerical one within 10%. Such an accuracy of the day-night effect measurements is yet to be achieved in the future experiments, such as LENA [9].
Further, the theoretical analysis predicts next-to-leading-order corrections to the day-night effect, which may arise due to the stationary phase points during the year, at which the nighttime neutrino oscillation phase freezes. These points occur on the winter and the summer solstices and are spectacular for the fact that their contributions do not vanish after (arbitrarily) long observations. These contributions are very small for low-energy neutrinos, as well as in the temperate latitudes, however, in the tropical latitudes, they are quite distinguishable. One can see these stationary point contributions in Fig. 3b (high-energy segment of Fig. 3), especially for the Tropical curve (latitude χ = +23.5 • ) demonstrating oscillations relative to the other curves plotted for the more temperate latitudes. With the present energy resolutions reaching 0.5 MeV [41], only the event rates are yet too small to observe this effect.
It is also worth mentioning that the good agreement of our analytical predictions with the numerical simulations in a surprisingly wide range of neutrino energies (i.e. oscillation lengths) is also a byproduct of a number of specific features of the actual Earth's density profile. For instance, the large density jumps are lying very deep inside the Earth; the crust is quite thin (and, as mentioned, is not a cliff for beryllium neutrinos), however, its density is quite low; the layers' widths are incomparable, etc. The framework developed in the present paper, however, is able to reveal the situations in which the 'universality' of the prediction (63) will not hold; one needs only the general features of the density distribution to make such a conclusion.

Valleys
In the valleys, we have a smooth and bounded functionθ(y) and a rapidly oscillating matrix exponential e 2iσ3ψ(y) = cos 2ψ(y) + iσ 3 sin 2ψ(y). Then, applying the double integration by parts, we arrive at x x0θ (y)e 2iσ3ψ(y) dy = (1 − D y )θ(y) 2iσ 3 λω(y) e 2iσ3ψ(y) x x0 + x x0 D 2 yθ (y) · e 2iσ3ψ(y) dy, (A2) Let ∆η be the total variation of the density parameter η on the valley. Then the variation of the effective mixing angle θ(x) has the order O (∆η), and all the gradients in the above expressions are suppressed by powers of the small parameter δ = ℓ osc /L = π/λL, where L is the width of the valley and ℓ osc is the oscillation length. Namely, Using these estimations, one can readily show that then the matrix norm of the remainder integral in Eq. (A2) x x0 D 2 yθ (y) · e 2iσ3ψ(y) dy ≤ This expression can also be derived using mathematically strict stationary phase technique [35].

Cliffs
On the cliffs, the change of the effective mixing angle is of the order O (η), while the spatial scale of this change is quite small. Therefore, the change of the phase of oscillations ∆ψ ≪ 2π, and we can expand the exponential in the Now, using the boundedness of the total variation of the mixing angle on the cliff, we obtain x x0θ (y)e 2iσ3ψ(y) dy = e 2iσ3ψ(x0) x x0 (1 + 2iσ 3 λ(y − x 0 ))θ(y)dy + O (ηδ 2 ).
This finally leads to the cliff approximation for the evolution operator

The Earth's crust
As one could see from the main flow of the paper, for neutrinos with energies E ∼ 1 − 10 MeV, the Earth's crust (and, for E ∼ 10 MeV, the upper mantle) cannot be considered either a valley or a cliff, since its width is comparable with the oscillation length, and the parameter δ is of the order of unity. Here, however, another approximation is useful, which takes into account the small density variation ∆θ over these layers, as well as their bounded thickness.
Note that the form of this approximation (A24) coincides with the cliff approximation (A17), up to the coefficient substitution ∆θ → β, µ → α. In particular, the cliff approximation is restored in the δ → 0 limit.