Testable Baryogenesis in Seesaw Models

We revisit the production of baryon asymmetries in the minimal type I seesaw model with heavy Majorana singlets in the GeV range. In particular we include"washout"effects from scattering processes with gauge bosons, Higgs decays and inverse decays, besides the dominant top scatterings. We show that in the minimal model with two singlets, and for an inverted light neutrino ordering, future measurements from SHiP and neutrinoless double beta decay could in principle provide sufficient information to predict the matter-antimatter asymmetry in the universe. We also show that SHiP measurements could provide very valuable information on the PMNS CP phases.


Introduction
It is well known that minimal extensions of the Standard Model that accommodate massive neutrinos, such as the type I seesaw models, could also explain the observed matterantimatter asymmetry in the universe [1]. The two new ingredients that make this possible are the existence of new particles that are not in thermal equilibrium sometime before the electroweak phase transition (T EW 140GeV) and the presence of new CP-violating interactions of these particles. Two basic scenarios have been shown to work. In the first one, the heavy Majorana singlets decay out of equilibrium generating a lepton asymmetry that sphaleron processes recycle into a baryonic one. These neutrinos have masses well above the electroweak scale, typically M 10 8 − 10 9 GeV [2,3] and M 10 6 GeV when B − L is almost conserved [4], while in resonant leptogenesis masses in the TeV scale are possible [5]. For a comprehensive review and references of this very well studied scenario see ref. [6]. In the second scenario, the heavy Majorana singlets have masses below the electroweak scale, and therefore their Yukawa couplings are small enough that one or more of these states might not reach thermal equilibrium by the time the electroweak phase transition takes place. A lepton (and baryon) asymmetry can be generated when the states are being produced, i.e. at freeze-in. The sterile states get populated via Yukawa interactions, but the coherence between collisions is essential to produce a CP asymmetry, via the interference of CP-odd and CP-even phases. This is why this mechanism is often referred-to as baryogenesis from neutrino oscillations. It was first proposed by Akhmedov, Rubakov and Smirnov (ARS) in their pioneering work [7] and pursued, with important refinements in references [8,9]. A list of recent references is [10][11][12][13][14][15].
In [15], we studied this second scenario and explored the available phase space for successful leptogenesis in the minimal models with two or three extra singlets, N = 2, 3. In particular we considered an accurate analytical approximation, where we could identify the relevant CP invariants, and that helped us explore the full parameter space. The case with N = 2 is effectively equivalent to the popular νMSM [8], which is a N = 2 + 1 model, where the lighter neutrino plays the role of dark matter and decouples from the problem in the generation of the baryon asymmetry. The original ARS scenario on the other hand required the interplay of all the three species. The analytical approximation of [15] allowed to clarify these different scenarios.
The purpose of this paper is twofold. First, we will refine our previous study in various aspects. In [15] (like in most previous works) the collision terms only included the dominant top quark scatterings. As has been known for sometime [16,17], scatterings off gauge bosons, as well as the resumed decays and inverse decay processes, are also very important. These rates have been computed in [17,18] in the limit of vanishing leptonic chemical potentials. In the generation of lepton asymmetries, it is very important however to include the effect of the latter, since these will tend to washout the asymmetry. In section 2 we derive new kinetic equations including all the scattering processes considered in [17,18], that we have re-evaluated in the presence of small leptonic chemical potentials. Furthermore Fermi-Bose statistics is consistently used through-out.
The second important improvement concerns our scans of parameter space. In our previous study we speeded-up the scan using the analytical approximation. This forced us to avoid some regions in parameter space, to ensure that the approximation was good enough. We have now optimized significantly the numerical solution of the kinetic equations, in particular addressing the stiffness problem. The analytical approximation is no longer needed, and therefore the ad hoc constraints on parameter space are avoided. We use a Bayesian approach to extract posterior probabilities on the relevant observables of the model, from a prediction of the measured baryon asymmetry, using the Multinest package. In section 3, we present the results of these scans of parameter space.
In the second part of the paper, we address the question: to what extent it would be possible to predict quantitatively the baryon asymmetry, within the minimal model N = 2, if the heavy neutrino states would be discovered in future experiments, such as SHiP. In section 4, we derive approximate analytical formulae valid in the range within SHiP reach, which demonstrate the complementarity of the different measurements: mixings and masses of the extra states from direct searches, neutrinoless double beta decay and the CP phase measurable in neutrino oscillations. The numerical study confirms these expectations and allows us to answer the question in the affirmative if nature is kind enough to provide us with positive signals at SHiP and an inverted neutrino ordering. Furthermore we show how such SHiP measurements could constrain the CP-violating phases of the PMNS matrix.

Kinetic equations
The Lagrangian of the model is given by: where Y is a 3×3 complex matrix and M a symmetric matrix. One convenient parametrization is in terms of the eigenvalues of the Y and M matrices, together with two unitary matrices, V and W . In the basis where the Majorana mass is diagonal, M = Diag(M 1 , M 2 , M 3 ), the neutrino Yukawa matrix is given by: Without loss of generality, using rephasing invariance, we can reduce the unitary matrices to the form: Obviously not all the parameters are free, since this model must reproduce the light neutrino masses, which approximately implies the seesaw relation: where v = 246 GeV is the vev of the Higgs. A very convenient parametrization that takes this constraint into account is the Casas-Ibarra one [19], where the Yukawa matrix can be written in terms of the light neutrino masses and mixings as where m light is a diagonal matrix of the light neutrino masses, U PMNS (θ 12 , θ 13 , θ 23 , δ, φ 1 , φ 2 ) is the PMNS matrix that describes the light neutrino mixing, M is the diagonal matrix of the heavy neutrino masses, and R is a complex orthogonal matrix, that depends generically on one (three) complex angle(s) z ij for N = 2 (N = 3). The kinetic equations that describe the production of sterile neutrinos in the early Universe have been studied in many previous works, see for example [11][12][13]. In this work we have rederived these equations with the following refinements with respect to our previous work [15]: • Fermi-Dirac or Bose-Einstein statistics is kept throughout • Collision terms include 2 ↔ 2 scatterings at tree level with top quarks and gauge bosons, as well as 1 ↔ 2 scatterings including the resummation of scatterings mediated by soft gauge bosons as obtained in refs. [16][17][18] • Leptonic chemical potentials are kept in all collision terms to linear order • Include spectator processes As usual we assume that all the spectator particles are in kinetic equilibrium. On the other hand, we neglect the effects of the top quark and Higgs chemical potentials. These effects are expected to be smaller than the effect of thermal masses in 2 ↔ 2 processes that we are neglecting. Note that, in contrast with the effects of the lepton chemical potential, the former do not bring in any new flavour structure. The starting point to derive the equations is the Raffelt-Sigl formalism [20], where the sterile neutrino density satisfies the equation: and Γ a N (k) and Γ p N (k) are the annihilation and production rates of the sterile neutrinos. The result can be written as where ρ F (y) = (exp y + 1) −1 is the Fermi-Dirac distribution and µ α is the leptonic chemical potential normalised by the temperature. γ N contain the contributions from all 2 → 2 processes that produce an N : and 1 ↔ 2 processes: φ →lN including resummed soft-gauge interactions. All these contributions have been computed for vanishing leptonic chemical potential in [10,17,18]. We have followed their methods including the effects of a lepton chemical potential to linear order. Defining and γ (1) with ρ F (y) ≡ dρ F (y) dy , the functions γ (i) N get contributions from quark (Q), gauge scattering (V) and the 1 → 2 resummed processes (LPM): The functions γ LP M has non-trivial temperature dependence due to the runnings of the coupling constants 2 . In Fig. 1 the three functions are plotted, where the two lines labeled LPM curve correspond to two temperatures 10 4 GeV and 10 10 GeV, while where ρ B (y) = (exp y − 1) −1 is the Bose-Einstein distribution. Inserting these functions in the kinetic equation we get: where µ ≡ Diag(µ α ). The equation for the antineutrino (opposite helicity state) density is the same but changing µ → −µ and Y → Y * . It is often useful to consider instead the evolution of the CP conserving and violating combinations: The equations for these combinations are: Finally we need the equations that describe the evolution of the leptonic chemical potentials. This is obtained from the equation that describes the evolution of the conserved charges in the absence of neutrino Yukawas, that is the B 3 − L α numbers. These numbers can only be changed by the same out of equilibrium processes that produce the sterile neutrinos: T − µ α , it is possible to relate the integrated rates of these equations to those of the sterile neutrinos and their densities: 18) or in terms of ρ ± : where I α is the projector on flavour α, and we have neglected terms of O(µρ − ).
The relation between the leptonic chemical potentials and the approximately conserved charges, B/3 − L α , is given for T ≤ 10 6 GeV by [21] where we have defined µ B/3−L β by the relation: Introducing finally the expansion of the universe and changing variables to the scale factor x = a and y = ka, the time derivative of the distribution functions change to: M Planck is the Hubble expansion parameter. Assuming a radiation dominated universe with constant number of relativistic degrees of freedom g * (T 0 ) 106.75 for T 0 ≥ T EW , then xT = constant that we can fix to one.

Momentum averaging
In principle the equations should be solved for all momenta of the sterile neutrinos, but it is a good approximation [10] to assume ρ ± (x, y) = r ± (x)ρ F (y), with r ± (x) independent of momentum. This allows to integrate explicitly over y and we get the momentum-averaged equations: The momentum averaged rates are:    Table 1. Coefficients in the momentum averaged rates. The LPM ones have been evaluated at T 1 = 10 4 GeV and T 2 = 10 10 GeV. with , (2.27) and the coefficients are given in the table 1.
For the couplings g, g we evaluate them at the scale πT : while the top Yukawa running is obtained numerically from the one loop renormalization group equations.
In Fig. 2 we compare the time evolution of the asymmetry that we obtain with the new equations and the old equations of [15] that only included top scattering processes and Maxwell-Boltzmann statistics. As expected the larger scattering rates induce a larger asymmetry at short times, but also a stronger washout at late times.
In [15], we identified four independent CP rephasing invariants that can contribute to this asymmetry in the general case with N = 3 as: where V, W are the matrices parametrizing the neutrino Yukawa matrix, eq. (2.1). In the minimal scenario with N = 2 only the first two invariants can contribute. We considered a convenient analytical approximation, based on a perturbative expansion in the mixing angles of these matrices, that allowed us to solve the differential equations analytically, neglecting non-linear terms. It is straightforward to apply the same method to the new equations. As an example we give the result for the asymmetry when only the CP invariant I 1 survives (i.e. for φ i3 =φ i3 = 0). Defining with M * P ≡ M Planck 45 4π 3 g * (T 0 ) , and neglecting the running of the couplings, the result is and The integrals are  Assuming parameters for which only the simplest invariant, I is non-vanishing, comparison of the analytical result of eq.(2.34) (dashed), the full numerical one (red), the numerical one neglecting the variation of the running of the couplings (blue) and that neglecting also nonlinear terms (green). The parameters have been chosen as , and sufficiently small V, W mixings. Fig. 3 shows the analytical result compared with the numerical solution to the equation for sufficiently small mixing angles.

Baryon asymmetry
The observed baryon asymmetry is usually quoted in terms of the abundance, which is the number-density asymmetry of baryons normalised to the entropy density. After Planck this quantity is known to per cent precision [22]: The baryon abundance is related to that of B − L by [23,24] and where we take g * = 106.75 (which ignores the contribution to the entropy of the sterile states). Our estimate for the baryon asymmetry is therefore

Numerical Results
The numerical solution of the kinetic equations is challenging, because in most part of the parameter space, they are stiff since the oscillation time is much shorter than the collision one. In [15] we performed an exploration of the parameter space viable for leptogenesis for the models N = 2, 3 employing the perturbative analytical approximation. This required to constrain certain regions of parameter space where the perturbative solution could fail. We want to improve on this scan by going beyond the perturbative estimate of the asymmetry and using the full numerical solution of the equations. We have solved eqs. (2.25) using the publicly available code SQuIDS [25,26]. The code is designed to solve the evolution of a generic density matrix in the interaction picture. The interaction picture is useful because it removes the short time scale, i.e. the oscillation scale, from the numerical integration, but fast oscillatory coefficients then appear in the terms involving the off-diagonal elements of the density matrix. In order to optimize the code, at some large enough time, we switch to a fully decoherent evolution (when the exponents of all the oscillatory terms are larger than 10 5 ). The decoherent evolution is already included in the last version of the SQuIDs code [26]. Using this approximation the solution speeds up the computation by a factor more than a hundred and the result agrees with the full solution with a relative error smaller than O(1%).
Using these optimizations we get a computational time for the full numerical solution of order minutes, which allows us to do a Bayesian parameter estimation from the loglikelihood: (3.1) For this, we use a nested sampling algorithm implemented in the public package Multi-Nest [27][28][29] and the Markov Chain sample analysis tool GetDist [30] to get the posterior probabilities. The number of random starting points is 5000. The scan is performed using the Casas-Ibarra parameters of eq. (2.5). We fix the light neutrino masses and mixings to the present best fit points in the global analysis of neutrino oscillation data of ref. [31] for each of the neutrino orderings (normal, NH, and inverted, IH), and leave as free parameters: the complex angle(s) of the R matrix, the CP phases of the PMNS matrix, the lightest neutrino mass as well as the heavy Majorana masses. For N = 2 these are six independent parameters, while for N = 3, there are thirteen free parameters.
In this work we consider the simplest case of N = 2, which can be obtained from the N = 3 model in the limit where one of the sterile neutrinos is effectively decoupled, that we can assume without loss of generality to be N 3 . This can be achieved with the choice of parameters: for the IH(NH), where P is the 123 → 312 permutation matrix (only necessary for the NH).
This is then the model that has been considered in most previous works on the subject [8-10, 13, 32], where the number of constrained parameters is reduced to six: only one complex angle in R, z ≡ θ + iγ, two CP phases, δ and φ 1 in the PMNS matrix, and two Majorana neutrino masses, M 1 , M 2 . Figures 4 and 5 show, for IH and NH, the posterior probabilities of the spectrum of the two relevant states, M 1 , M 2 , the active-sterile mixings of the first heavy state |U α4 | 2 (those of the second state are almost identical), the neutrinoless double beta decay effective mass |m ββ | and the baryon asymmetry Y B . An important consideration are the priors. We have considered flat priors in all the Casas-Ibarra parameters except the masses where we assume a flat prior in log 10  Fig. 4 correspond to the two options. The significant differences between the two posteriors show the effect of allowing or not for fine-tuning in the degeneracy of the two heavy states. Even though the contours are typically larger if more fine-tuning is allowed, we find interesting solutions with a mild degeneracy, which tend to imply smaller M 1 , M 2 , larger values of the active-sterile mixing parameters and a sizeable non-standard contribution to neutrinoless double beta decay, which obviously imply much better chances of testability. Figs. 6 zoom in the most interesting results from this study: the mild level of fine-tuning of the blue contours (neither a strong degeneracy is required, nor a very large deviation from the naive seesaw scaling of the mixings), correlated with a relatively large mixing, and a sizeable amplitude for neutrinoless double beta decay. We will come back to this point in the next section. In table 2 we show the 68% probability ranges for the relevant parameters as extracted from the 1D posterior probabilities, for the two neutrino orderings and the two prior choices. The ranges for the mixings of the second heavy state, |U α5 | 2 , are basically the same.  Table 2. For the minimal model N = 2: 68% posterior probability ranges of log 10 (param) assuming flat prior in log 10 (M 2 (GeV)) (M) or log 10 (∆M 12 (GeV)) (∆M).
In Figures 7 we zoom in the probability plots for the heavy neutrino mixings versus mass and compare them with present [33], [34] and future constraints from DUNE [35], SHiP [36] and FCC [37]. Constraints from Big Bang Nucleosynthesis are very restrictive in the low mass range, particularly below the pion threshold [38,39]. In ref. [15] similar figures were shown from a scan of parameter space assuming only flat priors in log 10 M 1  The blue and red contours correspond respectively to the assumption of a flat prior in log 10 M 1 and log 10 M 2 and to a flat prior in log 10 M 1 and log 10 (∆M 12 ). The star is the test point used for the SHiP study of the next section. and log 10 |M 2 − M 1 |. We note that the regions we show here are the result of a full numerical treatment, that avoids any constraint in parameter space and successfully explain the baryon asymmetry within its small 1% uncertainty. The most important addition is however that of the blue contours that use flat priors in log 10 M 1 and log 10 M 2 , and therefore avoid too large fine-tuning. These solutions point to a region of parameter space within SHiP reach as the most probable one. It is interesting that the sensitivity of SHiP and DUNE to the e or µ channels will cover to a large extent the blue regions. When a larger degree of degeneracy in the masses is allowed (red regions), the right baryon asymmetry can be obtained also for larger masses, up to 10-100 GeV, but this high mass region will be harder to test experimentally (for recent work see also [40]).
In the N = 3 case, there are 13 unknown parameters and the exploration of parameter space is significantly more challenging. This case will be considered elsewhere.  Figure 6.
Posterior probabilities for the amplitude of neutrinoless double beta decay (left), electron mixing (middle) and α=e,µ,τ |U α4 | 2 M 1 (right) versus the mass degeneracy.  Comparison of the posterior probability contours at 68% and 90% on the planes mixings with e, µ, τ versus masses, with the present (shaded region) and future constraints from DUNE, FCC and SHiP for NH (up) y IH (down).

Predicting the baryon asymmetry in the minimal N = 2 model
A very relevant question is whether the baryon asymmetry could be predicted in this scenario if the heavy sterile neutrinos are within reach of future experiments, such as the SHiP experiment. We will analyse this question in the simplest case N = 2 where the number of unknown parameters is minimal. Obviously the situation for the N = 3 case will be much more difficult.
The SHiP experiment will be capable of detecting heavy neutrinos in the few GeV range provided their mixings are sufficiently large. In particular significantly larger than what the naive seesaw scaling |U ai | 2 ∼ m i /M i would suggest. In the Casas-Ibarra parametrization of eq. (2.5), this implies that the entries of the R matrix need to be significantly larger than one, and therefore the imaginary part of the complex angle needs to be sizeable.
In order to understand the dependence of Y B on the different parameters, it is useful to consider the perturbative results of [15]. The CP asymmetries responsible for the baryon number generation, ∆ CP , in the weak washout regime, can be generically written as For the N = 2 case, when φ i3 = 0 and y 3 = 0, this quantity can be written as [15] Y B = ∆ CP = y 1 y 2 (y 2 2 − y 2 1 ) (y 2 2 − y 2 1 )I where the invariants I are defined in eq. (2.32). This is indeed the dependence obtained from the solution of the kinetic equations in the perturbative approximation obtained in [15], in the weak washout limit. For our new kinetic equations, the perturbative result is shown in eq. (2.34) (where only the I (2) 1 contribution has been kept). We can read from eq. (2.34) in the limit of weak washout, Γ ± t 1, and using eq. (2.42): In contrast with eq. (2.34), this approximation fails at large times since it is valid only for Γ ± t 1. It typically overestimates the asymmetry, but should give qualitatively the right dependence on the parameters.
What we need however are the expressions in terms of the Casas-Ibarra parameters. The relations are typically very complicated, but we can identify a few small parameters and perturb in them: where γ, assumed positive, is the imaginary part of the complex angle of the R matrix that needs to be large to avoid the naive seesaw scaling of the active-sterile mixings 3 . Defining 3 Note that γ can also be negative, but there is an approximate symmetry γ → −γ, that would lead to very similar results by expanding in e − |γ| 2 in this case.
the result for the heavy-light mixing for IH can be written as The leading order however is not precise enough. For IH, it is clear from eq. (4.9) that depending on the value of φ 1 a significant suppression of the leading terms in the numerator or denominator can take place, therefore at this point the NLO corrections are relevant and these bring a new undetermined parameter, δ, as can be seen from eqs. (4.7), which is the CP phase that can be measured in neutrino oscillations! In this case, the measurement of the ratio of mixings at SHiP cannot resolve φ 1 and δ simultaneously and a degeneracy between these two phases remains. We will come back to this interesting observation in the following section. Clearly the δ phase could be determined in future neutrino oscillation experiments.
At leading order in the expansion, the CP asymmetry in this regime can be approximated by (sin 2θ cos 2θ 12 − cos φ 1 cos 2θ sin 2θ 12 ) × sin 2 2θ 23 + (4 + cos 4θ 23 ) sin φ 1 sin 2θ 12 + O( ) , where in the case of NH we have included NLO corrections because the LO cancel for maximal atmospheric mixing. We see that for both neutrino orderings, it depends strongly on the real part of the Casas-Ibarra angle θ, which does not affect any other of the observables above. In particular, independently of the value of φ 1 , there is always a value of θ that can make the asymmetry vanish. For instance for the IH result the value can be approximated by IH : tan 2θ cos φ 1 tan 2θ 12 , (4.12) therefore the uncertainty in θ forbids to set a lower bound on the asymmetry, although an upper bound can be derived. Therefore even if the sterile states would be discovered at SHiP and their mixings to electrons and muons accurately measured, the asymmetry can only be predicted up to this angle. Furthermore as we have seen in order to explain the baryon asymmetry in the N = 2 case, a significant level of degeneracy between the two states is needed. The dependence on this quantity enters in the function g(M 1 , M 2 ). Although we have not found a detailed study of the expected uncertainty in the invariant mass at SHiP, given the momentum resolution for muons and pions quoted in [41], we expect that this uncertainty could be better than 1%. If the degeneracy cannot be measured, a large uncertainty in the CP asymmetry will result also from this unknown.
Interestingly neutrinoless double beta decay also depends on both unknowns: θ and ∆M 12 . The effective neutrino mass in neutrinoless double beta is given approximately by [42,43] |m ββ | IH ∆m 2 atm c 2 13 c 2 12 + e 2iφ 1 s 2 12 1 + Note that the non-standard contribution from the heavy states is very sensitive to the mass degeneracy. Furthermore the interference between the light and heavy contributions is very sensitive to the angle θ, and it is precisely in the region around 1 GeV where the heavy and light contributions can be comparable, and can effectively interfere [43,46,47]. There is therefore the possibility that neutrinoless double beta decay could provide the missing information to predict the baryon asymmetry. On the left plot of Fig. 8 we show |m ββ | as a function of the angle θ for IH and some choice of parameters that are within the range of SHiP and assumed known. |m ββ | has been computed exactly using the nuclear matrix elements for 76 Ge from ref. [44].The bands are the standard 3ν contributions to |m ββ | for NH/IH. If the uncertainty in |m ββ | could be better than the spread within the standard IH region, a determination of θ could result from this measurement. On the right plot we show the dependence of Y B (computed exactly) on the same angle. In Fig. 9 we show the correlation between |m ββ | and Y B as we vary θ. Ideally a precise determination of |m ββ | could predict the baryon asymmetry up to a global sign. In practice, this would require a very good control of the nuclear matrix elements which is very challenging. As a proof of principle we have studied the posterior probabilities for a hypothetical measurement of SHiP of M 1 and M 2 and their respective couplings to electrons and muons for IH. The point chosen is indicated by a star in Fig. 4. We did not look for a very special nor optimal point, simply that it was within the range of SHiP and could explain the baryon asymmetry. The corresponding triangle plots are shown in Fig. 10 for two values of the assumed errors in this experiment: relative errors on masses and mixings of (1%, 10%) and the much more optimistic (0.1%, 1%). We furthermore considered a combination of SHiP, with the more optimistic errors, with a determination of δ in future neutrino oscillation experiments such as HyperK and DUNE. We have assumed σ δ 10 • as derived from the studies in references [48,49]. The |m ββ | versus Y B plot is zoomed in in Fig. 11.
Clearly the determination of δ results in a more clear correlation between |m ββ | and Y B . In fact the resulting width of the contour can be understood from the propagation of the error in δ on the determination of φ 1 . This is shown in Fig. 12, where we compare the posterior probability with the three curves obtained in the following way: 1) fixing the parameters to those of the test point and changing only θ (solid line), 2) the same as 1) except δ which is fixed to δ test − σ δ , and correspondingly φ 1 , γ tuned to keep the mixings unchanged. There are two solutions for φ 1 and the two dashed lines correspond to the two values (see Fig. 13), according to the analytical approximate formulae. The region encompassed by these lines is roughly similar to the 68% and 90% regions. For NH, the expectations are much more pessimistic, since the SHiP measurement would have a hard time to pin down φ 1 , even if δ is known, and therefore two unknowns φ 1 , θ would remain. They could be determined in principle from a measurement of Y B and |m ββ | but not from one single measurement and therefore the baryon asymmetry will be  Figure 11. Posterior probabilities in the |m ββ | vs Y B plane from a putative measurement at SHiP, assuming 0.1%, 1% uncertainties in the masses and mixings (red) or the latter with an additional measurement of δ in DUNE and HyperK (blue). The grey band is the standard 3ν expectation.
very difficult to predict in this case. The measurement of |m ββ | for NH would nevertheless be futuristic since the value would be roughly a factor 10 smaller, which is beyond the reach of the next generation of neutrinoless double beta decay experiments.

U PMNS phases from SHiP and neutrinoless double beta decay
We have seen that the ratios of electron and muon mixings that could be measured at SHiP could give very interesting information of the phases of the PMNS matrix. We stress that this is independent of whether the baryon asymmetry can be explained or not. This relies on the fact that the mixings are large enough so that they can be discovered at SHiP. In this case, the ratio of the electron to muon mixings are dominantly controlled by the phases of the PMNS mixing matrix, as can be seen in eqs. (4.7) and (4.8).
In Fig. 13, we compare the correlation expected on the (φ 1 , δ) plane from a putative measurement of this ratio using the analytical formulae of eqs. (4.7) for the test point of the previous section, to the posterior probabilities obtained with the most competitive SHiP uncertainties assumed in the previous section. Clearly the analytical formulae work very well and demonstrate the potential of SHiP in constraining the CP violating phases of the mixing matrix.
This has the following interesting consequence. If neutrinoless double beta decay could be measured with sufficient precision and the effect of the unknown θ would be small (this would happen for example in a more degenerate situation or for larger heavy masses which suppress the heavy contribution to neutrinoless double beta decay), the combination of this  what is the reach of the combination of SHiP and neutrinoless double beta decay on δ is very interesting and deserves a dedicated study. Reversely, if a measurement of δ could be obtained from neutrino oscillation measurements, a putative measurement of SHiP could provide a more precise prediction of the neutrinoless double beta decay amplitude if the neutrino ordering is inverted. This would be an extra motivation for improving the nuclear matrix element determination.
The possibility of measuring also the tau mixing at SHiP could help to resolve the degeneracy. Fig. 14 shows the isocurves of |U e4 | 2 /|U µ4 | 2 and |U e4 | 2 /|U τ 4 | 2 in the (φ 1 , δ) plane. The test point we used did not have sensitivity to the tau mixing according to the expectations for SHiP, but if this measurement could be improved this would also be useful to reduce the (δ, φ 1 ) degeneracy.  Figure 14.

Conclusions
We have studied the production of the matter-antimatter asymmetry in low-scale O(GeV) seesaw models. We have improved our previous study [15] by including the washout processes from gauge interactions and Higgs decays and inverse decays, quantum statistics in the computation of all rates, as well as spectator processes. This together with a more efficient numerical treatment of the equations have allowed us to perform the first bayesian exploration of the full parameter space that explains the observed baryon asymmetry in the context of the minimal model, where only two singlets play a role in the generation of the baryon asymmetry.
We have demonstrated that successful baryogenesis is possible with a mild heavy neutrino degeneracy, and more interestingly that these less fine-tuned solutions prefer smaller masses M i ≤ 1GeV, which is the target region of SHiP, and significant non-standard contributions to neutrinoless double beta decay. We have also demonstrated the complementarity of future putative measurements from SHiP, neutrinoless double beta decay and searches for leptonic CP violation in neutrino oscillations, in the quantitative prediction of the baryon asymmetry within the minimal model. If singlets with masses in the GeV range would be discovered in SHiP and the neutrino ordering is inverted, the possibility to predict the baryon asymmetry (maybe up to a sign) looks in principle viable, in contrast with previous beliefs. Unrelated to whether the baryon asymmetry is explained or not, we have also shown that a measurement of the electron and muon mixings of heavy species within SHiP range could precisely determine, in the minimal model, a combination of the two phases of the U PMNS matrix.