Leptogenesis in $SO(10)$ Models with $A_4$ Modular Symmetry

We study the prediction for leptogenesis in two renormalizable supersymmetric $SO(10)\times A_4$ modular models in which the neutrino mass is dominantly generated by the type I seesaw mechanism. The evolution of the lepton asymmetries are described in terms of the three-flavored density matrix equations for three heavy Majorana neutrinos, where both vanishing initial condition and thermal initial condition of the right-handed neutrinos are considered. We also present an analytical approximation based on the Boltzmann equations. We find regions of parameter space compatible with the measured fermion masses and mixing parameters as well as the baryon asymmetry of the Universe. The predictions for the light neutrino masses, the effective mass in neutrinoless doble beta decay and the leptonic CP violation phases are discussed.


Introduction
The origin of neutrino mass and mixing remains one of the most important questions in particle physics, along with other unresolved questions such as the matter-antimatter asymmetry of the Universe, and the possible unification of the strong and electroweak interactions.
The excess of matter over antimatter in the Universe has been firmly established. Although the standard model (SM) qualitatively satisfies the Sakharov conditions [1], the CP violation in SM is too small to explain the baryon asymmetry. A dynamical generation of baryon asymmetry requires new physics beyond SM. The discovery of neutrino masses and mixing in neutrino oscillation experiments is a great progress in particle physics. The smallness of neutrino masses can be naturally explained in the type I seesaw mechanism in which a number of heavy Majorana neutrinos are added to the SM, and a direct cosmological implication is leptogenesis [2]. The observed baryon asymmetry can be produced by the out-of equilibrium, CP and B −L violating decays of the right-handed neutrinos, see Refs. [3,4] for recent review papers. In the type I seesaw model with three right-handed neutrinos, 18 additional parameters are introduced and the leptogenesis can successfully generate the observed baryon asymmetry if the neutrino Yukawa couplings and right-handed neutrino masses are freely chosen. It is difficult to test leptogenesis due to the large number of free parameters of in a generic seesaw model.
Masses and mixing angles of elementary fermions are known with good precision at present. The mass eigenvalues of quarks and leptons have a large hierarchy, the flavor mixings in quark and lepton sectors are drastically different from each other: small quark mixing angles versus large lepton mixing angles. Understanding the origin of the flavor structure of SM is another of the greatest challenges in particle physics. The fundamental principal behind the observed patterns of fermion masses and mixing is still elusive, and the approach of flavour symmetry has been extensively studied in the last several decades [5,6]. Many candidate groups for flavor symmetry have been discussed, and many models which can describe some pieces of the flavor puzzle were constructed. The flavor symmetry relates three generations of matter fields so that the resulting quark and lepton mass matrices have certain structure and the total number of free parameters in fermion mass matrices are reduced.
The flavor symmetry can not be exact symmetry and it has to be broken explicitly or spontaneously to accommodate the non-trivial mixing patterns of quarks and leptons. As a consequence, the fermion mass matrices are usually expanded as power series of the symmetry breaking parameters. The breaking of flavour symmetries typically relies on several scalar multiplets called flavons whose vacuum expectation values should be oriented along certain directions in flavor space. In addition, extra cyclic group generally is necessary to forbid the unwanted operators and realize the correct vacuum alignment [6]. Thus the construction of flavon potential is intricate and make the flavor symmetry models rather complicated.
The origin of all quark and lepton masses and mixing may be related with each other. This suggests to combine the Grand unified theories (GUTs) with flavor symmetry. The electromagnetic, weak, and strong forces are merged into a single force at high energies in GUTs. The quark and lepton fields are embedded into one or few gauge multiplets, thus GUTs specifically predict relations among the fermion masses. Inclusion of flavor symmetry allows to predict the mixing parameters of quark and leptons. However, combining traditional flavor symmetry with GUTs also requires the introduction of flavons and vacuum alignment in order to break the flavor symmetry, see Ref. [65] for a review. Therefore there is a strong motivation for introducing modular symmetry in the context of GUTs such that the complication of flavor symmetry breaking can be removed and the predictive power of the GUTs models can be improved considerably. The possible combinations of SU (5) GUTs and modular symmetry groups have been discussed in the literature, and several modular SU (5) GUT models have been built at level 2 [10,66], level 3 [14,67,68] and level 4 [69][70][71].
Besides the simple SU (5) GUT, the SO(10) GUT is another prototype of grand unified theory. The unification of matter is even more complete in SO(10) GUT, since the 16-dimensional spinor representation of SO(10) can accommodate all the known SM chiral fermions as well as the right handed neutrino. This makes SO(10) GUTs particularly attractive candidates for explaining the origin of neutrino mass and mixing, the unification of the strong and electroweak interactions, and the matter-antimatter asymmetry of the universe, simultaneously. It is notable that both the righthanded neutrino mass matrix and neutrino Yukawa coupling are constrained by the SO(10) GUT symmetry, and they are related with the quark and lepton mass matrices. The SO(10) GUTs employing an conventional flavor symmetry has been studied so far. For example, models with SO(10) × A 4 [72][73][74][75] and SO(10) × S 4 [76][77][78][79][80] have already been investigated. As in SU (5) GUTs, the cumbersome vacuum alignment required in conventional flavor symmetry can be eliminated by the use of modular symmetry.
Following this line of reasoning, SUSY SO(10)×A 4 modular models were constructed in [81]. We found that the models involving in addition to the Higgs fields in the 10 and the 126, also a Higgs field in the 120, proved to give a successful description of quark and lepton (including neutrino) masses and mixing, with many such models being found with sums of modular weights of up to 8 or less. The neutrino masses were generated by the type-I and/or type II seesaw mechanisms [82][83][84][85][86] arising from the SU (2) L singlet and/or triplet components of the Higgs in the 126 representation. However the question of leptogenesis was not addressed for these models.
In this paper, we study the prediction for leptogenesis in two of the renormalizable SUSY SO(10) × A 4 modular models from [81] in which the neutrino mass is dominantly generated by the type I seesaw mechanism 1 . The evolution of the lepton asymmetries are described in terms of the three-flavored density matrix equations for three heavy Majorana neutrinos, where both vanishing initial condition and thermal initial condition of the right-handed neutrinos are considered. We also present an analytical approximation based on the Boltzmann equations. We find regions of parameter space compatible with the measured fermion masses and mixing parameters as well as the baryon asymmetry of the Universe. The predictions for the light neutrino masses, the effective mass in neutrinoless doble beta decay and the leptonic CP violation phases are discussed. This paper is organized as follows. In section 2 we briefly review the modular invariant SO(10) models and the predictions for the fermion mass matrices. In section 3 we present the formalism of density matrix to describe the evolution of the lepton asymmetries and the number densities of right-handed neutrinos. Moreover, the decay and washout of three right-handed neutrinos are considered in turns in the framework of Boltzmann equations, and an approximation formula for the final lepton asymmetry is given. In section 4 we present our numerical analysis for the concerned benchmark models, the results of the fit and the predictions of the models are given. Section 5 concludes the paper. In the appendix A we give two examples where the analytical approximations disagree with the numerical results.

Fermion masses in SO(10) GUTs with A 4 modular symmetry
It is known that all the fermionic multiplets of each SM generation, plus one right-handed neutrino singlet, fit exactly into the 16 dimensional spinorial representation of the SO(10) grand unification group. As a consequence, in renormalizable SO(10) models, the Higgs fields that contribute to fermion masses are in the SO(10) representations 10, 126 and 120 which are denoted by H, ∆ and Σ respectively. The most general form of the Yukawa superpotential for renormalizable SO(10) is given by with i = 1, 2, 3, and ψ i stand for the three generations of matter fields which are in the 16 dimensional representation of SO (10). The Yukawa coupling matrices Y 10 and Y 126 are symmetric in generation space, whereas Y 120 ij are antisymmetric. Their explicit forms can be constrained by the flavor symmetry which relates the three generations of fermions. The GUT Higgs H and ∆ contain one pair of up-type and down-type Higgs doublets while Σ has two pairs of such doublets. It is assumed that only one linear combination of up-type and one of down-type Higgs doublets are light, and the remaining combinations acquire GUT scale masses. The vacuum expectation value (VEV) of the light Higgs doublets breaks the electroweak symmetry and generate the fermion masses. The mass matrices of quarks and charged leptons can be written as [81,[87][88][89][90][91] where v u and v d refer to the VEVs of the MSSM Higgs fields H u and H d respectively, and some ratios of VEVs have been absorbed into the Yukawa matrices Y 10 , Y 126 and Y 120 [81]. Moreover, the Higgs multiplet ∆ also contains the Pati-Salam multiplets (10, 3, 1) and (10, 1, 3) which can generate the Majorana masses for the left-handed and right-handed neutrinos. Hence generally the light neutrino mass receives contributions from both type I and type II seesaw mechanism, with The dimensionless parameters r 1,2,3 and c e,ν in Eqs. (2,4) are determined by the Clebsch-Gordan coefficients of SO(10) and the mixing among the Higgs fields. In the present work, we consider the scenario of type I seesaw dominance which is the limit of v L → 0.
At weight 6, we have three modular multiplets which transform in 1 and 3 of A 4 , Although the SO(10) GUT symmetry unifies all the fermions of each generation into a single representation 16, the fermions in different generations are not related so that the Yukawa coupling matrices Y 10 , Y 126 and Y 120 in Eq. (1) can be arbitrary complex symmetric and antisymmetric ψ 1,2,3 H Σ ∆ SO(10) 16 10 120 126 matrices respectively. In order to understand the observed flavor structure of quarks and leptons, the A 4 modular symmetry is imposed in the SO(10) GUT. We assign the three generations of matter fields ψ 1,2,3 to a A 4 triplet since the flavor symmetry would effectively be the Z 3 subgroup generated by T rather than A 4 for the singlet assignment. The GUT Higgs multiplets H, ∆ and Σ all transform as 1. Taking into account the A 4 modular symmetry, the SO(10) Yukawa coupling of Eq. (1) becomes where k F , k 10 , k 126 and k 120 are the modular weights of ψ, H, ∆ and Σ respectively. Notice that one has to sum over the contributions of all independent modular multiplets at the relevant weights.
Using the multiplication rules in Eq. (6), the Yukawa matrices Y 10 can be read off as follows, and Y 126 takes a similar form with α i replaced by γ i . Because the coupling with Σ is antisymmetric in the generation indices, only the antisymmetric triplet contraction (ψψ) 3 A contributes to the Yukawa matrix Y 120 ,

Benchmark models
In the following, we present two typical SO(10)×A 4 modular models given in our previous work [81], then study the predictions for leptogenesis. The first benchmark model denoted as BM1 is characterized by the modular weights (2k F + k 10 , 2k F + k 120 , 2k F + k 126 ) = (4, 2, 4), and it is the non-minimal model 2 of [81]. Thus the modular invariant superpotential is given by The second model denoted as BM2 is specified by the modular weights (2k F + k 10 , 2k F + k 120 , 2k F + k 126 ) = (4, 6, 4), and it is the non-minimal model 3 of [81]. We can straightforwardly read out the superpotential as follows, In the case that the light neutrino mass is dominated by the type-I seesaw contribution, the best fit values of the free parameters for the above two models are listed in table 3 and table 4, and we see that the experimental data can be accommodated very well.

Prediction for baryon asymmetry via leptogenesis
The baryon asymmetry of the Universe is a puzzle in particle physics, and its value can be inferred from observations in two independent ways. The first is from the big bang nucleosynthesis. With the assumption of only three light neutrinos, the predictions for the abundances of the light elements D, 3 He, 4 He, and 7 Li depend on a single parameter, that is the difference between the number of baryons and anti-baryons normalized to the number of photons: where the subscript 0 means "at present time". The range of η which is consistent with the abundances of D and 4 He is determined to be [92] η = (6.143 ± 0.190) × 10 −10 , at 1σ confidence level. This range corresponds to and it is related to the parameter η through Y ∆B η/7.04. Furthermore the baryon asymmetry produces acoustic oscillations in the power spectrum of the cosmic microwave background [93]. Observing these oscillations gives an even tighter bound on the value of baryon asymmetry, It is impressive that the above two unrelated measurement approaches give consistent range of baryon asymmetry. In this work, we will use the conservative value of Eq. (20) from the big bang nucleosynthesis. Generally speaking, both type I and type II seesaw mechanisms contribute to the light neutrino mass in SO(10) GUT, the contribution of type I seesaw is assumed to be dominant over that of type II in the present work. Then the CP violating out-of-equilibrium decay of the heavy Majorana neutrinos can generate a lepton asymmetry which is subsequently partially converted to a baryon asymmetry through the B + L violating sphaleron processes [2]. In this section, we shall study whether the measured value of the baryon asymmetry of the Universe Y ∆B = (8.72±0.27)×10 −11 in Eq. (20) can be correctly generated through thermal leptogenesis in the above benchmark SO(10) × A 4 modular models. The evolution of the lepton asymmetry is described by the density matrix equation.

Density matrix equation
In the simplest scenario of leptogenesis, the lepton and anti-lepton final states from the decays of the heavy neutrino N i are in coherent superposition of flavor, where i = 1, 2, 3 and α = e, µ, τ , and the projection coefficients at tree-level are given by Notice that we work in the basis where both charged lepton and right-handed neutrino mass matrices are diagonal. This is a good approximation for high temperature T 10 12 GeV when the charged lepton Yukawa interactions are negligible. However, the coherent evolution of the states | i and |¯ i would be broken down and different lepton flavors are distinguishable when these interactions are in thermal equilibrium. Then the left-handed leptons can be rapidly converted to right-handed components through scattering with the Higgs doublet. To be more specific, the τ and µ Yukawa interactions come into thermal equilibrium at the temperature T ∼ 10 12 GeV and T ∼ 10 9 GeV respectively.
It has been shown that the density matrix equation can provide an accurate description for the time evolution of the lepton asymmetry [94][95][96][97][98][99][100]. The approach of the density matrix accounts for the quantum decoherence effect and it allows to quantitatively describe the transitions between the one-flavored, two-flavored and three-flavored regimes. Comprehensive analysis of leptogenesis has been performed in the framework of density matrix [101][102][103]. The formalism of the density matrix for leptogenesis has been given in Ref. [100], we shall follow closely the notation of [100] in the following. For three decaying heavy Majorana neutrinos, the most general form of the density matrix equations are given by which describe how a given flavor of lepton is washed out. The quantity N N i is the number of the heavy neutrino N i in a portion of comoving volume which contains one photon at the temperature z 0. Hence the equilibrium number density N eq N i for the Boltzmann statistics is given by and K 2 (z) the modified Bessel functions of the second kind.
are the known number densities for the B/3 − L α asymmetry, and the off-diagonal elements N B−L αβ describe the degree of coherence between the flavor states. The rescaled decay rate D i of the right-handed neutrino N i is given by where H is the Hubble expansion rate. The washout parameters W i read Furthermore, the last two terms in Eq. (24) describe the effect of charged lepton Yukawa interactions which can induce the transition of left-handed leptons to right-handed leptons. They describe decoherence via the damping of the off-diagonal terms of the B − L asymmetry matrix N B−L and they are determined as where M Pl 1.22 × 10 19 GeV is the Planck mass, the Yukawa couplings y µ 6.07 × 10 −4 and y τ 1.02 × 10 −2 are fixed by the µ and τ masses m µ 105.66 MeV and m τ 1776.86 MeV respectively. When the temperature goes below 10 12 GeV, the τ Yukawa interactions come into thermal equilibrium and leads to decoherence of τ lepton states. This implies the transition from an unflavoured to two-flavoured regime. When the temperature drops below 10 9 GeV the similar effects arise from the µ Yukawa interactions. Analogously, the electron Yukawa dependent damping term should be considered for M 1 < 10 5 GeV. Finally the CP asymmetry ε (i) αβ generated by N i decay is of the following form [95,98,100,104,105] with .
The diagonal components ε (i) αα of the CP asymmetry matrix are exactly the usual flavored CP asymmetries, which is then converted into the baryon asymmetry of the Universe, η = 28 79 where the factor 28/79 accounts for the partial conversion of the B − L asymmetry into the baryon asymmetry by sphaleron process, and the second factor 1/27 describes the dilution of the baryon asymmetry due to the change of the photon density between leptogenesis and recombination. In the following numerical analysis, we will use the Python package ULYSSES [106] to solve the coupled density matrix equations shown in Eq. (24).

Analytic approximation within the Boltzmann equations
In our benchmark models, the right-handed neutrino masses are hierarchical with M 3 M 2 10 12 GeV while 10 9 GeV M 1 10 12 GeV. Therefore the decays of N 3 and N 2 are in the singleflavored regime while N 1 decay is in the two-flavored regime. In the following, we shall analyze the sequential decay of the three heavy Majorana neutrinos and evolution of lepton symmetry in terms of Boltzmann equations. In the first stage around the temperature T ∼ M 3 , a B − L asymmetry is generated from the N 3 decay. In the second (third) stage at T ∼ M 2 (M 1 ) another B − L asymmetry is generated from the N 2 (N 1 ) decays, while the inverse decay process of N 2 (N 1 ) become effective and washout the asymmetry of N 3 decay to some level. In each stage we can use the Boltzmann equations to obtain an approximate analytical expression for the B − L asymmetry.
In the first stage for M 3 T is the freeze-out temperature of the N 3 inverse decay and z B 3 2 + 4K 0.13 3 e −2.5/K 3 = O(1 − 10) [107], the rates of the τ and µ Yukawa interactions are much smaller than the expansion rate of the Universe. Consequently the charged lepton lepton state produced from N 3 decay is in coherent superposition and the flavor states are indistinguishable. In the single flavor approximation, the time evolution of the number densities of N 3 and B − L can be described by the following semi-classical Boltzmann equations: The solution to the above equations is given by [108] where κ(K 3 ) is the so-called efficiency factor and its value depends on the initial condition of decaying heavy neutrinos. For thermal initial abundance of right-handed neutrinos, it is well approximated by [108], For the vanishing initial abundance of right-handed neutrinos, the approximation formula for the efficiency factor is [108] with The lepton state | 3 produced by N 3 decay can be regarded as coherent superposition of a 2 parallel component and of a 2 orthogonal component which is denoted as | 2 . To be more concrete, | 2 is the projection of | 3 on the plane orthogonal to | 2 such that 2 | 2 = 0, and its explicit expression is where the probabilities p 32 ≡ | 2 | 3 | 2 and p 32 ≡ | 2 | 3 | 2 with p 32 + p 32 = 1. From Eq. (22), we can read out the general expression of p ij as Correspondingly the same decomposition can be made for the asymmetry N ∆ 3 , These two asymmetries can be used as initial condition at the beginning of the N 2 production and decay. When the temperature goes down to T ∼ M 2 , the N 2 decay and inverse processes break the coherent evolution of | 3 which becomes an incoherent mixture of | 2 and | 2 . The | 2 component undergoes the washout from N 2 while the orthogonal component | 2 is unwashed. The kinetic evolution equations are given by Hence the asymmetries at the temperature T T B 2 take the form, Since 10 9 GeV M 1 10 12 GeV, the τ Yukawa interaction becomes effective before the onset of N 1 washout processes. Let us indicate with M 1 T 10 12 GeV that value of temperature below which one can approximate the coherence of |τ and the orthogonal component completely damped. Thus the lepton quantum state can be treated as an incoherent mixture of three components: |τ , |τ 2 and |τ2 , where |τ 2 and |τ2 are projections of | 2 and | 2 on the e − µ plane, respectively. Hence the lepton states |τ 2 and |τ2 read as The corresponding probabilities are Thus the B − L asymmetry before N 1 decay can be decomposed into the following three parts: The asymmetries in the tau flavor and the orthogonal components |τ 1 and |τ1 experience different washout from inverse decay, The quantum state |τ 1 is the projection of | 1 on the e − µ plane and |τ1 is the quantum state orthogonal to |τ 1 in e − µ plane, with p 1τ = |C 1τ | 2 . One can use the two-flavored Boltzmann equations to describe the time evolution of the asymmetries during N 1 leptogenesis, where ε 1τ 1 ≡ ε 1e + ε 1µ . Thus the frozen values of the asymmetries at the end of N 1 washout are given by where K 1τ 1 ≡ K 1e + K 1µ . Notice that decomposing N ∆τ 2 (T T ) and N ∆τ2 (T T ) along the directions of |τ 1 and |τ1 gives rise to N ∆τ 1 (T T ) and N ∆τ1 (T T ) as follows, with We would like to mention that the factors C 1α , C 2α , C 3α depend on the neutrino Yukawa, as shown in Eq. (23). The total final asymmetry is the sum of the three terms It turns out that our benchmark models are in the strong washout regime with K 3 , K 2 , K 1τ , K 1τ 1 1, thus we neglect the contributions which undergo the washout exponential suppression of at least one right-handed neutrino. The final lepton asymmetry will be then dominated by the unwashed terms and it is approximately given by Notice that this approximation formula is applicable to both thermal and vanishing initial abundance of right-handed neutrinos, and the corresponding efficiency factors are given in Eq. (37) and Eq. (38) respectively.

Numerical results
In this section, we shall perform a detailed numerical analysis to explore whether the renormalizable SO(10)×A 4 modular models of section 2.2 can explain the experimentally measured values of masses and mixing parameters of both quarks and leptons as well as the matter-antimatter asymmetry of the Universe. We perform a χ 2 analysis to optimize the values of the free parameters. The χ 2 function is defined as   [109], with the SUSY breaking scale MSUSY = 500 GeV and tan β = 10, where the error widths represent 1σ intervals. The values of lepton mixing angles, leptonic Dirac CP violation phases δ l CP and the neutrino mass squared difference are taken from NuFIT 5.1 [110] for normal ordering neutrino masses with Super-Kamiokande data.
Here µ i and σ i denote the experimental central values and the 1σ uncertainties of the observables respectively, and their values are obtained by evolving their low energy values to the GUT scale with the renormalization group equations, as shown in table 2. For the neutrino oscillation data, we assume normal ordering neutrino mass spectrum so that the renormalization group induced corrections to the neutrino masses and mixing parameters can be negligible for small tan β, and we use the most recent results from the NuFit collaboration [110].
The mechanism of modulus stabilization is still an open question, consequently we treat the complex modulus τ as a free parameter in the fundamental region D = τ ∈ C Imτ > 0, |Reτ | ≤ 1 2 , |τ | ≥ 1 , which can represent all possible values of τ in the upper half complex plane up to a modular transformation. In the numerical analysis, the magnitudes of all coupling constants are limited in the region [0, 10 4 ] and their phases are freely varied in the range [0, 2π]. The overall scales α 1 v u , α 1 r 1 v u and α 2 1 v 2 u /v R are fixed by the experimentally measured values of top quark mass, bottom quark mass and solar neutrino mass squared difference ∆m 2 21 . For each set of given values of the input parameters, the quark and lepton mass matrices in Eqs. (2,3) can be straightforwardly diagonalized and subsequently mixing angles and CP violation phases can be extracted in the usual way. We numerically minimize the χ 2 function by using the minimization algorithms TMinuit [111] developed by CERN to determine the optimum values of the input parameters. In the present work, we will focus the scenario of type I seesaw dominance, and the out of equilibrium decay of right-handed neutrinos is used to generate the baryon asymmetry of the Universe.
For each model, we present the point in parameter space which minimizes the χ 2 , as a result of a numerical minimization procedure. We shall consider two typical initial conditions of heavy Majorana neutrinos: thermal initial abundance and vanishing initial abundance. We show the numerical results in table 3 and table 4 for vanishing and thermal initial conditions respectively. It is notable that the value of τ is close to the residual modular symmetry fixed point τ ST ≡ −1/2+i √ 3/2. Since both models at the best points are in the strong washout regime, the prediction for the baryon asymmetry is not so sensitive to the initial condition of right-handed neutrinos. Thus the values of each input parameter have little difference for the concerned two initial conditions. We also present the predictions for the mass ratios of charged fermions m e /m µ , m µ /m τ , m u /m c , m c /m t , m d /m s , m s /m b and m b /m τ , the mixing angles θ l 12 , θ l 13 , θ l 23 , θ q 12 , θ q 13 and θ q 23 as well as the quark CP violation phase δ q CP , which are compatible with experimental data. Moreover, from the fitted values of the parameters, we can derive the predictions for the still unmeasured observables including the leptonic Dirac CP violation phase δ l CP and Majorana CP phases α 21 and α 31 , the light neutrino masses m 1,2,3 , the effective Majorana neutrino mass m ββ for neutrinoless double decay and the right-handed neutrino masses M 1,2,3 at the best fit point. The right-handed neutrino masses have a hierarchical pattern M 3 M 2 M 1 which are in the range of 10 10 − 10 13 GeV. The effective Majorana mass m ββ is determined to be 0.532 (1.059) meV and 0.622 (0.550) meV in the models BM1 and BM2 respectively for vanishing (thermal) initial condition, they are much below the current most stringent bound m ββ < (61 − 165) meV given by the KamLAND-Zen Collaboration [112]. Future tonne-scale neutrinless double decay experiments can reduce the sensitivity on m ββ to few meV such as m ββ < (4.7 − 20.3) meV from nEXO [113] and m ββ < (9 − 21) meV from LEGEND-1000 [114]. Consequently our predictions for m ββ are even out of reach of tonne-scale detectors. As shown in tables 3 and 4, we have split the χ 2 function into four different contributions χ 2 total = χ 2 l + χ 2 q + χ 2 bτ + χ 2 Y B , where χ 2 l , χ 2 q , χ 2 bτ and χ 2 Y B stand for the contributions induced by the deviations of the lepton sector observables, the quark sector observables, the mass ratio m b /m τ and the baryon asymmetry Y B from their central values respectively. We see that the largest deviation arises from the quark sector in both models, and χ 2 total is dominated by χ 2 q although the predictions for the quark masses and mixing parameters are in the experimentally allowed region.
Regarding the baryon asymmetry, the concrete form of the neutrino Yukawa coupling matrix and right-handed neutrino masses are fixed at the best fit point, we numerically solve the density matrix equation of Eq. (24) with the help of ULYSSES to obtain the frozen value of B − L asymmetry. It is remarkable that the numerical results agree well with the analytical approximations. We see that the observed matter-antimatter asymmetry can be accommodated. In general, the decays of three righthanded neutrinos N 1,2,3 all contribute to the final baryon asymmetry, and asymmetries generated by N 2,3 are partially washed out by the N 1 process. In the model BM1 for both initial conditions and the BM2 model with vanishing initial abundance, the leptogensis is dominated by N 1 decay. The contributions of N 2 and N 3 are as important as the that of N 1 in the BM2 for thermal initial abundance. Moreover, we plot the evolution of Y B with the temperature T for the two benchmark models in figure 1, and the corresponding numerical results of Boltzmann equations are displayed for comparison.
Furthermore, we comprehensively scan the parameter space around the best fit points of table 3  and table 4, and we require all the fermion masses and mixing parameters are in the experimentally preferred 3σ regions. The experimentally allowed values of the complex modulus τ and the correlations between the masses and mixing parameters are shown in figures 2, 4, 6 and 8 for the benchmark models with vanishing and thermal initial conditions. The corresponding predictions for the baryon asymmetry Y B with respect to sin 2 θ l 23 , δ l CP , α 21 , α 31 are displayed in figures 3, 5, 7 and 9. We would like to mention that the contribution of Y B to χ 2 is not included in these figures. We see that imposing successful leptogenesis yields specific predictions for the fermion masses and mixing parameters.

Conclusion
In this paper we have studied the prediction for leptogenesis in two renormalizable supersymmetric SO(10) × A 4 modular models in which the neutrino mass is dominantly generated by the type I seesaw mechanism. The evolution of the lepton asymmetries are described in terms of the three-flavored density matrix equations for three heavy Majorana neutrinos, where both vanishing initial condition and thermal initial condition of the right-handed neutrinos are considered. We also presented an analytical approximation based on the Boltzmann equations. We found regions of parameter space compatible with the measured fermion masses and mixing parameters as well as the baryon asymmetry of the Universe. The predictions for the light neutrino masses, the effective mass in neutrinoless doble beta decay and the leptonic CP violation phases were discussed.
The right-handed neutrinos naturally appear in SO(10) GUT, all the chiral fermions of one generation plus one additional right-handed neutrino are unified into a single 16 dimensional spinor representation of SO (10). Thus SO(10) GUT provides a natural explanation for the tiny neutrino mass through the type I seesaw mechanism. From the view of cosmology, it is well known that the CP violating decay of the right-handed neutrinos can generate the matter-antimatter asymmetry of the Universe, this is the so-called leptogenesis. Moreover, the mass and Yukawa coupling of the right-handed neutrinos are closely related to those of quarks and charged leptons, and their concrete values can be determined by requirement that the measured masses, mixing angles and CP violation    phases of both quarks and leptons should be accommodated. It is attractive that the fermion masses and mixing together with baryon asymmetry of the Universe can be explained in a given SO(10) model. Modular symmetry is an important progress to address the flavor puzzle of the SM, the vacuum alignment problem in conventional flavor symmetry models is greatly simplified, and the complex modulus τ could be the unique source of flavor symmetry breaking. As a consequence, the modular symmetry models could be quite predictive. In the present work, we have analyzed the prediction for leptogenesis in two benchmark renormalizable SO(10) models with A 4 modular symmetry where the neutrino masses dominantly arise from type I seesaw mechanism. The three generations of matter fields are assigned to triplet of A 4 modular symmetry, and the Higgs fields H, Σ, ∆ in the 10, 120 and 16 dimensional representations of SO(10) are invariant under A 4 . As a consequence, the structure of the SO(10) × A 4 modular models are fully specified by the modular weights of the fermion fields and Higgs fields, as shown in table 1. The modular weights are (2k F + k 10 , 2k F + k 120 , 2k F + k 126 ) = (4, 2, 4) and (2k F +k 10 , 2k F +k 120 , 2k F +k 126 ) = (4, 6, 4) respectively in the two concerned benchmark models.
We used the density matrix equations to describe the kinetic evolution of the lepton asymmetry. In comparison with Boltzmann equations, the density matrix equations give a quantitative description of lepton flavor effects in leptogenesis, and allows to precisely calaulate the final lepton asymmetry for arbitrary right-handed neutrino mass spectrum. The right-handed neutrino masses are hierarchical with M 3 M 2 10 12 GeV and 10 9 GeV M 1 10 12 GeV in our benchmark models. Hence the decays of N 3 and N 2 are in the single-flavored regime while N 1 decay is in the two-flavored regime. In order to understand the physical process of leptogenesis, we assume the three heavy Majorana neutrinos decay in sequence. The lepton asymmetries generated from the decay of each right-handed neutrino and the washout effects are governed by the Boltzmann equation, the analytical approximate solutions are given. The approximation agrees with the numerical results of density matrix equations well if the interplay of the different right-handed neutrinos and the interference of lepton flavor are insignificant.
The superpotential of the Yukawa is strongly constrained by modular symmetry, it involves 7 couplings α 1,2,3 , β 1 and γ 1,2,3 in the BM1 model, and an additional parameter β 2 in the BM2 model. Moreover, the quark and lepton mass matrices depend on 5 parameters r 1,2,3 , c e,ν which describe the light Higgs combinations. Notice that the Higgs sector of the renormalizable SO(10) models with all 10, 120 and 126 dimensional Higgs multiplets are complicated and it reduces the predictive power of SO(10) models. The 126-plet Higgs ∆ contains a right-handed scalar triplet ∆ R whose vacuum expectation value ∆ R = v R gives mass to the right-handed neutrinos. The parameters α 1 , r 1 and v R can be taken real without loss of generality. Including the complex modulus τ , the BM1 and BM2 models depend on 25 and 27 real parameters respectively. We fit all these free parameters using the 20 observables in table 2. In particular, we take into account the matter-antimatter asymmetry of the Universe through leptogenesis besides the measured masses and mixing parameters of quarks and leptons. As shown in table 3 and table 4, the experimental data can be well accommodated and it is remarkable that all the coupling constants are order one. Furthermore, we can obtain predictions for unmeasured quantities such as the Majorana CP phases α 21 , α 31 , the effective mass m ββ of neutrinoless double decay and the right-handed neutrino masses M 1,2,3 . m ββ is predicted in the range 0.5-1.1 meV, it is far below the sensitivity of current and future experiments. Hence our models would be ruled out if a positive signal of neutrinoless double beta decay is reported in future.

A Two examples of analytical approximation failure
The formalism of density matrix allows for a more general description of leptogenesis than the semi-classical Boltzmann equations. The off-diagonal elements of the density matrix describe the degree of coherence between the flavour states, so that the decoherence effects are incorporated into the equations. Moreover, the density matrix equations can describe the dynamical process of the transition between flavour-regimes, and allows to calculate the lepton asymmetry in the intermediate regimes where the single-flavored and two-flavored treatments are inadequate. Hence it is found that the Boltzmann equations fail to describe correctly the generation of baryon asymmetry in many cases [101][102][103]. Although the analytical approximation agree well with the numerical results of the density matrix equation in section 4, we note that this not always true. In particular, when the effects of decoherence cannot be neglected or there are non-trivial interplay between the decays and the washout processes of the three right-handed neutrinos. We give two such benchmark points in table 5 at which our analytical estimates fail to give the correct value of baryon asymmetry, and even the sign of Y B is wrong.
For the model BM1 with the values of the free parameters given in table 5, the decay parameters of N 2 and N 3 are determined to be K 2 = 90 and K 3 = 98 respectively. Therefore the freezeout temperature of N 3 is approximately T B 3 M 3 /9 2M 2 , the effects of N 2 are not negligible before the freeze-out of N 3 inverse decay, as can be seen from the left-top panel of figure 10. In addition, the asymmetries generated by N 2 and N 3 are only partially washed out by the N 1 decay and inverse decay processes and in fact they dominate the final asymmetry. Moreover, the mass of N 1 is 2.308 × 10 11 GeV, consequently it is in the transition region between the single-and two-flavored leptogenesis. This might contribute to the failure of the analytical estimates, since the transitions between the different flavour regimes are not taken into account in the Boltzmann equations and the baryon asymmetry η could change its sign during the transitions [103]. Analogously, for another model BM2, there is also non-trivial interplay between the decay and the washout processes of the heavy right-handed neutrinos. The contribution of N 1 is found to be subdominant, and the N 1 mass is 8.367 × 10 9 GeV which is close to the transition scale between the two-flavor and three-flavor regimes. Furthermore, from figure 10 we see that the decoherence of flavor states is not infinitely fast as assumed in the Boltzmann equations. This can also lead to deviations of analytical approximation from the numerical results of density matrix.  Table 5: Benchmark points at which the analytical approximation of section 3.2 fails for the concerned models. The vanishing initial abundance of right-handed neutrinos is assumed in BM1 model while thermal initial abundance is assumed in BM2. Figure 2: The values of the complex modulus τ compatible with experimental data and the correlations between the neutrino mixing angles, CP violation phases, quark mass ratios and mixing parameters for the model BM1 with vanishing initial condition. The vertical and horizontal dashed lines are the 3σ bounds taken from [110]. Note that the contribution of YB to χ 2 is not included here. Figure 3: The numerical results for the correlation between the baryon asymmetry YB and sin 2 θ l 23 , δ l CP , α21, α31 in the BM1 model with vanishing initial condition. The horizontal band denotes the experimentally allowed region of YB. The star " " stands for the best fitting point in table 3. Note that the contribution of YB to χ 2 is not included here. Figure 4: The values of the complex modulus τ compatible with experimental data and the correlations between the neutrino mixing angles, CP violation phases, quark mass ratios and mixing parameters for the model BM2 with vanishing initial condition. The vertical and horizontal dashed lines are the 3σ bounds taken from [110]. Note that the contribution of YB to χ 2 is not included here. Figure 5: The numerical results for the correlation between the baryon asymmetry YB and sin 2 θ l 23 , δ l CP , α21, α31 in the BM2 model with vanishing initial condition. The horizontal band denotes the experimentally allowed region of YB. The star " " stands for the best fitting point in table 3. Note that the contribution of YB to χ 2 is not included here. Figure 6: The values of the complex modulus τ compatible with experimental data and the correlations between the neutrino mixing angles, CP violation phases, quark mass ratios and mixing parameters for the model BM1 with thermal initial condition. The vertical and horizontal dashed lines are the 3σ bounds taken from [110]. Note that the contribution of YB to χ 2 is not included here. Figure 7: The numerical results for the correlation between the baryon asymmetry YB and sin 2 θ l 23 , δ l CP , α21, α31 in the BM1 model with thermal initial condition. The horizontal band denotes the experimentally allowed region of YB. The star " " stands for the best fitting point in table 4. Note that the contribution of YB to χ 2 is not included here. Figure 8: The values of the complex modulus τ compatible with experimental data and the correlations between the neutrino mixing angles, CP violation phases, quark mass ratios and mixing parameters for the model BM2 with thermal initial condition. The vertical and horizontal dashed lines are the 3σ bounds taken from [110]. Note that the contribution of YB to χ 2 is not included here. Figure 9: The numerical results for the correlation between the baryon asymmetry YB and sin 2 θ l 23 , δ l CP , α21, α31 in the BM2 model with thermal initial condition. The horizontal band denotes the experimentally allowed region of YB. The star " " stands for the best fitting point in table 4. Note that the contribution of YB to χ 2 is not included here. Figure 10: The evolution of the baryon asymmetry YB with temperature for the local minimum shown in table 5. The left and right panels are the BM1 with vanishing initial abundance and the BM2 with thermal initial abundance respectively. The evolution of the B − L asymmetry number densities N B−L αβ , α, β = e, µ, τ displayed in the lower panel is described by the density matrix equations in Eq. (24).