Physics potential of the ESS$\nu$SB

The ESS$\nu$SB project proposes to base a neutrino ''Super Beam'' of unprecedented luminosity at the European Spallation Source. The original proposal identified the second peak of the oscillation probability as the optimal to maximize the discovery potential to leptonic CP violation. However this choice reduces the statistics at the detector and penalizes other complementary searches such as the determination of the atmospheric oscillation parameters, particularly the octant of $\theta_{23}$ as well as the neutrino mass ordering. We explore how these shortcomings can be alleviated by the combination of the beam data with the atmospheric neutrino sample that would also be collected at the detector. We find that the combination not only improves very significantly these drawbacks, but also enhances both the CP violation discovery potential and the precision in the measurement of the CP violating phase, for which the facility was originally optimized, by lifting parametric degeneracies. We then reassess the optimization of the ESS$\nu$SB setup when the atmospheric neutrino sample is considered, with an emphasis in performing a measurement of the CP violating phase as precise as possible. We find that for the presently preferred value of $\delta \sim -\pi/2$, shorter baselines and longer running time in neutrino mode would be optimal. In these conditions, a measurement better than $14^\circ$ would be achievable for any value of the $\theta_{23}$ octant and the mass ordering. Conversely, if present and next generation facilities were not able to discover CP violation, longer baselines and more even splitting between neutrino and neutrino modes would be preferable. These choices would allow a $5 \sigma$ discovery of CP violation for around a $60\%$ of the possible values of $\delta$ and to determine its value with a precision around $6^\circ$ if it is close to $0$ or $\pi$.


Introduction
After the discovery of a non-zero θ 13 [1][2][3][4][5] the emerging picture from the last decades of neutrino oscillation searches consolidates a structure for the PMNS matrix [6][7][8][9][10] describing lepton flavor mixing strikingly different from its CKM counterpart in the quark sector, making the Standard Model favour puzzle even more intriguing. Far from the hierarchical structure described through the tiny mixing angles of the CKM, large mixing angles characterize the lepton mixing. The "atmospheric" mixing angle θ 23 is presently compatible with maximal mixing as well as with a large but non-maximal value in either the first or the second octant. Similarly, the "solar" mixing angle θ 12 is around 33 • and only θ 13 ∼ 8 − 9 • is relatively small and its value is still comparable in magnitude to the Cabibbo angle, the largest in the CKM. The large mixing opens the window to the present and next generation of neutrino oscillation experiments to tackle new questions that could provide answers to fundamental open problems.
Present experiments such as T2K [11,12] and NOνA [13] have started to provide the first hints on the potentially CP-violating phase δ. The discovery of the violation of the particle-antiparticle symmetry in the lepton sector would be extremely suggestive, given that CP-violation is a necessary ingredient to explain the matter over antimatter excess to which we owe our existence and that the CKM contribution has been shown to be insufficient [14,15] for this purpose. Similarly, present neutrino oscillation experiments already show some preference for normal ordering (positive ∆m 2 31 ) with respect to inverted ordering. This parameter is a fundamental input to combine with the searches for the neutrinoless double beta decay process in order to probe the Majorana nature of neutrinos. Finally, present experiments as well as their successors T2HK [16] and DUNE [17] will also provide even more precise measurements of the oscillation parameters that could hold the key to discriminate among different flavour models addressing the flavour puzzle.
The European Spallation Source (ESS) at Lund provides an opportunity to build a newgeneration, long-baseline neutrino oscillation experiment with an unprecedented neutrino luminosity through an upgrade of the ESS Linac [18]. Its 2.5 GeV protons would lead to a rather low energy neutrino flux, between 200 and 600 MeV. This energy range is very well suited for a water Cerenkov detector of the MEMPHYS type [19,20]. In Ref. [18] a greenfield study optimizing the physics reach to leptonic CP-violation was performed for this ESS neutrino Super-Beam facility (ESSνSB). Interestingly, the outcome of this optimization study was that the best baseline at which to study the neutrino beam from the ESS facility at a MEMPHYS-type detector would be between 400 and 600 km. Two candidate mines that could host the detector were identified: Garpenberg at 540 km and Zinkgruvan at 360 km from the ESS site. This choice makes the ESSνSB design unique, as the neutrino flux observed by the detector mainly corresponds to the second maximum of the ν µ → ν e oscillation probability, with a marginal contribution of events at the first oscillation peak.
For the value of θ 13 = 8.6 • currently preferred [21] by Daya Bay [22] and RENO [23], the "atmospheric" term of the ν µ → ν e oscillation probability [24], which is governed by oscillations driven by the large frequency ∆m 2 31 and with an amplitude sin 2 2θ 13 , dominates over the sub-leading "solar" term driven by ∆m 2 21 with amplitude sin 2 2θ 12 at the first oscillation maximum. Thus, the interference between the two, which is the only term dependent on the yet unknown CP-violating phase δ, will also be a sub-leading contribution to the full oscillation probability at the first peak and potentially hidden by systematic uncertainties. Conversely, at the second oscillation maximum the slower "solar" oscillation has had more time to develop and thus the CP-violating interference term can give a significant contribution to the oscillation probability, thus increasing the sensitivity to CP violation [25].
The price to pay in order to observe the oscillation probability at its second maximum is high. Despite this being the optimal choice to maximize the dependence of the oscillation probability on the leptonic CP violating phase, the ratio of the oscillation baseline to the neutrino energy (L/E) needs to be a factor 3 larger compared to the first maximum. This implies roughly an order of magnitude less statistics than if the experiment had been designed at the first peak. Indeed, the neutrino flux decreases with L −2 from the beam divergence and the neutrino cross section and beam collimation increase with the neutrino energy. Despite the unprecedented neutrino luminosity from the upgraded ESS linac and the megaton-class MEMPHYS detector, only around 100 signal events for each beam polarity would be accumulated after 10 years data taking (2 years in neutrinos and 8 years in antineutrinos) at the 540 km Garpenberg baseline (see Fig. 7 of Ref. [18]). Conversely, the 360 km Zinkgruvan baseline has a 2.25 times larger neutrino flux. However, the neutrino spectrum for this baseline is rather centered at the first oscillation minimum while the first and second peaks are sampled by the high and low energy tails respectively. Overall this gives similar statistics at the second oscillation maximum when compared to the Garpenberg option, but also some additional statistics at the first peak and in between.
For the ESSνSB the increased dependence on the CP violating phase of the probability is well worth the loss of precious neutrino events at the second maximum. Indeed, it could provide unprecedented discovery potential to leptonic CP-violation or the most precise measurement of the corresponding phase after discovery, which could be instrumental in tackling the flavour puzzle. Moreover, as pointed out in Ref. [25] and as we will elaborate in later sections, this choice also makes the physics reach much more resilient against unexpected sources of systematic errors, since the signal, while small, has a leading dependence on the unknown parameters. Conversely, statistics will be the bottleneck of the ESSνSB physics reach and thus longer periods of data taking would greatly increase its capabilities.
On the other hand, other potential oscillation searches, different from the CP violation search, will be negatively impacted by the choice of the second oscillation maximum baseline. In particular the sensitivity to the octant of θ 23 is severely reduced by this choice. Indeed, this measurement mainly relies on the "atmospheric" term of the oscillation probability, which is leading at the first maximum instead, together with θ 13 information from reactor measurements and ∆m 2 31 and sin 2 2θ 23 from ν µ disappearance. Similarly the ν µ disappearance data and hence the precise determination of ∆m 2 31 and sin 2 2θ 23 are negatively affected by the choice of the second oscillation maximum. The lack of knowledge on the octant of θ 23 can lead to "octant degeneracies" [26] that in turn somewhat limit the CP discovery potential of the ESSνSB [27]. The sensitivity to the mass ordering is also limited at the ESSνSB given the small matter effects from the low energy and short baseline. However, since these matter effects are small, the resulting "sign degeneracies" [28] do not compromise the sensitivity to δ of the facility [18,27].
A very effective and convenient way of increasing both the octant and mass ordering sensitivity of a neutrino Super Beam experiment is to combine the signal from the neutrino beam with the huge atmospheric neutrino sample that can be collected at such a detector [29,30]. In the case of the ESSνSB this combination is particularly synergistic. Indeed, the atmospheric neutrino sample can provide not only significantly increased sensitivity to the octant and the mass ordering to solve parametric degeneracies, but also improved precision to ∆m 2 31 and sin 2 2θ 23 which is otherwise one of the main drawbacks of the setup.
In this work we will combine the observation of the ESSνSB flux tuned for the second maximum of the ν e appearance probability with the complementary atmospheric neutrino data, more strongly dominated by the first maximum and ν µ disappearance, and characterized by stronger matter effects. We will explore how the physics reach of the facility improves when beam data is considered together with the atmospheric neutrino sample and then review the optimization of the ESSνSB facility using both data sets. Finally, we will discuss which sources of systematic errors among the ones considered impact the final sensitivity more significantly.
This paper is organized as follows. In Section 2 we discuss the peculiarities of the neutrino oscillation probability and the appearance of parametric degeneracies when observing at the second oscillation maximum. In Section 3 we describe the experimental setup considered and the details of the numerical simulations performed. Section 4 describes the results of the simulations and in Section 5 we present our conclusions and summarize our work.

Measurements at the second oscillation peak
The determination of the oscillation parameters at beam experiments is, in general, hindered by the appearance of degenerate solutions, cf. e.g., Refs. [31][32][33][34][35]. These degeneracies have been extensively studied for the experimental setups of T2HK [36][37][38][39][40][41] and DUNE [17,36,40,[42][43][44][45][46][47][48][49][50][51] (and also their combination [52,53]). As stated in Section 1, the L/E range which the ESSνSB focuses on is different from those of other forthcoming experiments, 1 Therefore, here we will discuss the peculiarities of ESSνSB and the differences from other experiments in the determination of the oscillation parameters before presenting our numerical results. The ν e appearance oscillation probability in matter is given by [24] (see also [65][66][67]): where ∆ ij ≡ ∆m 2 ij /2E,J = c 13 sin 2θ 12 sin 2θ 23 sin 2θ 13 , A = √ 2G F n e is the matter potential with n e the electron density and G F the Fermi constant, andB ∓ ≡ |A ∓ ∆ 13 |. In this expression the only dependence in the CP violating phase δ appears in the last term, which is the interference between the "atmospheric" oscillation in the first term and the "solar" in the second. Since sin 2θ 13 ∼ 0.3 while ∆ 12 L ∼ 0.05 at the first oscillation peak, the "atmospheric" term tends to dominate the oscillation probability and the interesting CP interference is only subleading. Conversely, at the second oscillation maximum ∆ 12 L ∼ 0.1 so that the dependence on δ of the oscillation probability is much higher which allows to improve the sensitivity to this parameter [25]. This can be seen in Fig. 1 where the change in the probability upon changing the values of δ is much more significant at the second peak maximum compared to the first.
In Eq. (2.1) the leading dependence on the mass ordering comes from the "atmospheric" term, as it goes as the inverse of the square ofB ∓ . For E ∼ |∆m  resonance which will produce an enhancement in neutrinos against antineutrinos or viceversa depending on the mass ordering. For a typical average matter density of 3.0 g/cm 3 one finds that the approximate energy for this resonance to happen is E ∼ O(GeV). Given that the peak of the flux for ESSνSB happens at E ∼ O(100) MeV (see Fig. 1), the importance of the matter effects and hence of the sensitivity to the mass ordering for this facility is not expected to be significant. The bi-probability plots [28] shown in Fig. 2 help to illustrate the degeneracy problem at the ESSνSB experiment. Here all oscillation parameters other than δ, the octant of θ 23 , and the sign of ∆m 2 31 are fixed at the current best fit values [21], and the matter density along the neutrino baseline is assumed to be constant with an average density of 3.0 g/cm 3 . The baseline length L and the neutrino energies E are set to L = 540 km (ESS-Garpenberg) and E = {280, 380, 480} MeV. The ellipses show the variation of the appearance probabilities for the neutrino and antineutrino channels from changes in δ. The four ellipses in each plot correspond to the different choices of the octant of θ 23 and the mass ordering. When the ellipses overlap sharing the same region in the P (ν µ → ν e )-P (ν µ →ν e ) plane, the same oscillation probabilities can be obtained by changing δ, the octant of θ 23 and/or the mass ordering, implying the existence of degenerate solutions.
Let us first focus on the middle plot with E = 380 MeV where the oscillation probabilities are close to the second maximum, |∆m 2 31 |L/(4E) ∼ 3π/2. The centres of the ellipses are located on the CP conserving line P (ν µ → ν e ) = P (ν µ →ν e ), which reflects the fact that the matter effect, which could induce an explicit difference between the neutrino and antineutrino oscillation probabilities unrelated to the intrinsic CP violation from δ, is irrelevant for this energy and baseline. The major axes of the ellipses extend widely along the diagonal line orthogonal to the CP conserving line. This means that the CP violating term proportional to sin δ in Eq.(2.1) is very relevant in the oscillation probability for this energy and baseline, leading to the improved CP sensitivity at the second oscillation peak.
The "fake" CP violation effect due to the matter effect separates the two ellipses with opposite mass ordering at the first oscillation maximum, where T2HK focuses on, causing the δ-sign(∆m 2 31 ) degeneracy in the CP violation search, cf. the right most plot in Fig. 3. Conversely, the CP violation search at the second oscillation maximum is not noticeably affected by the matter effect [27,68]. Changing the value of θ 23 , the ellipses almost keep the same shape and move in parallel along the CP conserving line, which causes the δ-θ 23 degeneracy [33,34].
The vertices of the ellipses are located at δ = {π/2, −π/2}, where the oscillation probabilities do not change much with a change of δ. As a consequence, the precision in the determination of δ becomes worse close to the oscillation maxima [69]. In other words, since the two points with δ and π − δ on an ellipse are close to each other around δ = {π/2, −π/2}, it is hard to separate them [69]. Although at the probability level from Fig. 2 the expectation would be that this quasi-degeneracy effect occurs similarly at δ = π/2 and δ = −π/2, the numerical simulations we will report in Section 4 show that the ESSνSB suffers this effect more severely at δ = −π/2 than at δ = π/2. This is due to the significant difference in event rates between these two points. Indeed, for δ = −π/2, the oscillation probability for neutrinos is enhanced while the antineutrino one is suppressed. Since both the flux and the cross section are also smaller for antineutrinos, this strongly penalizes the measurement at δ = −π/2 since the antineutrino sample is essentially lost given that the event rate at the second oscillation peak is already necessarily small. On the other hand, at δ = π/2, the oscillation probability for neutrinos is suppressed, but the larger cross section and flux compensate for it and prevents such a big loss of sensitivity.
In the energy region that the ESSνSB focuses on, the oscillation phase changes rapidly. As a consequence, the shape and location of the ellipses changes very significantly even within the same energy bin. In Fig. 2, we also show the bi-probability plots with E =280 and 480 MeV where the oscillation probabilities are approaching the minima, which are also well-covered by the ESSνSB flux. The ellipses are not distributed symmetrically to the CP conserving line, which means that, contrary to the second peak, matter effects do have some impact on the oscillation probabilities. However, this impact is still subleading, given the rather low energy, and does not shift the energies where the extrema are located, cf. Fig. 1. As a result, the two ellipses for the different mass hierarchies are not separated in the entire energy region.
The drastic shape change of the ellipses when varying the energy is largely due to the ratio of the sin δ and the cos δ terms in the oscillation probability, see Eq. (2.1). The sin δ term is most significant close to the oscillation peak with |∆m 2 31 |L/(4E) 3π/2 for E 380 MeV. As the probabilities depart from the maximum, the major axes of the ellipses start following along the direction of the CP conserving line, which means that the cos δ term increases in importance as we approach the minima with |∆m 2 31 |L/(4E) π (right panel of Fig. 2 In the left and the right plots, the ellipses with different mass orderings intersect each other at points with different values of δ at different energies. Therefore, in principle, with precise enough measurements at various energies, one could determine the value of δ and the sign of ∆m 2 31 separately. However, the oscillations are too fast for the ∼ 100 MeV resolution achievable at these energies with a water Cerenkov detector to resolve and also the event rate at the second maximum is not large enough to perform a very fine binning. Thus, it is not possible to track the rapid oscillations in Fig. 1, although some mild sensitivity to the mass ordering can be achievable. A large overlap between the two ellipses with different mass orderings and different octants at the oscillation maximum (middle panel in Fig. 2), where most of the statistics is concentrated, suggests that the mass ordering sensitivity at the beam experiment is affected by the octant degeneracy.
The ellipses for different octants barely separate in the entire energy region, which implies a rather poor sensitivity to θ 23 in the appearance channel leading to octant degeneracies that can spoil both the determination of δ and of the mass ordering at the ESSνSB. Conversely, for experiments focusing on the first maxium the two ellipses for different octants are more separated [27], cf. the right panel in Fig. 3. Therefore, we will explore the impact of the addition of the atmospheric neutrino data collected at the far detector of the ESSνSB to the beam data since atmospheric neutrinos can provide both sensitivity to the θ 23 octant and the mass ordering helping to lift parametric degeneracies [29,30].
The mass ordering sensitivity from an observation of atmospheric neutrinos comes from the oscillation signals driven by ∆m 2 31 and the matter effect (first term in Eq. (2.1)) and therefore, it does not depend on the value of δ. On the other hand, the sensitivity is better for θ 23 in the second octant than the first octant, since the term is proportional to sin 2 θ 23 [70].
If the shorter baseline L = 360 km (ESS-Zinkgruvan) is instead considered, the neutrino flux at the high energy tail up to E ∼ 600 MeV covers the first oscillation maximum. This situation corresponds to the bi-probability ellipses presented in the right panel of Fig. 3, which show the same shape and position characteristic of other experiments located at the first oscillation maximum such as T2HK. The matter effect is not significant enough to completely separate the two mass orderings. In the relevant energy range (200-600 MeV), the oscillation probabilities go from the first maximum (right panel) to the first minimum (middle panels) and to the second maximum (left panel). The leftmost panel with E = 250 MeV, where the second oscillation peak would be located, looks very similar to that with E = 380 MeV in the case of L = 540 km. The ellipses for the different mass orderings are separated more clearly in the case of L = 360 km than L = 540 km in a large energy region, which leads to a slightly better sensitivity to the mass ordering even though the baseline is shorter. From the information at the first oscillation maximum, the ESSνSB with L = 360 km also has better sensitivity to θ 23 than the L = 540 km option, so that it is expected that the longer baseline option will benefit more from the addition of the atmospheric neutrino data, which helps to determine θ 23 and its octant.

Simulation and experimental details
The simulation of the ESSνSB data has been performed with the GLoBES software [71,72]. We have assumed that the neutrino beam will shine on a near and a far detector to reduce the systematic uncertainties [18]. The far detector is a 1 Mt MEMPHYS-like water Cerenkov detector [20], while the near detector has been assumed to be identical to the far detector in terms of efficiencies and background rejection capabilities with a fiducial mass of 0.1 kt. The response of the detectors has been implemented through migration matrices, both for the signal efficiency and the background rejection from Ref. [20]. A beam power of 5 MW with 2.5 GeV protons and an exposure of 1.7 × 10 7 operating seconds per year has been assumed [18]. The fluxes have been simulated explicitly at 1 km for the near detector [73], accounting for possible geometrical effects since the source cannot be considered point-like, as well as for 100 km (and consequently rescaled) for the longer baselines considered for the far detector [18]. The event rate peaks around O(100) MeV energies (see Fig.1), so the dominant contribution to the cross section will be in the quasi-elastic regime (QE). For the cross section we use the results from the Genie [74] tune G18 10a 00 000.
We have assumed a total running time of 10 years. Nonetheless, we will also study the  Table 1: Systematic uncertainties for a super beam as described in Ref. [36] for two different scenarios, the "Optimistic" one and the "Conservative" scenario where systematics are larger.
dependence of the physics reach on the relative running time spent in positive and negative focusing in order to optimize it for the measurement of CP violation. Likewise, although the preferred location of the far detector for the ESSνSB is the Garpenberg mine at 540 km [18], different baselines, with emphasis in the alternative Zinkgruvan option at 360 km, will be studied to address the optimal choice. Finally, we will also study how the CP discovery potential depends on the total exposure. Throughout all the simulations we adopt the same treatment of the systematic errors from Table 1 as in Ref. [36]. Unless otherwise specified, we will assume the "Optimistic" systematics from the first "Opt." column in Table 1 although we will also show how the results are affected when the more conservative ones in the second column "Cons." are considered instead. Correlated systematics between near and far detector will be assumed, uncorrelated between neutrinos and antineutrinos and different flavours. Additionally, to account for the uncertainty in the cross section between the near and far detector, arising from the different flavour composition of the beam (mainly ν µ in the near site and ν e for the signal in the far detector), a completely uncorrelated systematic is included for their ratio (last row of Table 1).
We have added to the resulting χ 2 a gaussian prior with the central values and 1σ errors from Ref. [21] for "solar" and "reactor" parameters. For the "atmospheric" parameters we set a prior on sin 2 2θ 23 and |∆m 2 31 | given that the octant for θ 23 and the mass ordering are still unknown. Since the determination of these two parameters comes primarily from atmospherics, when adding this sample to the beam data no prior has been added on θ 23 and ∆m 2 31 . The simulation of the atmospheric neutrino sample in MEMPHYS is the one used in the analysis from Ref. [30] where the neutrino fluxes at Gran Sasso from Honda calculations [75] were used. This is a conservative estimate as fluxes become larger at higher geomagnetic latitudes such as Garpenberg or Zinkgruvan.

Results
In Fig. 4 we show the impact on the CP discovery potential of the ESSνSB before (dashed lines) and after (solid lines) the inclusion of the atmospheric sample for the Zinkgruvan (360 km) and Garpenberg (540 km) options in the left and right panels, respectively. The plots represent the ∆χ 2 with which CP conserving values of δ = 0 or π can be disfavoured as a function of the true value of δ. The ∆χ 2 can be interpreted as the significance for exclusion of CP-conserving values (and hence evidence for CP violation) as long as the assumptions behind Wilks' theorem hold [76]. Deviations from these assumptions can be sizable for presently running experiments, but are expected to be smaller for next generation facilities [77]. Even though the sensitivity of the atmospheric neutrino dataset to δ is almost negligible, the improvement of the ESSνSB physics reach upon its inclusion is quite remarkable. The improvement is generally larger for the longer 540 km baseline than for the Zinkgruvan 360 km option. This is in line with the expectations discussed in Section 2 of the atmospheric sample being more complementary to the beam information at the longer baseline. Indeed, at the second oscillation maximum the ν µ disappearance oscillation is not sampled as efficiently as at the first peak and this deteriorates the determination of the atmospheric oscillation parameters θ 23 and ∆m 2 31 , which play an important role in the measurement of δ. Conversely, the 360 km baseline has higher statistics and some events also cover the first oscillation maximum such that the atmospheric oscillation information is less complementary and the gain upon its inclusion is less noticeable. From these results we can conclude that the ESSνSB setup combined with the atmospheric neutrino sample would be able to rule out CP-conserving values of δ for ∼ 60% (∼ 55%) of the possible values of δ at the 5σ level regardless of the octant and the mass ordering when observing at the 540 km (360 km) baseline. Figure 4 also shows that the gain in CP discovery potential is much more pronounced in some particular regions of the parameter space, especially for δ < 0 and θ 23 in the first octant or δ > 0 and the second octant. In these examples the dotted curves for beam only often show a kink that reduces the slope and the values of δ for which CP-violation could be discovered with high significance. Conversely, the corresponding solid curves with atmospheric data either do not display the kink or develop it at higher significance so that the resulting CPdiscovery potential is much larger. These kinks occur due to the presence of an unresolved octant degeneracy at a CP-conserving value of δ that prevents drawing conclusions regarding CP violation. When atmospheric data is added, the sensitivity to the octant improves and these degeneracies are either lifted or only show up at much higher significance.
This situation is illustrated in Fig. 5, where the allowed regions at the ∆χ 2 = 25 level are shown in the δ-sin 2 θ 23 plane. The left (right) panels assume the true values δ = −40 • (δ = 150 • ), sin 2 2θ 23 = 0.418 (sin 2 2θ 23 = 0.582) and normal ordering. As can be seen, when only the beam information is taken into account (blue curves), an octant degeneracy that spreads the allowed region towards CP conserving values appears. Conversely, the atmospheric data on their own (red curves) have no capability to determine δ at all, but can instead rule out the wrong octant of θ 23 . Thus, the combination of the two data sets (black curves) very significantly improves the CP discovery potential of the facility in these areas of parameter space. The dotted lines correspond to "sign" degeneracies with the opposite mass ordering to the one chosen as true value. In the right panel this degeneracy is also solved with atmospheric data while for the values of δ and θ 23 chosen in the left panel a small sign degeneracy remains between the 4 and 5σ level. Notice that an "intrinsic degeneracy" [31] at δ π − δ true also shows up at the 5σ level when only the beam information is taken into account. As for the "sign" degeneracy, the atmospheric neutrino data is enough to lift it for the parameters chosen in the right panel while a small remnant is present in the left. In any case, both the "intrinsic" and the "sign" degeneracies appear at δ π − δ true , given the comparatively small matter effects for the setup, and their allowed regions are smaller or comparable to that of the true solution so that only the "octant"degeneracy plays a significant role in reducing the CP-discovery potential when atmospheric data is not exploited to lift it. In Fig. 6 we show how the significance with which the ESSνSB would be able to disfavour the wrong octant of θ 23 as a function of the true value of θ 23 (blue lines). As already anticipated in Section 2, this capability improves dramatically upon the inclusion of the atmospheric neutrino sample (red lines) and thus the potentially dangerous "octant" degeneracies are lifted. The curves are almost identical for both mass orderings and for the Zinkgruvan and Garpenberg baselines.
The significance with which the ESSνSB would be able to disfavour the wrong mass ordering is shown in Fig. 7, where dotted (solid) lines correspond to beam only data (beam and atmospheric data). The left (right) panels correspond to the 360 km (540 km) baseline and upper (lower) panels are for the scenario in which the true ordering is normal (inverted). As can be seen the ESSνSB beam data allows to disfavour the wrong mass ordering at around the 3σ (2σ) level for the 360 km (540 km) baseline for any value of δ and the octant. When the atmospheric data is added, the sensitivity to the wrong ordering is boosted to the 4-5σ level or even higher for the particular case of normal ordering and second octant of θ 23 (sin 2 θ 23 = 0.582 from Ref. [21]) for which the signal in atmospheric neutrinos is enhanced, as expected from Eq.(2.1). For normal ordering (upper panels) the inclusion of the atmospheric neutrino data also change the shape of the curve, in particular a larger increase in the significance is seen around δ = 0 than for other values. This is due to the solution of the octant degeneracy since, as can be seen in the middle panel of Fig. 2 or the first panel of Fig. 3, for δ = 0 and normal ordering the ellipse with opposite octant and ordering has a significant overlap. In Fig. 8 we analyze the precision with which the ESSνSB experiment would be able to measure the CP-violating phase δ. In this figure we assumed the currently preferred option of normal ordering and second octant of θ 23 . In the upper panels we show the improvement in the 1σ allowed region with which δ would be constrained by adding the atmospheric neutrino sample (solid lines) to the beam information alone (dotted lines). As can be seen, both for the 360 km (left panel) and 540 km baseline (right panel), the precision with which δ could be determined has a very pronounced shape. For CP violating values of δ around ±90 • , the 1σ uncertainty in the measurement peaks leading to the poorest precision, while for δ around 0 or 180 • the most precise measurements would be achieved.
As discussed in Ref. [69], this structure follows from the dependence of the oscillation probability on δ shown in Eq.(2.1). At an oscillation peak |∆m 2 31 |L/(4E) = (2n − 1)π/2 and thus mainly sin δ is probed. Since the derivative of sin δ vanishes at δ = ±90 • , the precision with which δ can be determined is worst close to these values. In order to constrain δ around δ = ±90 • , measurements away from the oscillation maxima to determine cos δ would instead be necessary. These off-peak measurements are easier at the Zinkgruvan 360 km baseline since the statistics is higher and also the beam is not exactly centered at the maximum, while they for the current best-fit parameters [21]. In the upper panels we show the comparison between the precision obtained with (solid lines) and without (dashed lines) the atmospheric sample for a running time of 5 years in each focusing. In the lower plots we show the dependence of the precision on the relative running time in each mode, where t ν (tν) corresponds to the time the experiment would run in neutrino (antineutrino) mode, combining atmospheric and beam datasets.
are very challenging at Garpenberg since very few events away from the oscillation peak are expected. This explains why the reconstructed sensitivities around δ = ±90 • are much worse in the right panel compared to the left. Moreover, the double-peak structure that can be seen for δ = −90 • for 540 km corresponds to the "intrinsic" degeneracies depicted in Fig. 5 that merge into one bigger allowed region. Since, as seen in Fig. 5, the addition of atmospheric data can lift these degeneracies, in the solid lines where this information was included the difference between the two baselines is significantly reduced. Conversely, for δ = 0 or 180 • the measurement on peak is what allows to determine δ and, since this is better covered at the longer 540 km baseline, the precision is slightly better there. This fact also translates into the better CP-discovery potential observed for the 540 km baseline in Fig. 4. Since the error in δ is smaller around CP-conserving values, the 540 km option could get closer to these values but still allow to claim the discovery of CP violation with high significance. In the lower panels of Fig. 8, the impact of changing the relative running times in positive focusing (neutrino mode) and negative focusing (antineutrino mode) is shown. Since off-peak measurements are required for δ = ±90 • , statistics are crucial and easier to accumulate in neutrino mode, since fluxes and cross sections are larger, and thus the best precision would be obtained by devoting longer periods of data taking to positive focusing. Conversely, around δ = 0 or 180 • the complementarity between the neutrino and antineutrino samples pays off and more even splits of the running time provide better sensitivity.
Since the ESSνSB would be a next-generation facility, its measurement strategy can profit from the previous hints by preceding oscillation experiments and adapt the splitting between neutrino and antineutrino modes depending on what value of δ data point to. If such a strategy is followed and the best splitting between neutrino and antineutrino modes is adopted for each value of δ, the precision presented in Fig. 9 would be obtained. If the mass ordering is confirmed to be normal and θ 23 lies in the second octant as present data prefer, the precision with which the ESSνSB facility would determine δ ranges from 16 • (13 • ) for δ ∼ −90 • to 6 • (7 • ) for δ ∼ 0 or δ ∼ 180 • for 540 km (360 km).
From Figs. 4 and 9 one can conclude that if the experiments preceding the ESSνSB do not find any evidence for CP-violation, the best option would be the 540 km baseline and a more or less even split of the neutrino and antineutrino running times. Indeed, this choice would minimize the errors with which δ would be determined around CP-conserving values and allow to increase the CP-discovery potential. On the other hand, if the previous set of experiments determine δ to be close to maximally CP-violating, then the best scenario for the ESSνSB would be the shorter 360 km baseline and increased neutrino run time to determine δ with the best precision possible.
In Fig. 10 we show the impact of individual systematic uncertainties on the fraction of values of δ for which CP violation could be discovered (∆χ 2 ≥ 25). The sources of uncertainty considered, summarized in Table 1, are the flux uncertainties for the signal (δφ S ) and background (δφ B ), the cross section systematic (δσ), the neutral current background (δN C B ), and the uncertainty on the ratio of the electron and muon flavour neutrino cross section (δσ e /σ µ ). The plot shows that the systematic uncertainties that most significantly   Table 1, red squares correspond to assuming that particular uncertainty to be 5 times larger and blue triangles to reducing the uncertainty by a factor of 5.
affect the performance of the ESSνSB are the ones related to the background components of the beam, since for these the determination at the near detector is more challenging. Namely, δφ B , δN C B as well as δσ e /σ µ since the only ν e present at the near detector that would allow to fix this parameter are those from the intrinsic background contamination of the beam. Among these, the strongest impact on the sensitivity is due to the cross section ratio since, not only it is difficult to constrain, but it is also most relevant to the signal at the far detector, which consists of ν e . Indeed, reducing or increasing this particular source of systematic error has the biggest impact on the physics reach. The impact is in any event limited, since the main bottleneck to the performance when observing at the second oscillation peak is statistics. In particular, a reduction of this systematic by a factor of 5 improves the CP fraction by ∼ 2% (no impact forν) while the same factor in the opposite direction worsens the sensitivity by ∼ 9% (∼ 4%).
The importance of these systematic errors in the physics reach is crucially dependent on the baseline of the experiment. In the left panel of Fig. 11 we show the fraction of all the possible values of δ for which it would be possible to rule out δ = 0 or δ = 180 • with a ∆χ 2 = 25 or higher significance. The upper blue line is for the more optimistic systematics from Table 1 and the lower red one for the more conservative values. As can be seen, the fraction of values of δ at which a 5σ discovery would be possible, peaks between 400 km and 700 km in both cases. But this peak is much more pronounced when the more conservative values are assumed for the systematic uncertainties. Indeed, for larger values of the systematics, the shorter baselines are strongly penalized since the dependence of the oscillation probability is subleading around the first peak and easily hidden by the systematics. Conversely, if very small systematic errors can be achieved, then the main limiting factor would be statistics and shorter baselines would perform better. Thus, by measuring at the second oscillation  Table 1. In the right panel we show the CP fraction for the Garpenberg (L = 540 km) and Zinkgruvan (L = 360 km) mines, assuming the current best fit values for the oscillation parameters and the "Optimistic" systematics for increasing total exposure. maximum the ESSνSB setup becomes much more resilient to sources of systematic errors unaccounted for than when observing only at the first peak. In the right panel of Fig. 11 we show how the fraction of values of δ for which CP violation would be discovered at the 5σ level by the ESSνSB beam and atmospheric data increases with the exposure. As expected from an observation at the second oscillation peak, statistics is the main factor controlling the final reach of the experiment. Indeed, for 5 years data taking the CP fraction is around 46%, by 10 years it increases to 62% and reaches 70% for 20 years of exposure. The slope only flattens significantly after 25 years.

Conclusions
In this paper we have performed an exhaustive analysis of the physics reach of the ESSνSB facility exploring its capability to determine all the presently unknown neutrino oscillation parameters such as the mass ordering and the octant of θ 23 but with a focus on the discovery of leptonic CP violation and a precision measurement of δ, which are the main declared goals of the experiment. For the first time we combined the atmospheric neutrino sample that would also be observed at the facility with the beam information and studied the complementarity between the two data sets. We studied how the physics reach of the facility could be optimized by exploring different baselines and focusing on the two candidate sites of Zinkgruvan at 360 km and Garpenberg at 540 km. We have also explored how the time split between neutrino and antineutrino modes can be exploited to improve the physics reach.
We conclude that the inclusion of the atmospheric data set can significantly increase the ESSνSB physics reach. Due to the peculiarities of observing the oscillation probability at the second oscillation maximum we find that this combination is particularly synergistic. The atmospheric neutrino sample not only significantly increases the sensitivity to the mass ordering, like for other similar facilities [29,30], but it is also very effective in improving the constraints on ∆m 2 31 and θ 23 and its octant. These measurements are especially challenging for the beam alone when sitting at the second maximum, given the low statistics, particularly in antineutrinos and in the ν µ disappearance channel. However, the determination of δ can be affected by correlations with θ 23 [34] and degeneracies with the wrong octant and thus the atmospheric information is also crucial to increase the CP discovery potential of the ESSνSB indirectly. We find this complementarity is somewhat more pronounced for the longer 540 km baseline since there the flux is more centered at the second oscillation peak and the statistics are smaller so it benefits more from the information gained from the atmospheric neutrino data.
Regarding the optimal baseline, we find the choice is rather dependent of the actual value of δ. For δ ∼ ±90 • a precise measurement needs events away from the oscillation maximum. In this sense the shorter 360 km baseline is better since the statistics for off-peak events are higher and this leads to a more precise measurement. Conversely, if δ is close to CP conserving values and the previous set of measurements have not been able to claim the discovery of CP-violation, the longer 540 km baseline would allow to cover a larger part of the parameter space. Indeed, after 10 years of data taking, the fraction of values of δ for which a 5σ discovery would be possible is 56% for Zinkgruvan and 62% for Garpenberg.
As for the splitting of the data taking time between neutrino and antineutrino modes, the optimal strategy also depends on the value of δ. This fact could be exploited since previous and present data at the time of the measurement should already show a strong preference for some part of the parameter space. Thus, the running strategy can be adapted to the situation optimizing the precision with which this measurement can be performed. In particular we find again that given the need of going beyond measurements at the peak for δ ∼ ±90 • , statistics is much more relevant and maximizing the time in neutrino mode translates to the best precision for these values. Conversely, close to CP-conserving values of δ, the information from events on-peak is most relevant and the complementarity between neutrino and antineutrino modes pays off so that a more even split of the running time would provide the best precision.
Finally we explored the possible bottlenecks for the physics reach of the facility exploring how it is affected by varying the values of the different systematic errors considered as well as the total exposure. As expected, the choice of observing the oscillation probability at its second maximum significantly reduces the impact of the systematic errors. We find that around the first oscillation peak the fraction of values of δ for which a 5σ discovery is possible is reduced by more than a factor 2 when considering the more conservative values of Table 1. On the other hand, at the second peak the reduction is only by a factor around 1.2. Among the different sources of systematic uncertainties considered, the most important is the possible difference in the ratio of the electron to muon neutrino cross sections. This uncertainty is difficult to constrain from near detector information since the flux is mainly composed of ν µ , but the far detector signal consists of ν e . Conversely, the observation at the second maximum considerably reduces the number of events and statistics play a much more relevant role. At the longer 540 km baseline, the fraction of values of δ allowing for a discovery would go from 47% to 62% and 70% for data taking periods of 5, 10, and 20 years, respectively. of the ESSνSB design study, in particular Monojit Ghosh, for comments on our manuscript. This work made extensive use of the HPC-Hydra cluster at IFT. This work is supported in part by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreements 674896-Elusives, 690575-InvisiblesPlus, and 777419-ESSnuSB, as well as by the COST Action CA15139 EuroNuNet. MB, EFM, and SR acknowledge support from the "Spanish Agencia Estatal de Investigación" (AEI) and the EU "Fondo Europeo de Desarrollo Regional" (FEDER) through the project FPA2016-78645-P; and the Spanish MINECO through the "Ramón y Cajal" programme and through the Centro de Excelencia Severo Ochoa Program under grant SEV-2016-0597. MB also acknowledges support from the Göran Gustafsson foundation.