Nonunitarity of the lepton mixing matrix at the European Spallation Source

If neutrinos get mass through the exchange of lepton mediators, as in seesaw schemes, the neutrino appearance probabilities in oscillation experiments are modified due to effective nonunitarity of the lepton mixing matrix. This also leads to new CP phases and an ambiguity in underpinning the ''conventional'' phase of the three-neutrino paradigm. We study the CP sensitivities of various setups based at the European spallation source neutrino super-beam (ESSnuSB) experiment in the presence of nonunitarity. We also examine its potential in constraining the associated new physics parameters. Moreover, we show how the combination of DUNE and ESSnuSB can help further improve the sensitivities on the nonunitarity parameters.


I. INTRODUCTION
The discovery of neutrino oscillations [1,2] has brought neutrinos to the center of particle physics. The current experimental data mainly converge into a consistent global picture in which the oscillation parameters are pretty well determined. However, three challenges still remain, namely, to determine the CP phase, the atmospheric octant and the ordering of the neutrino mass spectrum [3,4]. These will be the target of a number of future experiments, such as DUNE [5]. A fourth item must be added to this list, namely probing the robustness of the interpretation, such as testing the unitarity of the lepton mixing matrix. This is crucial because it undermines the efforts of underpining the CP phase δ CP [6,7].
This task is well justified also on theory grounds. Indeed, one of the most attractive ways to generate neutrino mass is through the mediation of heavy neutral leptons. While these emerge in many gauge extensions of the standard model, they can be postulated directly at the SU(3) c ⊗ SU(2) L ⊗ U(1) Y level, as the neutrino mass generation mediators.
This, in fact, provides the most general realization of the seesaw mechanism and many of its variants [8]. For generality here we focus exclusively on this case, namely, the standard SU(3) c ⊗ SU(2) L ⊗ U(1) Y seesaw mechanism. The resulting lepton mixing matrix is in general quite complex when compared with Cabibbo-Kobayashi-Maskawa (CKM) mixing. First, lepton mixing contains extra phases that can not be eliminated by field redefinitions [8] and are therefore physical [9], crucially affecting lepton number violating processes. However, they do not affect conventional oscillations, so we will ignore them in what follows. On the other hand the lepton mixing matrix must in general take into account the admixture of the heavy lepton seesaw mediators with the light active neutrinos [10]. These are usually neglected, as the smallness of neutrino masses indicated by neutrino experiments suggests a very high seesaw scale. ton number symmetry is restored at low, instead of high, values of the lepton number violating scale. The models are natural in t'Hooft sense, and lead to small, symmetryprotected neutrino masses. Such "low-scale" seesaw realizations include the inverse [17,18] as well as the linear seesaw mechanisms [19][20][21]. In all of these we expect potentially sizeable unitarity violation in the leptonic weak interaction. This paper is dedicated to probing such effects at the European Spallation Source neutrino Super-Beam (ESSnuSB) experiment. As far as we can tell, this is the first study of this kind for ESSnuSB. In addition, we also review the sensitivities of the DUNE experiment to the nonunitarity parameters and show how DUNE and ESSnuSB can play a complementary role to each other. Sensitivity studies of nonunitarity at other future neutrino facilities can be found in Refs. [7,[22][23][24][25][26][27][28][29].
We briefly describe the theoretical framework for unitarity violation in the charged current (CC) leptonic weak interaction in Sec. II, and the matter three-neutrino oscilation probabilities with and without unitarity violation in Sect. III. Next, in Sec. IV, we describe the experimental setups of interest, and also present the details of the simulation we have performed. Our results are given in Sec. V and include our calculated ESSnuSB and DUNE sensitivities to nonunitary (NU) neutrino mixing in Sec. V A, the CP violation discovery potential in the presence of unitarity violation is given in V B, and the CP reconstruction capabilities both for the standard phase as well as the seesaw phase of α 21 in V C. Finally we briefly summarize in Sec. VI.

II. THEORETICAL FRAMEWORK
In the standard 3 × 3 oscillation picture, the neutrino mixing matrix is described symmetrically by a product of three mixing matrices where each ω i,j describes an effective 2 × 2 complex rotation, characterized by a mixing angle and its phase. This symmetrical form complements the original description [8] by specifying the most convenient factor ordering. In explicit form, the standard 3 × 3 leptonic mixing matrix is given by where s ij = sin θ ij , c ij = cos θ ij , and ϕ ij is the corresponding phase. One sees the appearance of two extra physical phases with no counterpart in the quark sector: the so-called Majorana phases [8]. Note that the above parameterization of the neutrino mixing matrix is equivalent to the oscillation-sensitive part of the PDG form [30] with ϕ 13 − ϕ 12 − ϕ 23 ≡ δ CP [31] so that when ϕ 12 = ϕ 23 = 0 one has ϕ 13 = δ CP .
Apart from the presence of these new physical [9] phases, the leptonic CC interaction will in general also contain the mixing of neutral heavy leptons that mediate neutrino mass generation, as in the so-called type-I seesaw mechanism. These two facts make the mixing of massive neutrinos substantially richer in structure than that which describes the quark weak interactions [8]. As a result, in this general neutrino framework, U n×n can be expressed as the product of the new physics (NP) piece, times the Standard Model (SM) piece Thinking in terms of the seesaw mechanism it is convenient to express the full U matrix as four submatrices.
Notice that we have a block, N , relating the light neutrino sector with the three active neutrino flavors. Here V will be a (n − 3) × (n − 3) submatrix, while S and T will be, in general, rectangular matrices.
Clearly, in this general case, the full unitarity condition will take the form Therefore, the 3 × 3 matrix N describing the mixing of light neutrinos will no longer be unitary. One can show [33] that in the most general case N can be parametrized as where the diagonal α's are real and close to 1, while the off-diagonals are small but complex. Indeed, for any number of additional neutrino states, we will have for the diagonal entries of this matrix that with no sum over j. For small mixings, the non-diagonal entries are given as From this last expression, one can see that and similar equations for the other two non-diagonal terms. Using the triangle inequality, one can now derive the consistency relations 2 , The muon neutrino appearance probability will be given as where N is given in terms of the α's as in Eq. (6).
2 For a general derivation of this expression without the assumption of small mixing angles see [34].
For the case of vacuum oscillations, the parameters characterizing unitarity violation in the µ-e sector are α 11 , α 22 , and α 21 . In the presence of matter effects, the appearance probability could also involve the third neutrino type, since the charged and neutral current potential will modify the effective form of the matrix N . The charged current potential for the nonunitary case will be given by where G F is the Fermi constant and N e the number density of electrons in the medium. The matrix product will be written, in terms of the α's, as [7]: (13) The corresponding potential for neutral currents will be where the matrix product (N N † ) 2 , at leading order in the non-diagonal α's, takes the form [7]    Neglecting cubic terms in α 21 , sin θ 13 , and ∆m 2 21 , one finds that, in the vacuum case limit, the main contribution to the conversion probability will be given by where P 3×3 µe denotes the usual three-neutrino conversion probability, P 3×3 µe = 4 cos 2 θ 12 cos 2 θ 23 sin 2 θ 12 sin 2 ∆m 2 21 L 4E ν + cos 2 θ 13 sin 2 θ 13 sin 2 θ 23 sin 2 ∆m 2 31 L 4E ν + sin 2θ 12 sin θ 13 sin 2θ 23 sin where P I µe is the interference term P I µe = −2 sin(2θ 13 ) sin θ 23 sin with φ 21 = arg(α 21 ).

III. THREE NEUTRINO OSCILLATION PROBABILITIES
Before coming to our numerical results, in this section we discuss the behaviour of the appearance (anti)neutrino probabilities. To this end, we show in the left (right) panel of Fig. 1 the ν µ → ν e (ν µ →ν e ) oscillation probabilities as a function of the neutrino energy. The upper, middle, and lower panels correspond to 540 km, 360 km, and 200 km baselines, respectively. We show the conversion probability in the standard unitary framework as a solid line, while the dashed one represents the nonunitary case. For the unitary case, we consider the values of the neutrino oscillation parameters given in Table I with δ CP = −90 • . For the nonunitary case, besides these values, we fix the nonunitary α-parameters to be α 11 = 0.97, α 22 = 0.99, α 33 = 1, |α 21 | = 0.02, φ 21 = 90 • , |α 31 | = 0, and |α 32 | = 0. In Fig. 1 we also include in green color (arbitrary units) the unoscillated (anti)neutrino ν µ flux times the (anti)ν e -nucleus cross section [35,36]. This provides a bird-eye view of the real analysis at the event level. The main region of interest extends from 0.1 GeV to 0.7 GeV, with a peak around 0.30 GeV. One can see the relative location of this peak with respect to the second oscillation maximum for the three baselines of 540 km and 360 km and the 200 km case, see also Table II. The true expected event number for the appearance signal must, of course, take into account the convolution of the cross section, the appearance probability and the unoscillated neutrino energy spectrum at the detector, which depends on its distance to the source. Thus, the total unoscillated neutrino flux for the 200 km case will be approximately three times that for a 360 km baseline, and seven times the flux of the 540 km case. As a result, the shortest the baseline, the higher the expected number of events, as we will see in Fig. 3.

IV. EXPERIMENTAL SETUP AND SIMULATION
In this section, we briefly discuss the experimental specifications of the ESSnuSB and DUNE setups used in this work, followed by a description of our simulation procedure.

A. Experimental setup options
The ESSnuSB project is a proposed accelerator neutrino experiment sourced at Lund (Sweden), where the ESS linac facility is currently under construction. The original ESSnuSB proposal was to use of a very intense proton beam of 2 GeV energy and an average beam power of 5 MW, resulting in 2.7 × 10 23 protons on target (POT) per year (208 effective days) [37][38][39][40]. Here we will adopt this configuration. It is expected that the future linac upgrade can increase the proton energy up to 3.6 GeV. The neutrino and antineutrino fluxes arising from the 2 GeV proton beam peak around 0.25 GeV [41]. These (anti)neutrinos will be detected by a 500 kton fiducial mass Water Cherenkov detector similar to the MEMPHYS project [42,43]. Since the baseline of the far detector has not been finalized yet, we have considered in this work three possible baselines [37], which are 200 km, 360 km, and 540 km respectively. It has been shown in [37] that if the detector is placed in any of the existing mines in between 200 km to 600 km from the ESSnuSB site Lund, a 3σ evidence of CP violation could be achieved for 60% coverage of the full δ CP range. Our simulation matches the event numbers of Table 3 and all other results given in [37]. In all the numerical results presented here, we have assumed 2 years of neutrino and 8 years of antineutrino running with an optimistic assumption of uncorrelated 5% signal normalization and 10% background normalization error for both neutrino and antineutrino appearance and disappearance channels, respectively. For more details about the accelerator facility, beamline design, detector and baseline positions of this setup, see [37]. Note that, while working on this paper, an updated analysis from the collaboration has come out [44]    ity is expected to start taking data around the year 2030.
DUNE is a future long-baseline accelerator-based neutrino experiment with a baseline of 1300 km from the source at Fermilab to the far detector placed deep underground at the Sanford Laboratory site in South Dakota. DUNE will use a 40 kton LArTPC detector and a 120 GeV proton beam with 1.2 MW beam power resulting in 1.1 × 10 21 POT/year. For the numerical simulations, we have followed the experimental configurations provided by the collaboration in the Technical Design Report (TDR) [45,46], assuming equal runtime of 3.5 years in neutrino and antineutrino mode, which results in 336 kton-MW-year exposure for the TDR setup. More details on the systematic errors, efficiencies and energy resolutions can be found in Refs. [46,47].

B. Simulation procedure
In order to assess the statistical sensitivity of the ESS-nuSB facility to neutrino oscillations, we have made use of the built-in χ 2 function of the GLoBES package [48,49], which incorporates the systematic errors through the pull terms [50]. To perform the nonunitarity analysis we have used the modified version of [51]. The total χ 2 is a sum of all the contributions coming from different channels, Unless stated otherwise, the benchmark choices of the standard three-neutrino unitary oscillation parameters and their marginalization status in our analysis are given in Table I. Our benchmark choices closely follow the current global fit analysis [3]. Following the same analysis, we have adopted a 1% uncertainty on the atmospheric mass-squared splitting ∆m 2 31 and a 3.2% uncertainty on the reactor mixing angle sin 2 θ 13 . We have freely marginalized over the atmospheric parameter sin 2 θ 23 from 0.35 to 0.65.
For the case of ESSnuSB, we have considered a lineaveraged constant matter density ρ = 2.8 g/cm 3 following the PREM profile [52,53]. For DUNE, we have also assumed the same matter density but with a 5% uncertainty due to the longer baseline. For definiteness, we have assumed the currently preferred case of normal mass ordering (NO) throughout all of our analyses. Whenever appropriate, we have also marginalized over the NU parameters, along with their associated CP phases, implementing the current 3σ bounds shown in Table III. These come essentially from short-baseline oscillation searches such as from NOMAD [54,55] and CHORUS [56,57] and the long-baseline experiments T2K [58], NOvA [59] and, most importantly, MINOS/MINOS+ [60].

V. RESULTS
In this section we discuss in detail the numerical findings of our analyses, where we explore the sensitivity of the ESSnuSB facility to the nonunitary neutrino mixing, as well as the impact of nonunitarity on the measurement of the standard three-neutrino oscillation parame-ters, with emphasis on the CP-violating phase, δ CP .

A. Probing nonunitary neutrino mixing at ESSnuSB
We start our discussion from Fig. 2, where we show the ESSnuSB sensitivity to nonunitarity in the (|α 21 |, δ CP ) plane. Left, middle and right panels show the results for 540, 360, and 200 km baselines, respectively. The light, medium, and dark green contours in each panel correspond to the 1σ, 2σ, and 3σ allowed regions for 2 degrees of freedom (d.o.f.) i.e., ∆χ 2 = 2.3, 6.18, and 11.83, respectively. In this figure we have assumed the standard unitary framework as the true hypothesis, and then we have fitted the nonunitary hypothesis against it. Normal mass ordering has been assumed all over the analyses. The true data have been generated assuming the benchmark choices of the standard unitary oscillation parameters in Table I with δ CP (true) = −90 • , while for the reconstruction we have fixed the solar oscillation parameters and marginalized over θ 13 , θ 23 , and ∆m 2 31 . In addition, we have also marginalized over the NU parameters α 11 , α 22 , and α 33 within their allowed 3σ ranges as given in Table III. We have also freely varied φ 21 from −π to +π. We have assumed zero values of the other non-diagonal NU parameters, |α 31 | and |α 32 |. The red vertical lines in each panel indicate the current 3σ upper limit on |α 21 | from neutrino data. One sees that the 200 km baseline gives the best 1σ, 2σ, and 3σ sensitivities on |α 21 |, in comparison to the 360 km and 540 km baseline option. Quantitatively, the attainable upper limits for |α 21 | at 1σ are 0.044, 0.024, and 0.02 for 540 km, 360 km, and 200 km baselines, respectively. Conversely, we have checked that, there is basically no sensitivity to the new CP phase φ 21 for any of the baselines. On the other hand, the measurement on δ CP is not affected much by the presence of nonunitarity. One can also see that the uncertainty on the measurement of δ CP is the lowest for 200 km and the highest for 540 km.   Fig. 3, see text for more details.
To understand better the previous result one can select a point allowed at 1σ, 2σ, and 3σ for a 540 km, 360 km, and 200 km baseline, respectively. We show in Fig. 3 the appearance event spectra for the neutrino (upper-left panel) and antineutrino mode (upper-right panel). For each of the three different baselines we choose a benchmark point (magenta star in Fig. 2) given by the NU parameter |α 21 | = 0.04 and the standard CP phase δ CP = −90 • . For this benchmark choice, the other NU parameters arising from the χ 2 marginalization take the values shown in Table IV. Note that the optimized values of the standard oscillation parameters are not shown here explicitly. We use all of them as input to generate the event spectra in Fig. 3. These are shown as dashed lines, while the unitary scenario is represented with solid lines. In order to better compare the standard unitary (UN) with the NU case, we show in the lower panels of Fig. 3 the absolute differences of the number of events divided by the statistical uncertainty in each bin i.e., ∆ = |N UN − N NU |/ √ N UN . This provides a crude measurement of the statistical significance of our unitarity test. One sees that the best sensitivities to the standard parameter δ CP and the NU parameter |α 21 | are achieved for the 200 km baseline, as already shown in Fig. 2. Notice also that the 360 km baseline does somewhat better than the 540 km baseline. This is also reflected in the left panels of Fig. 4.
In Fig. 4, we show the one-dimensional projection of the ∆χ 2 as a function of the NU parameters |α 21 | (left panel), α 11 (middle panel), and α 22 (right panel). The red, blue, and green curves in the upper panels show the sensitivities corresponding to the three ESSnuSB baselines, 540 km, 360 km, and 200 km, respectively, while the black curve corresponds to the DUNE sensitivity. In the lower panels, the same color code is used to represent the sensitivities for the combined analysis of DUNE and the ESSnuSB for baselines of 540, 360, and 200 km.
As before, this figure is obtained assuming the standard unitary framework as the true hypothesis, and testing the NU framework against it, i.e. ∆χ 2 = χ 2 NU − χ 2 UN . The true values of the oscillation parameters are taken from Table I, fixing the two solar parameters θ 12 and ∆m 2 21 at their benchmark choices, and marginalizing over θ 13 , θ 23 , and ∆m 2 31 . Since the exact value of the standard CP phase δ CP is currently not accurately known, we have marginalized over its true and test values within its full range.
For the sensitivity to |α 21 |, shown in the left panels of Fig. 4, we have freely varied the NU parameters α 11 , α 22 , and α 33 within their allowed 3σ ranges as given in Table III and the phase φ 21 from −π to +π. The other non-diagonal NU-parameters |α 31 | and |α 32 | have been set to zero. From the upper left panel we can see that the best sensitivity for the ESSnuSB comes from the 200 km baseline, followed by the 360 km and 540 km baselines. The expected 3σ upper limits on |α 21 | corresponding to 540 km, 360 km, and 200 km baselines are 0.16, 0.075, and 0.04, respectively. These limits would be independent and complementary to those given in Table III and illustrate the ESSnuSB potential in probing new physics. For the case of DUNE, the upper limit at 3σ would be 0.046. One sees the sensitivities expected at ESSnuSB are competitive and, for the smaller baseline, it performs slightly better than DUNE [7]. In the lower left panel of Fig. 4, we show how the combined sensitivities of DUNE and ESSnuSB, improve with respect to the individual sensitivities. The 3σ upper limits corresponding to the DUNE+ESS540, DUNE+ESS360, and DUNE+ESS200 combinations would be 0.04, 0.03, and 0.022, respectively. It is encouraging to see that the combination of DUNE and the 200 km baseline of ESSnuSB would be able to set a constraint on |α 21 | which is better than the current 3σ upper bound given in Table III.
The middle panel of Fig. 4 shows the sensitivities to the diagonal NU parameter α 11 . The analysis method is very similar to the previous one as far as the standard unitarity-parameters are concerned. However, for the NU parameters, we marginalize over the parameters α 22 , |α 21 |, and its associated new CP-phase φ 21 . As before, the off-diagonal parameters α 31 and α 32 have been set to zero and α 33 was set to unity. One can see that the 3σ sensitivities on α 11 corresponding to the 540, 360, and 200 km baselines are lower than for DUNE, i.e. 0.87 versus 0.95. From the lower middle panel one sees that the combined sensitivities corresponding to DUNE+ESS540, DUNE+ESS360, and DUNE+ESS200 are very similar and approximately equal to 0.955, which would somewhat improve the current sensitivity on α 11 . The sensitivities on α 22 are shown in the right panels of Fig. 4. In this case, we follow exactly the same steps as in the α 11 analysis but marginalizing over α 11 . Our results show that the best performance is expected from DUNE, followed by the ESSnuSB with baselines 360 km, 540 km, and 200 km. The expected sensitivities for the individual setups do not improve over the current lower bound. However, one sees that the combined sensitivity of DUNE and any ESSnuSB baseline gives a slightly lower sensitivity compared to the current 3σ lower bound. One concludes that all these sensitivities, specially those coming from the combined analyses, are encouraging and will play an important role especially when combined with the already available data. For the sake of convenience  we quote all our bounds as Table V. Note also that these NU parameter sensitivities at the ESSnuSB may be complemented by additional information coming from future experiments such as T2HK or JUNO [61,62].

B. CP violation discovery potential
Now we turn back to the "conventional" CP violation discovery potential within our generalized nonunitary framework. The CP violation (CPV) discovery potential of the ESSnuSB and DUNE setup is summarized in Figs. 5, 6 and 7. Our results are given both for the unitary and the nonunitary framework. As in the standard δ CP sensitivity study [5], the CP-violating hypothesis is tested against a CP-conserving scenario through [7] This way, we obtain the significance with which one can reject the test hypothesis of no CP violation. We have assumed five nonzero NU parameters: the three diagonal ones (marginalized over the true and test  values), plus one non-diagonal parameter, either |α 21 | (left panels) or |α 31 | (right panels), with the associated complex CP phase (φ 21 or φ 31 ), also marginalized. Concerning the standard three-neutrino parameters (both in the unitary and nonunitary scenarios), we have marginalized over the two mixing angles θ 13 and θ 23 and the "atmospheric" mass-squared splitting ∆m 2 31 . The left (right) column of Fig. 5 represents the CPV discovery potential of ESSnuSB for the three baseline choices and three different values of the NU parameter |α 21 | (|α 31 |). The black solid line in each plot corresponds to the standard three-neutrino unitary framework. Upper, middle, and lower panels represent the results obtained for 540, 360, and 200 km baselines, respectively. The red, green, and blue dashed curves in the left plots represent the CPV discovery sensitivities in the NU framework corresponding to different values of |α 21 |, 0.01, 0.02, and 0.03, respectively (the latter, relatively large value, is taken mainly for comparison). One sees that the CPV sensitivity in the standard unitary neutrino oscillation framework always lies around 8σ for δ CP (true) = ±90 • for all three baselines, see the solid black line in all the panels. This agrees with the results presented in Refs. [37,44,63]. All baselines have more or less similar sensitivities, except for the fact that the δ CP range over which CPV can be established for 540 km and 360 km is slightly bigger than for 200 km. This fact is also confirmed in Ref. [37]. However, we will see the merits of the 200 km baseline in what follows. As far as the NU framework is concerned, two of our benchmark values, |α 21 | = 0.01, and 0.02 lie within the current 3σ limit, whereas |α 21 | = 0.03 lies slightly outside the current allowed limit. 4 .
We stress that the CPV discovery sensitivity is degraded with respect to the unitary case. This is expected, due to the presence of new phases associated to unitarity violation [6]. Clearly, the CPV discovery sensitivity decreases with the increasing values of |α 21 |, specially for 540 km, leading to a minimum sensitivity of 4.6σ for δ CP (true) = ±90 • . The deterioration of the CPV sensitivity is smaller for 360 km and, with a minimum sensitivity of around 5.7σ for δ CP (true) = ±90 • . For the 200 km baseline, the deterioration further reduces, with a minimum sensitivity of 6.1σ for all three benchmark choices, at δ CP (true) = ±90 • . All in all, one sees that the degrading in CP sensitivity is not as large as one might expect, showing the robustness of the oscillation picture with respect to unitarity violation. The best sensitivities to δ CP and the nonunitary parameter |α 21 | are achieved for a 200 km baseline.
We now turn to the right panels of Fig. 5. There, we repeated the same analysis for the non-diagonal parameter |α 31 | and its associated CP phase, φ 31 . The three benchmark choices considered for |α 31 | are 0.03, 0.07, and, 0.10 5 . In this case, one finds a mild deterioration of the CPV sensitivity in comparison to the unitary framework for all baselines, so the impact of |α 31 | is not significant for 540 km, and negligible for 360 km and 200 km baselines. This is due to the fact that |α 31 | does not appear in the vacuum appearance probability, Eq. (16), and also because of the lower matter effects for ESSnuSB with respect to DUNE [7,64]. As expected, one finds a negligible impact on the CPV sensitivity arising from unitarity violation in this case. Note, however that this will not be the case for DUNE and, as a result, the combined analysis with ESSnuSB will be useful to recover the maximum CPV sensitivity in this case.
Next, in Fig. 6, we show the CPV sensitivities for DUNE. As before, the black solid lines in each panel correspond to the CPV sensitivity in the standard unitary framework, showing that DUNE can establish CPV discovery above 5σ for δ CP (true) = ±90 • . The red, green, and blue dashed curves in the left panel represent the CPV discovery sensitivities in the NU framework corresponding to different values of |α 21 |: 0.01, 0.02, and 0.03, respectively. Similarly, in the right panel, these lines correspond to the CPV sensitivities for |α 31 | equal to 0.03, 0.07, and 0.10. The simulation method followed in both the panels is the same as for Fig. 5. From the left panel, we see that the CPV sensitivity in the presence of NU is degraded with respect to the unitary framework, according to the value of |α 21 |. As a result, for the maximal values of δ CP (true) one can only ensure a minimum CPV sensitivity of 3.6σ. Similarly, from the right panel, we can see that the CPV sensitivity decreases with increasing |α 31 | with a minimum sensitivity of 3.6σ for δ CP (true) = ±90 • .
Comparing Figs. 5 and 6 one sees that, for the conventional unitary scenario the CPV sensitivity of ESSnuSB is better than that of DUNE. Also the deterioration due to NU effects is smaller than for DUNE, where the matter effect is large. Hence, in this case, ESSnuSB will play a complementary role to DUNE in achieving a robust CP violation sensitivity in the presence of nonunitarity.
So far we have analyzed the effect of the non-diagonal NU parameters on the CPV sensitivity of the experimental setups under consideration. Now we will investigate the role of the diagonal nonunitary parameters, α 11 , α 22 , and α 33 . The corresponding CPV sensitivities in this case for ESSnuSB and DUNE are shown in the four panels of Fig. 7. The results for the standard unitary case (solid black line) for both the ESSnuSB and DUNE have been already described in Fig. 5 and Fig. 6. The red, green, and blue dashed lines in each panel correspond to the CPV sensitivities for fixed values of α 11 = 0.94, α 22 = 0.99, and α 33 = 0.90, respectively, in the data as well as in the theory. The procedure we have followed for this analysis is slightly different from what we have used for the non-diagonal parameters. In order to evaluate the CPV sensitivity in the NU framework for α 11 = 0.94, we have marginalized over α 22 , |α 21 |, and φ 21 within their allowed 3σ ranges. All other off-diagonal parameters are assumed to be zero, and we have also fixed α 33 = 1 in the data as well as in the theory. The analysis method for α 22 is analogous to the previous one, where we have marginalized over α 11 , |α 21 |, and φ 21 . On the other hand for α 33 , we marginalized over α 11 , |α 31 |, and φ 31 , fixing α 22 = 1.
One sees that ESSnuSB performs better than DUNE for the overall CPV discovery sensitivity. Of course, in all cases these sensitivities are degraded substantially with respect to those of the unitary case. In fact, the CPV discovery for DUNE falls below 4σ for any value of δ CP (true). As seen in Fig. 7 the best result is obtained for the 200 km baseline, followed by 360 km of ESSnuSB.
In short, the ESSnuSB will play a crucial role in securing a robust CPV discovery sensitivity in the presence of unitarity deviations close to the current upper bound.

C. CP reconstruction
In our simplest scenario there are two relevant CP phases, the standard three-neutrino Dirac phase δ CP and the phase φ 21 associated to nonunitarity 6 . One can therefore have four "CP conserving" cases, when either of them equals 0 or π. Likewise, four cases in which one has "maximal" CP violation, defined by having the modulus of any of them equal to π/2. In this section, we discuss how well the European spallation source setups can reconstruct the standard CP phase δ CP as well as the nonunitarity phase φ 21 for a few selected benchmarks.
The results of our analysis are shown in Fig. 8. The upper two panels correspond to the two CP conserving cases (0, 0) and (π, π), and the lower two panels to the two CP violating scenarios (−π/2, −π/2) and (π/2, π/2). The red, green, and cyan contours in each panel correspond to 540, 360, and 200 km baselines of the ESSnuSB experiment, respectively, whereas the orange contours represent the sensitivity expected in DUNE. All contours correspond to the 2σ level for 2 d.o.f. For this analysis, we have fixed |α 21 | = 0.02, which lies within the current 3σ allowed boundary. All the other off-diagonal NU parameters have been kept fixed to zero. We have marginalized over the mixing angles θ 13 and θ 23 and the atmospheric mass-squared splitting, ∆m 2 31 , as well as over the true and test values of the NU parameters α 11 , α 22 , α 33 within their allowed ranges.
We have checked that the expected 1σ uncertainties on δ CP (φ 21 ) for the CP-violating scenarios are 16 • (100 • ) for 540 km, and 13 • (70 • ) for 360 km. For the CP-conserving benchmarks, the uncertainties on δ CP are 10 • for 360 km, and 12 • for 540 km respectively, whereas no sensitivity on φ 21 is found for these two baselines. On the other hand, for the 200 km baseline, the typical 1σ level uncertainty on δ CP (φ 21 ) is 10 • (40 • ). One sees that the best performance would be obtained for the shortest baseline, 200 km. For comparison, we have also projected the sensitivity of the DUNE experiment in the same plot. For DUNE, the typical 1σ level uncertainty on the reconstructed CP phase δ CP (φ 21 ) is 21 • (45 • ).
In short, one can see that the δ CP reconstruction capability of ESSnuSB does not get too much impaired by the presence of unitarity violation, while the sensitivity to the NU phase φ 21 is competitive with that in DUNE and even better for 200 km baseline. Certainly, a future combined analysis of the actual results of these experiments would improve the standard and new CP-phase reconstruction sensitivities.

VI. SUMMARY AND OUTLOOK
Here we have explored the physics potential of the proposed European Spallation Source facility in the presence of nonunitarity of the lepton mixing matrix, as generally expected within the seesaw paradigm. We have also explored the DUNE physics potential in this framework as well as the combined result of both DUNE and ESSnuSB experiments. First, we presented in detail the theory framework of neutrino oscillations with effective nonunitary neutrino mixing, discussing in Fig. 1 the resulting neutrino and antineutrino appearance oscillation probabilities. Throughout the paper we have assumed normal mass ordering, and considered three reference baseline choices of 540, 360, and 200 km for the ESSnuSB setup. In Fig. 2 we have presented the sensitivity contours in the (|α 21 |, δ CP ) plane. The promising results for the 200 km baseline were understood in terms of the expected ν µ → ν e andν µ →ν e appearance event spectra. These are given in Fig. 3 as a function of the reconstructed neutrino energy. We found encouraging ESSnuSB and DUNE sensitivities for the off-diagonal NU parameter |α 21 |, and the two diagonal NU parameters α 11 and α 22 , as seen in Fig. 4. One also sees how the combined analysis of DUNE and ESSnuSB can help improving the sensitivity on the above mentioned NU parameters beyond the current 3σ bounds 7 .
More remarkable, perhaps, is the ESSnuSB CPV discovery potential, which is better than that of DUNE, as illustrated in Figs. 5, 6, and 7. One appreciates also a milder degrading of the ESSnuSB CP violation sensitivities with respect to the standard unitary mixing scenario in comparison with DUNE, where the effect is stronger. ESSnuSB would therefore contribute to establishing the robustness of CP determination against small departures from unitarity arising, say, from the seesaw mechanism. Finally, we have also obtained a promising CP reconstruction potential, both for the standard CP phase of the three-neutrino paradigm, as well as for the phase associated to nonunitarity, see Fig. 8. Altogether, within the generalized nonunitary neutrino mixing framework, we have found that the proposed ESSnuSB facility is competitive and complementary to DUNE, not only for leptonic CP violation studies, but also for probing new physics parameters associated to unitarity violation.