Heavy neutrino mixing in the T2HK, the T2HKK and an extension of the T2HK with a detector at Oki Islands

We study the discovery potential for the mixing of heavy isospin-singlet neutrinos in extensions of the Tokai-to-Kamioka (T2K) experiment: The Tokai-to-Hyper-Kamiokande (T2HK), the Tokai-to-Hyper-Kamiokande-to-Korea (T2HKK), and a plan of adding a new detector at Oki Islands to the T2HK. We parametrize the mixing of heavy neutrinos in terms of a non-unitary mixing matrix acting on active flavors. It is shown that in the T2HK and T2HKK, the sensitivity to heavy neutrino mixing deteriorates for some values of $CP$-violating phases in the standard and the non-unitary mixing matrices, but the deterioration is drastically mitigated if a detector at Oki Islands is added to the T2HK. We also consider the feasibility of measuring the mass hierarchy and the standard $CP$-violating phase $\delta_{CP}$ in the presence of heavy neutrino mixing, by fitting experimental data with the standard oscillation parameters only, i.e., under the assumption of unitary mixing. It is revealed that such measurement can be performed to some accuracy in the T2HKK and the extension of the T2HK with a detector at Oki Islands, owing to the fact that data with largely varying $L/E$ are collected in these experiments, while it cannot be in the T2HK.


Introduction
It is a viable possibility that isospin-singlet neutrinos mix with active flavors ν e , ν µ , ν τ after electroweak symmetry breaking. Notably, such a mixing is inherent in the Type-I seesaw model [1], and although the mixing angle is undetectably small in most of the parameter space, it can be sizable when there is a cancellation in the seesaw mass. To see this, suppose, for example, there are three isospin-singlet neutrinos with a Majorana mass M N , having a Dirac mass term m D with active neutrinos. Since the tiny neutrino mass is derived as −m D M −1 N m T D , one can express m D as [2] where U P M N S denotes the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix [3,4], m 1 , m 2 , m 3 are the tiny neutrino mass, and R 3 is a 3-dimensional rotation matrix with complexvalued rotation angles α, β, γ. Unless m 1 , m 2 , m 3 are degenerate, one can arbitrarily enhance m D and the mixing angles between isospin-singlet and active neutrinos by taking large imaginary values for α, β, γ. In this case, a cancellation is taking place in the components of −m T D M −1 N m D . Large mixing between isospin-singlet and active neutrinos is also realized in extensions of the Type-I seesaw model, such as the inverse seesaw model [5].
In the presence of isospin-singlet neutrino mixing, transitions among active flavors require for their consistent description extra parameters beyond the standard oscillation parameters, the three mixing angles θ 12 , θ 23 , θ 13 , one CP -violating phase δ CP and two mass differences ∆m 2 12 , ∆m 2 13 . Tension in fitting neutrino oscillation data with the standard parameters would thus be evidence for isospin-singlet neutrino mixing. If isospin-singlet neutrinos are "heavy", by which we signify that their mass is much above the neutrino beam energy, they do not contribute to oscillations and simply alter the mixing matrix U P M N S into a non-unitary matrix in the following fashion: where α ii 's are real and α ij 's (i = j) are complex numbers.
One purpose of this paper is to estimate the largest possible significance of heavy isospinsinglet neutrino mixing in the Tokai-to-Hyper-Kamiokande (T2HK) [6] and the Tokai-to-Hyper-Kamiokande-to-Korea (T2HKK) [7] (for early proposals, see Refs. [8,9]) experiments, as well as an experiment in which a new detector at Oki Islands is added to the T2HK (for early proposals on a detector at Oki, see Refs. [10,9]). Emphasis is on how the significance varies with CP -violating phases in the PMNS and the non-unitary mixing matrices and to what extent the sensitivity deteriorates with these phases. The secondary goal of this paper is to pursue the possibility that the neutrino mass hierarchy and the standard CP -violating phase in the PMNS matrix δ CP are measured in the presence of heavy neutrino mixing by fitting neutrino data with the standard oscillation parameters only, namely, without assuming heavy neutrino mixing. This is achievable if data allow us to discern contributions from the standard oscillation and those from the non-unitary mixing matrix. Such measurement is of practical importance when a hint of heavy neutrino mixing were discovered but were not conclusive. The mass hierarchy and δ CP measurement in the presence of heavy neutrino mixing has been investigated in a general context [11] and for the T2K, T2HK, NOνA and DUNE experiments [12,13,14,15].
The T2HK and the T2HKK are proposed extensions of the on-going T2K experiment [17], utilizing an off-axis ν µ +ν µ beam delivered from J-PARC. In the T2HK, there are two water Cerenkov detectors, each with 187 kton fiducial volume, located solely at Kamioka, where the baseline length is L = 295 km and the off-axis-angle is 2.5 • . In one plan of the T2HKK, one water Cerenkov detector with 187 kton fiducial volume is placed at Kamioka and the same detector is also at a site in Korea where the baseline length is L ≃ 1000 km and the off-axisangle is 1 • . We also consider a plan of placing, in addition to the T2HK, a small 10 kton water Cerenkov detector at Oki Islands, where the baseline length is L = 653 km and the off-axis-angle is 1.0 • . Merits of the above long baseline experiments for the search for heavy neutrino mixing are as follows: (i) In disappearance experiments, the effects of heavy neutrino mixing as parametrized by Eq. (2) cancel between measurements at the near and far detectors, while the T2K and its extensions are primarily appearance experiments. (ii) ν µ → ν e (ν µ →ν e ) transition amplitude arising from the non-unitary mixing matrix Eq. (2) interferes with the standard oscillation amplitude, enhancing the impact of non-unitary mixing on the transition probability. The merit (ii), however, bears two subtleties: • CP -violating phases in the PMNS and the non-unitary mixing matrices possibly lead to suppression of the interference.
• Non-unitary mixing α 21 in Eq. (2) can mimic matter effects and the effect of the standard phase δ CP , which distorts the measurement of the mass hierarchy and δ CP .
These issues become transparent by expressing the transition amplitude in the leading order of sin θ 13 , ∆m 2 21 , α 21 and matter effects as e −iδ CP sin θ 13 sin θ 23 + ∆m 2 21 L 4E sin(2θ 12 ) cos θ 23 sin(2θ 12 ) cos θ 23 e iδ CP sin θ 13 sin θ 23 + 1 2 where E denotes the neutrino energy and L the baseline length, and V cc is the potential due to charged current interaction inside matter.
We notice that the following interference terms, sin θ 13 sin θ 23 Im α 21 e i∆m 2 31 L/(4E)−iδ CP − 2 ∆m 2 21 L 4E sin(2θ 12 ) cos θ 23 Im(α 21 ), (6) are simultaneously suppressed at the first oscillation peak, ∆m 2 31 L/(4E) = π/2, if Arg(α 21 ) ≃ δ CP ±π/2 holds accidentally (note the relation ∆m 2 21 L/(4E) ≃ 0.1·2 sin θ 13 at the first oscillation peak). This suppression is mitigated if the energy distribution is not concentrated at the first oscillation peak. We thus infer that when Arg(α 21 ) ≃ δ CP ± π/2 holds, the T2HKK with a Korea detector placed at 1 • off-axis angle exhibits better sensitivity to α 21 than the T2HK, and that the extension of the T2HK with an Oki detector gives qualitative improvement. Of course, these anticipations depend on statistics, namely, competition with the loss of statistics in the T2HKK, or the size of an Oki detector.
We observe in Eqs. (3), (4) that when δ CP ≃ 0, π, the real part of non-unitary mixing Re(α 21 ) plays exactly the same role as matter effects, which afflicts the measurement of the mass hierarchy because it is probed through interference with matter effects. However, when measurements with different baselines and hence different V cc L's are combined, one resolves degeneracy between the non-unitary mixing and matter effects, thereby resurrecting sensitivity to the mass hierarchy. We thus expect that the mass hierarchy measurement in the T2HK is subject to strong influence from heavy neutrino mixing, whereas the T2HKK maintains sensitivity to the mass hierarchy even with heavy neutrino mixing. Also, the addition of a small detector at Oki to the T2HK leads to qualitative improvement in the mass hierarchy measurement.
Eqs. (3), (4) tell us that when δ CP ≃Arg(α 21 ), non-unitary mixing α 21 and the first term in the first line of Eqs. (3),(4), which is proportional to e iδ CP , interfere maximally. In such a case, the presence of α 21 affects the measurement of δ CP and causes a discrepancy between the true δ CP and the value measured by fitting data with the standard oscillation parameters only.
This effect is mitigated if the baseline is longer, because more information about the energy dependence of the standard oscillation amplitude helps distinguish it from non-unitary mixing α 21 . We thus expect that the discrepancy is smaller in the T2HKK than in the T2HK.
In the main body of this paper, we demonstrate, through a simulation, the following: (a) The sensitivity to heavy neutrino mixing is suppressed for Arg(α 21 ) ≃ δ CP ± π/2. Contrary to our anticipations, the T2HK shows better sensitivity than the T2HKK for any values of CP -violating phases, simply because of higher statistics. Also, adding a small detector at Oki to the T2HK does not make qualitative difference.
(b) In the presence of heavy neutrino mixing, the wrong mass hierarchy can be favored in the T2HK. The T2HKK maintains sensitivity to the mass hierarchy and measures it correctly despite heavy neutrino mixing. The extension of the T2HK with a small detector at Oki exhibits revived sensitivity to the mass hierarchy, even if statistics at Oki is quite subdominant compared to that at Kamioka.
(c) In the presence of heavy neutrino mixing, the true δ CP and the value obtained by fitting data with the standard oscillation parameters can have a sizable discrepancy that is of the same order as 1σ resolution of the δ CP measurement. The discrepancy tends to be larger in the T2HK than in the T2HKK. Adding a small Oki detector to the T2HK does not make any difference due to its limited statistics.
Our most prominent finding is related to (b), which is that adding a small 10 kton detector at Oki to the T2HK leads to revived sensitivity to the mass hierarchy.
This paper is organized as follows: In Section 2, we present the formalism for incorporating the effects of heavy isospin-singlet neutrino mixing with a non-unitary mixing matrix. In Section 3, we estimate the largest possible significance of heavy neutrino mixing in the T2HK, the T2HKK and the extension of the T2HK with a small detector at OKi, and further study the mass hierarchy and δ CP measurement in the presence of heavy neutrino mixing. Here, we describe our benchmark model, procedures for simulating the long baseline experiments, and analysis with a χ 2 fitting. Section 4 summarizes the paper.
2 Formalism for incorporating heavy isospin-singlet neutrino mixing where U αi is a N × N unitary matrix. In this paper, the N − 3 isospin-singlet neutrinos are assumed to be heavier than about 10 GeV. In long baseline neutrino oscillation experiments, only the light mass eigenstates ν 1 , ν 2 , ν 3 are produced and propagate. Therefore, the probability for a neutrino of active flavor α with energy E to transition to active flavor β after traveling a distance L inside matter is where U 3×3 denotes the α = e, µ, τ and i = 1, 2, 3 part of the mixing matrix U in Eq. (7), is because a phase shift due to time evolution, Edt, is equivalent to − E 2 − m 2 i dx for the mass eigenstate of mass m i . V nc (x) is the potential of the neutral current interaction inside matter divided by the neutrino velocity, and V cc (x) is that of the charged current interaction, which are evaluated to be where n n and n e respectively denote the neutron and electron number density, Y n and Y e are respectively the number of neutrons and electrons per one nucleon, which are typically 1/2, and ρ(x) is the matter density. [...] 3×3 in Eq. (10) signifies that one extracts the i, j = 1, 2, 3 part of the N × N matrix. It is convenient to parametrize U in terms of mixing angles, θ ij , and CP -violating phases, φ ij 1 , (i < j, i, j = 1, 2, ..., N) as [19,16] (for the case with N = 6, see also Ref. [20]) (Ω n,m ) kk = 1 for k = n, m, (Ω n,m ) nn = (Ω n,m ) mm = cos θ nm , (Ω n,m ) nm = sin θ nm e −iφnm , (Ω n,m ) mn = − sin θ nm e iφnm , (Ω n,m ) kl = 0 otherwise. (12) Note that the PMNS mixing matrix U P M N S corresponds to (Ω 2,3 Ω 1,3 Ω 1,2 ) with phases φ 12 , φ 23 set to 0 by the phase redefinition of |ν 1 , |ν 2 . We can thus express U in terms of the PMNS mixing matrix as matrices, respectively, α 11 , α 22 and α 33 are real numbers and α 21 , α 31 and α 32 are complex ones.
Since the mixing of isospin-singlet and active neutrinos is tiny, we make the approximation that we ignore terms in the second order of the mixing angles of active and isospin-singlet neutrinos. This gives (U † ) 3×3 = (U 3×3 ) † , and also We thus arrive at the following formula for the transition probabilities.
where we have made a further approximation with E 2 ≫ m 2 1 , m 2 2 , m 2 3 and discarded terms proportional to the unit matrix in H(x). The transition probabilities for antineutrinos, P (ν α → ν β ), are obtained by changing 3 Significance of heavy neutrino mixing and its effects on the mass hierarchy and δ CP measurement We estimate the largest possible significance of heavy neutrino mixing and its effects on the mass hierarchy and δ CP measurement, through the following simulation: First we introduce a benchmark model where α 11 , α 22 , |α 21 | in Eq. (2) are set at the current experimental bounds, in order to evaluate maximal impact of heavy neutrino mixing. Neutrino detection in the T2HK, the T2HKK, and the extension of the T2HK with a small detector at Oki is simulated based on the above benchmark, with Particle Data Group values employed for the standard parameters θ 12 , θ 23 , θ 13 , ∆m 2 21 , |∆m 2 31 |. Since no conclusive data on the mass hierarchy and δ CP exist, and the phase of α 21 is undetermined in the benchmark, we repeat the simulation for both mass hierarchies and for various values of δ CP and Arg(α 21 ). To assess merits of the extension of the T2HK with an Oki detector, we additionally consider, for comparative study, an experiment where the detector at Oki is instead placed at Kamioka. Finally, we perform a χ 2 fit of the simulation results under the assumption that no heavy neutrino mixing is present, namley, (13). The minimum of χ 2 is the square of the significance of heavy neutrino mixing. The mass hierarchy and the value of δ CP with which χ 2 is minimized correspond to their measured values obtained without assuming heavy neutrino mixing.

Benchmark model
Our analysis is based on the following benchmark model: We employ the parametrization of Eq. (13). Since ν µ → ν e process is most affected by α 11 , α 22 and α 21 , we further set At present, the most severe constraint on α 11 , α 22 and α 21 derives from experimental tests of lepton flavor universality and a mathematical inequality among α 11 , α 22 and α 21 . These are enumerated below: • One test of lepton flavor universality has been done in π ± decays. When the non-unitary matrix Eq. (15) enters into the neutrino mixing Eq. (7), the ratio of the branching fractions becomes where ǫ SM indicates a Standard Model correction taking into account final state radiation. The experimental value and the Standard Model prediction are respectively given by [21,22,23,27,28] Γ(π ± → e ± ν) from which we obtain the following bound: • Another test of lepton flavor universality comes from unitarity test of the Cabibbo-Kobayashi-Maskawa (CKM) matrix. In the presence of the non-unitary matrix Eq. (15), the Fermi constant measured in the muon decay becomes G µ = G F α 11 α 2 22 + |α 21 | 2 , whereas the Fermi constant measured in decays of hadrons involving eν becomes G β = G F α 11 . Deviation of G β /G µ from 1 mimics non-unitarity in the first row of the CKM matrix if |V ud | is measured in β-decays of nuclei and |V us | is measured in K → πeν decays, since G µ is regarded as the true Fermi constant in these measurements while it is G β that is involved in the processes, and |V ub | 2 is known to be smaller than experimental error of |V us | 2 . We therefore have Measurements of superallowed β-decays of nuclei yield the following experimental value [24]: Combining measurements of the branching ratios of only three processes K ± → π 0 eν, K L → π ± eν and K S → π ± eν [24], we obtain the following experimental value for |V us | times the kaon form factor, f + (0): The average of lattice calculations of the kaon form factor f + (0) in Refs. [25,26] is given by From Eqs. (20), (21), (22), we obtain the following bound: • The mathematical inequality among α 11 , α 22 and α 21 stems from the fact that the nonunitary mixing matrix N N U is derived from the product of mixing matrices as in Eq. (12). The inequality reads [27,28,29] We combine the 3σ experimental bounds in Eqs. (18), (23) and the mathematical inequality Eq. (24), and obtain the following 3σ bound 2 : To estimate the maximum significance of heavy neutrino mixing, we set the values of α 11 , α 22 and |α 21 | at the above 3σ bounds and therefore employ a benchmark model with the following non-unitary mixing matrix N N U : Since no experimental bound is reported on the phase of α 21 , we leave it as a free parameter.
We comment in passing on constraints from past neutrino oscillation experiments. Among those experiments, the NOMAD experiment [30], a short baseline experiment searching for ν µ → ν e andν µ →ν e transitions, gives the most severe bound, which is translated into the following: Our benchmark Eq. (26) satisfies the above bound. Finally, we mention constraints from lepton flavor violating processes. These constraints can always be evaded by tuning heavy neutrino mass m 4 , m 5 , m 6 , ..., m N around the W boson mass, because the amplitude of a lepton flavor violating process involving charged leptons ℓ α and ℓ β , A αβ , depends on N N U through the following combination: where m i , m j are neutrino mass eigenvalues, M W is the W boson mass, and F (x) is a function depending on the process, and in the second line, we make the approximation with m 1 , m 2 , m 3 ≪ M W as m 1 , m 2 , m 3 correspond to the tiny neutrino mass. It is evident that one can negate the term (N N U N † N U ) βα F (0) by taking appropriate O(1) values for m j /M W (j = 4, 5, 6, ..., ). We therefore neglect any constraints from lepton flavor violating processes.
In Appendix C, we consider another benchmark that satisfies the experimental bounds of Eqs. (18), (23) at 2σ level as well as the mathematical inequality Eq. (24).

Simulation
We simulate the neutrino propagation and detection based on the benchmark model Eq. (26).
The oscillation probability is calculated numerically with the formula Eq. (14) without making any approximation. For the standard oscillation parameters, we employ the central values of θ 12 , θ 13 , θ 23 , ∆m 2 21 , |∆m 2 32 | found in the Particle Data Group [24]. As the CP -violating phase δ CP and the sign of ∆m 2 32 have not been measured conclusively, and the phase of α 21 , Arg(α 21 ), is undetermined in the benchmark, we repeat the simulation for various values of δ CP and Arg(α 21 ) and for both cases with ∆m 2 32 > 0 and ∆m 2 32 < 0. The values of parameters used in the simulation are shown in Table 1.
We assume the baseline length L, fiducial volume, matter density along the baseline and neutrino beam off-axis angle [9,18] for the water Cerenkov detectors used in the T2HK, the T2HKK, the extension of the T2HK with a detector at Oki (T2HK+Oki), and the extension with the same detector at Kamioka (T2HK+Kami) as Table 2. The matter density is approximated to be uniform [31,18]. We assume that J-PARC operates with 1.3 MW beam power for 10 years in the neutrino mode and for another 10 years in the antineutrino mode, delivering 27 × 10 21 proton-on-target (POT) flux of neutrino-focusing beam and 27 × 10 21 POT flux of antineutrino-focusing beam, as Table 3. The energy distributions of ν µ andν µ in neutrino-focusing and antineutrino-focusing beams at each site are calculated based on Ref.
[32]. ν e andν e components in the beam are ignored in this analysis. In Appendix A, we present the energy distributions.
The cross sections for charged current quasi-elastic scattering between a neutrino and a proton, ν ℓ n → ℓ − p, and that between an antineutrino and a neutron,ν ℓ p → ℓ + n, (p and n denote proton and neutron, respectively, and ℓ denotes e or µ) are quoted from Ref. [33]. In Appendix B, we present the cross sections.
We assume the ideal measurement of ν e andν e and do not take into account acceptance and efficiency factors and finite energy resolution. Signals coming from processes other than charged current quasi-elastic scattering, as well as backgrounds due to neutral current quasielastic scatterings and other processes are not considered in the analysis. The total numbers of neutrino events at the detectors at Kamioka, Oki and in Korea in the T2HK, T2HKK, T2HK+Oki, and T2HK+Kami with the detector size of Table 2 and the   beam flux of Table 3 are tabulated in Table 4. We calculate the numbers of ν e andν e events in 0.05 GeV bins of the neutrino energy E in the range 0.4 GeV≤ E ≤3 GeV for various values of δ CP and Arg(α 21 ) and for both cases with ∆m 2 31 > 0 and ∆m 2 31 < 0. We remind that in real experiments, the neutrino flux is determined by the measurement of ν µ andν µ at the near detector, for which the standard oscillation is negligible but the non-unitary mixing reduces the ν µ andν µ flux by the following amount: In the analysis, we consider a virtual neutrino flux that is larger than the one found in Ref.
P denotes the transition probability Eq. (14), with E dependence made explicit. N n and N p are respectively the number of neutrons and protons in a water Cerenkov detector. σ denotes the cross sections for ν ℓ n → ℓ − p andν ℓ p → ℓ + n processes.

χ 2 analysis
We fit the numbers of events in the bins with neutrino-focusing and antineutrino-focusing beams, detected at Kamioka/Korea/Oki, N e,i,site andÑ e,i,site (i = 1, 2, ..., 52) (site=Kamioka, Korea, Oki), under the assumption that no heavy neutrino mixing is present, that is, N N U = I 3 . This is performed by minimizing χ 2 (Π), with respect to the standard oscillation parameters Π = (θ 12 , θ 13 , θ 23 , ∆m 2 21 , ∆m 2 32 , δ CP ). Here, χ 2 statistical represents statistical uncertainty, N unitary e,i,site (Π) is the number of ν e andν e events in a bin with a neutrino-focusing beam when N N U = I 3 , calculated as a function of the standard oscillation parameters Π, andÑ unitary e,i,site (Π) is the corresponding number with an antineutrinofocusing beam.
The following simplification is made for χ 2 parametrical (Π): We note that θ 12 , θ 13 , θ 23 , ∆m 2 21 , |∆m 2 32 | can be accurately measured in disappearance experiments withν e →ν e , ν µ → ν µ andν µ →ν µ , which are not influenced by heavy neutrino mixing as parametrized by N N U , because N N U equally affects the neutrino flux of the same flavor at the near and far detectors. (cos δ CP and the sign of ∆m 2 32 are also possibly measured in disappearance experiments, but the sensitivity is limited.) We therefore assume that the true values of θ 12 , θ 13 , θ 23 , ∆m 2 21 , |∆m 2 32 | are known prior to the experiments considered in this paper, and accordingly, we fix θ 12 , θ 13 , θ 23 , ∆m 2 21 , |∆m 2 32 | in χ 2 (Π) at the values used in the simulation and set χ 2 parametrical (Π) = 0. Finally, we have χ 2 systematic = 0 in the current analysis, because no acceptance, efficiency, energy resolution and background events are considered.
To summarize, χ 2 (Π) is approximated as χ 2 (Π) = χ 2 statistical ( δ CP , sgn(∆m 2 32 ) ) with θ 12 , θ 13 , θ 23 , ∆m 2 21 , |∆m 2 32 | given in Table 1. (32) We numerically evaluate the minimum of χ 2 (Π) Eq. (32), which corresponds to the square of the significance for heavy neutrino mixing. In Figure 1, we show contour plots of the minimum of χ 2 (Π), on the plane spanned by the CP -violating phase in the PMNS matrix δ CP and the phase in the non-unitary mixing matrix Arg(α 21 ), when the true mass hierarchy is normal. The four subplots correspond to the T2HK, the T2HKK, the extension of the T2HK with a detector at OKi, and that with the same detector at Kamioka. We study the impact of heavy neutrino mixing on the mass hierarchy measurement per-formed by fitting data without incorporating heavy neutrino mixing. For this purpose, we evaluate the difference between the minimum (with respect to δ CP only) of χ 2 (Π) for the wrong hierarchy and that for the true mass hierarchy, which corresponds to the square of the significance of rejecting the wrong mass hierarchy. We show in Figure 2 contour plots of the difference, on the same plane when the true mass hierarchy is normal. The black-filled regions in the left subplots are where the wrong mass hierarchy is favored over the true mass hierarchy, i.e., min χ 2 (wrong H) − min χ 2 (true H) < 0. We study the impact of heavy neutrino mixing on the δ CP measurement performed by fitting data without incorporating heavy neutrino mixing. To assess the impact, we compute the difference between the true value of δ CP , and the value that minimizes χ 2 (Π) of Eq. (32) with the true mass hierarchy. We show in Figure 3 contour plots of the difference, ∆δ CP ≡ (δ CP that minimizes χ 2 (δ CP , true sgn(∆m 2 32 )) − (true δ CP ) on the same plane when the true mass hierarchy is normal.  The following observations are made: Figures 1,4, we see that the significance of heavy neutrino mixing can be above 3σ for Arg(α 21 ) ∼ δ CP , but decreases considerably for Arg(α 21 ) ≃ δ CP ± π/2 in all the experiments, confirming that the sensitivity deteriorates due to suppression of interference between the standard oscillation amplitude and non-unitary mixing α 21 at the first oscillation peak, as can be read from Eqs. (5), (6). The regions of depleted sensitivity is slightly off the lines Arg(α 21 ) = δ CP ± π/2, because of the influence of the second, subleading interference term in Eqs. (5), (6). We also find that the significance does not ameliorate in the T2HKK and in the extension of the T2HK with a 10 kton detector at Oki. Rather, the T2HKK shows inferior performance due to the loss of statistics, which is seen in Table 4.   Figure 7 indicates that for δ CP ∼ 0, π, the T2HK has limited sensitivity to the mass hierarchy in the unitary case, and the presence of non-unitary mixing may aggravate it to the extent that the wrong mass hierarchy is favored. Such impact of non-unitary mixing is because α 21 cancels matter effects in the neutrino transition probability Eq. (3) when (δ CP , Arg(α 21 )) ∼ (0, 0), (π, π) holds in the normal hierarchy case, and when (δ CP , Arg(α 21 )) ∼ (π, 0) holds in the inverted hierarchy case. Although α 21 and matter effects add positively in the antineutrino transition probability Eq. (4) in the same cases, it does not change the situation because neutrino flux is larger than antineutrino flux in our setup (see Table 3 and Figure 9).
The T2HKK maintains strong sensitivity to the mass hierarchy in the presence of heavy neutrino mixing. Also, adding a 10 kton detector at Oki to the T2HK drastically ameliorates sensitivity to the mass hierarchy, in spite of the small size of the Oki detector. The advantage of the T2HKK and T2HK+Oki is because measurements with different matter effects resolve degeneracy between non-unitary mixing α 21 and matter effects.
We stress that even a small 10 kton detector at Oki negates the impact of heavy neutrino mixing on the mass hierarchy measurement and allows correct measurement.
(C) In Figures 3,6, we find that the value of δ CP that minimizes χ 2 under the assumption of standard unitary oscillation can deviate from the true δ CP by more than 0.15 radian in the T2HK, and by more than 0.1 radian in the T2HKK. The result should be compared to the resolution of δ CP . To this end, we present in Figure 8 the difference between δ CP that gives χ 2 = 1 and the true δ CP when non-unitarity is absent (N N U = I), (∆δ CP at 1σ) ≡ (δ CP that gives which quantifies 1σ resolution of the δ CP measurement in the standard oscillation case. Eq. (36) is calculated through a simulation performed with the same procedure.  Figure 8: The resolution of the δ CP measurement when non-unitarity is absent, quantified as the difference between δ CP that gives χ 2 = 1 and the true δ CP Eq. (36). The results are presented as functions of the true δ CP . The left and right plots correspond to the cases when the true mass hierarchy is normal and inverted, respectively. In each plot, the solid line indicates the sensitivity in the T2HK, and the dashed line the sensitivity in the T2HKK.
The 1σ resolution for δ CP is better than 0.1 radian in the T2HK, and better than 0.08 radian in the T2HKK, which is of the same order as the deviation of δ CP due to non-unitarity. This indicates that non-unitarity can have a non-negligible effect on the measurement of δ CP , and hence, if a hint of non-unitarity has been discovered, it must be incorporated in the data fitting in order to measure δ CP correctly.
The deviation due to non-unitarity is enhanced when (true δ CP ) ∼ Arg(α 21 ) holds. This is because interference between the part of the standard oscillation amplitude proportional to e iδ CP and the amplitude involving non-unitary mixing α 21 is maximized, and so is the influence of α 21 on the value of δ CP that minimizes χ 2 .
The deviation is smaller in the T2HKK than in the T2HK. This is due to the fact that the energy dependence of the part of the standard oscillation amplitude proportional to e iδ CP , which helps distinguish the standard and non-unitary amplitudes, is more accurately measured with a longer baseline. Hence, the δ CP measurement is less affected by non-unitarity in the T2HKK.

Summary
We have studied the discovery potential for the mixing of heavy isospin-singlet neutrinos in the Tokai-to-Hyper-Kamiokande (T2HK), the Tokai-to-Hyper-Kamiokande-to-Korea (T2HKK), and a plan of adding a small detector at Oki Islands to the T2HK, and further examined the feasibility of measuring the mass hierarchy and the standard CP -violating phase δ CP in the presence of heavy neutrino mixing by fitting data without assuming heavy neutrino mixing. The mixing of heavy neutrinos is parametrized with a non-unitary mixing matrix for active flavors. A benchmark that maximizes the non-unitary mixing and is consistent with the current experimental and theoretical bounds is employed to estimate the largest possible significance of heavy neutrino mixing and its impact on the mass hierarchy and δ CP measurement.
Through a simulation, we have revealed that the significance of heavy neutrino mixing can be above 3σ in all the experiments when the standard phase δ CP and a new CP -violating phase originating from heavy neutrino mixing, Arg(α 21 ), satisfy Arg(α 21 ) ≃ δ CP . The significance decreases considerably for Arg(α 21 ) ≃ δ CP ± π/2, due to suppression of interference between the standard oscillation amplitude and a non-unitary mixing term.
In the T2HK, the mass hierarchy measurement is so affected by heavy neutrino mixing that the wrong mass hierarchy is favored for (δ CP , Arg(α 21 )) ∼ (0, 0), (π, π) in the normal hierarchy case, and for (δ CP , Arg(α 21 )) ∼ (π, 0) in the inverted hierarchy case. This is because a nonunitary mixing term possibly cancels matter effects. In contrast, the T2HKK and the extension of the T2HK with a 10 kton detector at Oki show strong sensitivity to the mass hierarchy even with heavy neutrino mixing, because measurements with different baseline lengths at Kamioka and in Korea or at Oki resolve degeneracy between a non-unitary mixing term and matter effects. We thus conclude that (i) the mass hierarchy measurement in the T2HKK is highly stable against any effects of heavy neutrino mixing, and that (ii) although the mass hierarchy measurement in the T2HK is easily affected, the addition of a small detector at Oki negates effects of heavy neutrino mixing.
In the presence of non-unitarity, the value of δ CP measured without assuming non-unitarity deviates from the true δ CP . The deviation is enhanced when δ CP and Arg(α 21 ) are similar, because in such a case, the part of the standard oscillation amplitude proportional to e iδ CP interferes maximally with a non-unitary mixing term at oscillation peaks. The deviation can be of the same order as 1σ resolution of the δ CP measurement, which suggests that if a hint of non-unitary mixing has been discovered, one must employ an appropriate ansatz including the non-unitary mixing terms to fit neutrino oscillation data and measure δ CP correctly.

Appendix A
In Figure 9, we show the flux of ν µ andν µ in neutrino-focusing and antineutrino-focusing beams from J-PARC, detected at a water Cerenkov detector at Kamioka, Oki and in Korea if the neutrino oscillation were absent. The baseline length and beam off-axis angle assumed are given in Table 2. In each plot, the blue lines correspond to a neutrino-focusing beam and the red lines correspond to an antineutrino-focusing beam. The solid lines denote ν µ flux and the dashed lines denoteν µ flux.

Appendix B
In Figure 10, we show the cross sections for charged current quasi-elastic scatterings ν ℓ n → ℓ − p andν ℓ p → ℓ + n (ℓ = e, µ). The solid line corresponds to the cross section for ν ℓ , and the dashed line to that forν ℓ .

Appendix C
We consider an alternative benchmark with a smaller non-unitary mixing that satisfies the experimental bounds Eqs.
Based on the benchmark with Eq. (37), we repeat the same simulation and calculation as subsection 3.2 and 3.3, and obtain the results in Figures 11, 12, 13, 14, 15, 16.
The plots for minχ 2 , Figures 11, 14, show that the significance of non-unitary mixing is reduced below 3σ in most parameter regions in the T2HK and in the entire region in the T2HKK. The dependence on (δ CP , Arg(α 21 )) is similar to the plots for the original benchmark, Figures 1, 4.
The plot for {min χ 2 (wrong H) − min χ 2 (true H)} in the normal hierarchy case, Figure 12, tells us that in the T2HK, the wrong mass hierarchy can still be favored if the true hierarchy is normal. In contrast, Figure 15 gives that this is no longer the case if the true hierarchy is inverted.
The plots for ∆ CP , Figures 13, 16, show that the value of δ CP that minimizes χ 2 (Π) of Eq. (32) still possibly deviates from the true value by more than 0.1 radian in the T2HK, and more than 0.05 radian in the T2HKK, even with the alternative benchmark with the smaller non-unitary mixing. This augments our argument that the non-unitary mixing terms must be incorporated in the δ CP measurement when a hint of non-unitarity has been discovered.  Figure 11: The minimum of χ 2 (Π) Eq. (32), min χ 2 , on the plane of δ CP and Arg(α 21 ). The benchmark with Eq. (37) is assumed, and the true mass hierarchy is normal. The upper-left, upper-right and lower-right subplots correspond to the T2HK, the T2HKK, and the plan of the T2HK plus a 10 kton water Cerenkov detector at Oki, respectively. For comparative study, we show in the lower-left a subplot for a plan of the T2HK plus a 10 kton water Cerenkov detector at Kamioka. min χ 2 = 1, 4, 9 on the red solid, blue dashed, and purple dot-dashed contours, respectively.   Figure 13: The difference between the true value of δ CP , and the value that minimizes χ 2 (Π) of Eq. (32) with the true mass hierarchy. The benchmark with Eq. (37) is assumed, and the true mass hierarchy is normal. The upper-left, upper-right and lower-right subplots correspond to the T2HK, the T2HKK, and the plan of the T2HK plus a 10 kton water Cerenkov detector at Oki, respectively. For comparative study, we show in the lower-left a subplot for a plan of the T2HK plus a 10 kton water Cerenkov detector at Kamioka. |∆δ CP | = 0, 0.05rad, 0.1rad, and 0.15rad on the black solid, red sold, blue dashed, and purple dot-dashed contours, respectively.  37) is assumed, and the true mass hierarchy is inverted. The upper-left, upper-right and lower-right subplots correspond to the T2HK, the T2HKK, and the plan of the T2HK plus a 10 kton water Cerenkov detector at Oki, respectively. For comparative study, we show in the lower-left a subplot for a plan of the T2HK plus a 10 kton water Cerenkov detector at Kamioka. min χ 2 = 1, 4, 9 on the red solid, blue dashed, and purple dot-dashed contours, respectively.   Figure 16: The difference between the true value of δ CP , and the value that minimizes χ 2 (Π) of Eq. (32)∆δ CP with the true mass hierarchy. The benchmark with Eq. (37) is assumed, and the true mass hierarchy is inverted. The upper-left, upper-right and lower-right subplots correspond to the T2HK, the T2HKK, and the plan of the T2HK plus a 10 kton water Cerenkov detector at Oki, respectively. For comparative study, we show in the lower-left a subplot for a plan of the T2HK plus a 10 kton water Cerenkov detector at Kamioka. |∆δ CP | = 0, 0.05rad, 0.1rad, and 0.15rad on the black solid, red sold, blue dashed, and purple dot-dashed contours, respectively.