Modular Invariant Models of Leptons at Level 7

We consider for the first time level 7 modular invariant flavour models where the lepton mixing originates from the breaking of modular symmetry and couplings responsible for lepton masses are modular forms. The latter are decomposed into irreducible multiplets of the finite modular group $\Gamma_7$, which is isomorphic to $PSL(2,Z_{7})$, the projective special linear group of two dimensional matrices over the finite Galois field of seven elements, containing 168 elements, sometimes written as $PSL_2(7)$ or $\Sigma(168)$. At weight 2, there are 26 linearly independent modular forms, organised into a triplet, a septet and two octets of $\Gamma_7$. A full list of modular forms up to weight 8 are provided. Assuming the absence of flavons, the simplest modular-invariant models based on $\Gamma_7$ are constructed, in which neutrinos gain masses via either the Weinberg operator or the type-I seesaw mechanism, and their predictions compared to experiment.


Introduction
The puzzle of quark and lepton masses and mixing, left unanswered by the Standard Model (SM), may be addressed by introducing some family symmetry, which is generally non-Abelian and may be associated with a finite discrete group. Modular symmetry has been suggested as the origin of such a flavour symmetry, with neutrino masses as complex analytic functions called modular forms [1]. Finite non-Abelian discrete family symmetry emerges from the modular symmetry at various positive integer levels, with each level associated with a particular flavour group. At each level, the physical fields carry various modular weights which do not have to add up to zero in the coupling terms of the effective Lagrangian since the effective Yukawa couplings may be modular forms, which are holomorphic functions of a complex modulus field τ [1]. This may allow flavon fields to be removed, with higher-dimensional operators in the superpotential being completely determined by modular invariance and supersymmetry. The neutrino masses and mixing parameters may be predicted in terms of a few input parameters, although the predictive power of this framework may be reduced by the Kähler potential which is less constrained by modular symmetry [2].
The finite modular groups Γ 2 ∼ = S 3 [3][4][5][6], Γ 3 ∼ = A 4 [1,3,4,[7][8][9][10][11][12][13][14], Γ 4 ∼ = S 4 [13,[15][16][17][18][19] and Γ 5 ∼ = A 5 [18,20,21] have been considered. For example, simple A 4 modular models can reproduce the measured neutrino masses and mixing angles [1,8,12]. The quark masses and mixing angles may also be included together with leptons in an A 4 modular invariant model [22], and it has been shown how natural fermion mass hierarchies can arise as a result of a weighton field [23]. The modular invariance approach has been extended to include odd weight modular forms which can be decomposed into irreducible representations of the the homogeneous finite modular group Γ N [24], and the modular symmetry Γ 3 ∼ = T has been discussed, including the new possibility of texture zeroes [25]. Also modular symmetry may be combined with generalized CP symmetry, where the modulus transforms as τ → −τ * under the CP transformation [26][27][28][29][30]. The formalism of the single modulus has been generalized to the case of a direct product of multiple moduli [31,32], which is motivated by the additional extra dimensions in superstring theory, assuming toroidal compactification. Indeed, from a top-down perspective, modular symmetry naturally appears in string constructions [27,[33][34][35][36]. It has been realised that, if the VEV of the modulus τ takes some special value, a residual subgroup of the finite modular symmetry group Γ N would be preserved. The phenomenological implications of the residual modular symmetry have been discussed in the context of modular A 4 [10,13], S 4 [13,16] and A 5 [20] symmetries. If the modular symmetry is broken down to a residual Z 3 (or Z 5 ) subgroup in charged lepton sector and to a Z 2 subgroup in the neutrino sector, the trimaximal TM1 and TM2 mixing patterns can be obtained [10,16]. Moreover, the dynamics of modular symmetry could potentially be tested in present and forthcoming neutrino oscillation experiments [37].
In this paper we consider the level 7 finite modular group Γ 7 ∼ = P SL(2, Z 7 ), the projective special linear group of two dimensional matrices over the finite Galois field of seven elements. This is the smallest simple discrete group which contains complex triplets and sextet representations. It contains 168 elements and is sometimes written as P SL 2 (7) or Σ(168) [38,39]. The relationship of this group to some other family symmetries that have been used in the literature is discussed in [40][41][42][43][44]. It has been proposed as a finite modular group in [45], whose notations we follow. Thus, we consider for the first time level 7 modular invariant flavour models based on Γ 7 where the lepton mixing originates from the breaking of modular symmetry and couplings responsible for lepton masses are modular forms. The latter are decomposed into irreducible multiplets of the finite modular group Γ 7 . At weight 2, there are 26 linearly independent modular forms, organised into a triplet, a septet and two octets of Γ 7 . A full list of modular forms up to weight 8 are provided. Assuming the absence of flavons, the simplest modular-invariant models are constructed in which neutrinos gain masses via either the Weinberg operator or the type-I seesaw mechanism, and their predictions compared to experiment.
We organise the rest of this paper in the following. In Section 2, we describe properties of the modular group Γ and its finite subgroup Γ 7 . A full list of modular forms of level 7 and weight up to 10 are derived in Section 3, based on the SageMath algebra system [46]. We construct a class of flavon-less lepton flavour models and study their experimental constraints in Section 4. Summary is given in Section 5. The group theory of Γ 7 is listed in Appendix A. We refer to Appendices B and C for alternative methods based on the Dedekind eta function method [1] and the theta function method [20], respectively. How to find an independent set of higher weight modular forms from hundreds of constraints is discussed in Appendix D.
2 Modular symmetry and modular forms of level N = 7 In the following, we briefly review the modular symmetry and the its congruence subgroups. The special linear group SL(2, Z) is constituted by 2 × 2 matrices with integer entries and determinant equal to one [47,48]: The SL(2, Z) group acts on the upper half plane H = {τ ∈ C | τ > 0} as the linear fractional transformation, It is easy to see the following identity (γ(τ )) = τ |cτ + d| 2 , γ = a b c d ∈ SL(2, Z) .
Consequently the image γ(z) ∈ H for any γ ∈ SL(2, Z) and τ ∈ H. It is obvious that and therefore we identify Hence the action of γ and −γ act on the complex modulus τ is exactly the same, and it is sufficient to consider the projective special linear group P SL(2, Z) = SL(2, Z)/{I, −I}, which is the quotient of SL(2, Z) by ±I. The group P SL(2, Z) is also called the modular group in the literature, it is a discrete group with infinite elements and it can be generated by two transformations S and T [47] which fulfill the relations S 2 = (ST ) 3 = 1 .
The actions of S and T on H are given by For a positive integer N , the principal congruence subgroup of level N is defined as which is a normal subgroup of the special linear group SL(2, Z). Obviously Γ(1) ∼ = SL(2, Z) is the special linear group. It is easy to obtain which implies T N ∈ Γ(N ), i.e., T N is an element of Γ(N ). Taking the quotient of Γ(1) and Γ (2) by {I, −I}, we obtain the projective principal congruence subgroups Γ(N ) = Γ(N )/{I, −I} for N = 1, 2, and Γ(N > 2) = Γ(N ) since the element −I does not belong to Γ(N ) for N > 2. The quotient groups Γ N = Γ(1)/Γ(N ) are usually called inhomogeneous finite modular groups, and the homogeneous finite modular group is defined as Γ N = SL(2, Z)/Γ(N ) which is the double covering of Γ N [24]. The finite modular group Γ N for N ≤ 5 can be obtained from Γ(1) by imposing the condition T N = 1. Consequently the generators S and T of Γ N obey the relations The groups Γ N with N = 2, 3, 4, 5 are isomorphic to the permutation groups S 3 , A 4 , S 4 and A 5 respectively [45]. Note again that for this group, as for all groups with N > 5, at least one additional relation is necessary in order to render the group finite. It is easy to calculate which implies Hence the element −(ST 3 ) 4 is belong to Γ (7). Notice that −(ST 3 ) 4 and (ST 3 ) 4 are identified as the same element of the Γ(1) group, since they lead to the same linear fraction transformations. Therefore the finite modular group Γ 7 of level N = 7 be generated by two generators S and T which satisfy the following multiplication rules 1 The crucial element of the modular invariance approach is the modular forms f (τ ) of weight k and level N . They are holomorphic functions of the complex modulus τ with well-defined transformation properties under the group Γ(N ), The modular forms of weight k and level N form a linear space M k (Γ(N )) of finite dimension.
In the present work, we shall focus on even weight modular forms, i.e., k being an even number. Then it is always possible to choose a basis of M k (Γ(N )) such that the modular forms transform according to a unitary irreducible representation r of Γ N [1], where γ is a representative element of Γ N , and ρ r (γ) is the representation matrix of γ in the irreducible representation r. If the modular weight k is odd, the modular forms can be decomposed into irreducible representations of the homogeneous finite modular group Γ N [24]. In order to determine the proper basis, it is sufficient to apply Eq. (16) to the generators S and T which can generate all elements of γ N .
3 Modular forms of level N = 7 The general dimension formula for the linear space of modular forms of level N and weight k is given by 1 The multiplication rule of Γ7 is [45].
For N = 7, we have Hence the linear space of modular forms of level 7 and weight 2 has dimension 14 × 2 − 2 = 26.
One can obtain q-expansions for a basis b i of the space of lowest weight modular forms for Γ 7 from the SageMath algebra system [46]. They are given by with q = e 2πi τ and where fractional powers q n/7 should be read as q n/7 = e 2nπi τ /7 . The above lowest weight modular forms can be organized into a triplet transforming in the representation 3 of Γ 7 , a septet transforming in the representation 7 of Γ 7 , and two octets in the 8 of Γ 7 . To be more explicit, we have Notice that the weight two modular multiplets Y (2) 8a (τ ) and Y (2) 8b (τ ) are not unique, and in principle they can taken to be any two linearly independent combinations of Y (2) 8a (τ ) and Y (2) 8b (τ ). Higher weight modular multiplets can be obtained from tensor products of the lowest weight multiplets Y 8a and Y (2) 8b . The missing 1 and 6 representations arise at weight 4. Even though one can form 351 products Y i Y j where some vanishing modular forms easily seen from the Clebsch-Gordan coefficients are not counted, the space of modular forms of weight 4 (and level 7) has dimension 14k − 2 = 54. Therefore, there are 297 constraints between the Y i Y j , which we list in Appendix D. The 54 linearly independent modular forms of weight 4 can be arranged into the following multiplets of Γ 7 : 3b , Y 3c , Y 3a , Y 3b , Y 6a , Y 6b , Y 6c , Y 6d , Y Table 1: Summary of modular forms of level 7 up to weight 8, the subscript r denotes the transformation property under Γ7 modular symmetry. Here Y (2) 8a and Y (2) 8b stand for two weight 2 modular forms transforming in the representation 8 of Γ7. The same convention is adopted for other modular forms.
. (32) 4 Lepton models based on Γ 7 modular symmetry In this section, we shall construct some typical models for neutrino masses and mixing based on the Γ 7 modular symmetry. We will not introduce any flavon field, the flavor symmetry is broken when the complex modulus τ obtains a vacuum expectation value. The Higgs doublets fields H u,d are assumed to be singlets of Γ 7 with vanishing modular weights. The three right-handed (RH) charged leptons E c 1,2,3 transform as singlet 1 under Γ 7 modular group nevertheless they are distinguished by the different modular weights k 1,2,3 . We assign the three generations of left-handed (LH) lepton doublets L and the three right-handed neutrinos N c to two triplets 3 and3 with the weights k L and k N . We shall employ potentially the lowest weight modular forms as much as possible in order to reduce free parameters.

Charged lepton sector
If the left-handed lepton fields L are embedded into the triplet 3, modular forms in the represen-tation3 should be invoked in the charged lepton mass terms. The superpotential for the charged lepton Yukawa coupling reads as Notice that there are two linearly independent weight ten modular forms Y (10) 3a and Y (10) 3b transforming as3 at level N = 7. The phases of the coupling constants α, β and γ 1 can be absorbed into the lepton fields while γ 2 is generally a complex parameter. Modular invariance of the superpotential W e in Eq. (33) requires the modular weights should fulfill the conditions After the value of τ is fixed by certain modulus stabilization mechanism, the charged lepton mass matrix takes the following form where we denote Y 3,2 , Y 3,3 ) T , and the expressions of Y 3,1 , Y 3,2 and Y 3,3 are given in Appendix D. Similar notations are adopted for Y ρ L k 1,2,3 + k L Charged lepton mass matrices Applying the decomposition rules of Γ 7 , we obtain the light neutrino mass matrix as 6a,6b stands for Y (4) 6a,6b and Y (6) 6a,6b for k L = 2 and k L = 3, respectively. If the light neutrino masses are generated by the type I seesaw mechanism, the superpotential for the neutrino masses can be generally written as where Y D and Y N denote the modular form multiplets, and they are required to ensure modular invariance. The explicit forms of Y D and Y N are determined by the weight and representation assignments for L and N c . The Majorana mass term for the heavy neutrinos N c is similar to the Weinberg operator in Eq. (41), and Y N should be modular form multiplets transforming as 6 under Γ 7 . The mass matrix for the Majorana neutrinos N c reads as 6a,6b for k N = 3. Now we analyze the neutrino Yukawa couplings.
In this case, modular invariance requires Y D should be transform in 3 or 6. In the case of k L + k N = 2, and the Dirac neutrino mass matrix take the following form Since this m D is a 3 × 3 anti-symmetric matrix with zero determinant, the light neutrino mass matrix given by seesaw formula is at most of rank 2, such that at least one light neutrino is massless. If k L + k N = 4, the neutrino Yukawa couplings would involve three independent terms, and we have Similar in the case of k L + k N = 6, we have (48) It contains three complex input parameters.
For this type of assignment, L and N c can form an invariant singlet for k L + k N = 0, and consequently m D has a simple structure, The lowest non-vanishing weight k L +k N = 2 gives rise to the following neutrino Yukawa coupling The Dirac neutrino mass matrix is given by which contains one more complex parameter than the above case of k L + k N = 0. A number of free coupling constants are introduced for higher modular weights, the expression of m D would be too lengthy to display. For the assignment L ∼ 3 and N c ∼3, the matrix m D can be obtained The models in which neutrino masses are generated through Weinberg operator and seesaw mechanism are denoted as W1,2 and S1,2,3,4 respectively.
from Eq. (49) by performing transposition. The different possible forms of the neutrino mass matrices for Weinberg operator and type I seesaw mechanism are summarized in table 3. We see that the light neutrino mass matrix of S 2 contains more free parameters than other possible cases. Hence we will not perform numerical analysis for the S 2 model in the following.

Benchmark models
In the following, we consider two scenarios: the modular symmetry only acts on the neutrino sector and the charged lepton matrix is diagonal in the first scenario and the modular symmetry acts on both charged lepton and neutrino sector in the second scenario.

Modular symmetry on neutrino sector
The possibility that the dynamics of flavor in the charged lepton and neutrino sectors are different can not be excluded [18]. For instance, certain flavon may be involved in the charged lepton sector and neutrino sector is dictated by modular symmetry [18]. For simplicity, we assume that charged lepton sector is diagonal. In this case, there are six different models which are shown in table 3.
We find that the model S 1 can not accommodate the experimental data, the model S 2 contains more input parameters than other models, consequently we will not discuss it. Then the rest four models in table 3 only depend on the following four inputs τ, τ, |g 2 /g 1 |, arg (g 2 /g 1 ) , and an overall parameter which can be fixed by the mass squared difference ∆m 2 21 = m 2 2 − m 2 1 . The five dimensionless observable quantities: only depend on the four input parameters in Eq. (52), where ∆m 2 3 = m 2 3 − m 2 1 > 0 for normal ordering (NO) and ∆m 2 3 = m 2 3 − m 2 2 < 0 for inverted ordering (IO) [49]. In order to quantitatively assess how well a model can describe the experimental data on the five dimensionless observable quantities in Eq. (53). We define a χ 2 function to estimate the goodness-of-fit of a set of chosen values of the input parameters, where O i denote the global best fit values of the five observable quantities in Eq. (53), and σ i refer to the 1σ deviations of the corresponding quantities, and P i are the theoretical predictions for the five physical observable quantities for the input parameters taking certain values. Here the contribution of the Dirac phase δ CP is also included in the χ 2 function. For each value of the input parameters, one can obtain the predicted values P i and the corresponding χ 2 , then one can find out the lowest χ 2 . After performing a detailed numerical analysis for the three mixing angles, Dirac CP phase and ∆m 2 21 /∆m 2 3 , we find that only models W 1 with k L = 3 and S 4 with k N = 3 for NO case and models W 2 with k L = 3 and S 3 with k N = 3 for IO case can give results in agreement with the experimental data. As an example, we only show the results of NO case. For model W 1 with k L = 3, we find the minimum of χ 2 is χ 2 min = 1.937, and the best fit values of the free parameters are τ = 0.448, τ = 0.915, |g 2 /g 1 | = 0.129, arg (g 2 /g 1 ) = 0.566π, g 2 1 v 2 u /Λ = 56.981 meV .
where m β is the effective mass probed by direct kinematic search in beta decay and m ββ refers to the effective Majorana mass in neutrinoless double beta decay. The latest result from KATRIN is m β < 1.1 eV at 90% CL [50].  Figure 1: The correlation between δCP and sin 2 θ23, the left and right panels are for the models S4 with kN = 3 and W1 with kL = 3 respectively in the charged lepton diagonal basis. The gray bands denote the experimentally preferred 1σ and 3σ ranges of δCP and sin 2 θ23 adapted from [49].
EXO-200 [51], and they would potentially be tested in next generation experiments. As regards the experimental bound on neutrino mass sum, the result sensitively depends on the cosmological model and the experimental data considered. Combining the Planck TT, TE, EE, lowE polarization spectra, baryon acoustic oscillation (BAO) data with the CMB lensing reconstruction power spectrum, the Planck collaboration gives i m i < 120 meV at 95% confidence level [52]. However, if only the BAO data and the CMB lensing reconstruction power spectrum are taken into account in the data analysis, this bound becomes i m i < 600 meV [52]. For the above two models, we find the neutrino mass sum i m i is 207.649 meV and 376.853 meV respectively which are consistent with the Planck's looser constraint i m i < 600 meV.
We perform a comprehensive numerical scan over the free parameters of the above two models. We find that the three mixing angles can take any values in their 3σ ranges. The two Majorana CP phases are restricted to the ranges α 21 /π ∈

Modular symmetry on both neutrino and charged lepton sector sectors
Models mass matrices ρ L , ρ N c modular weights Combining the different possible constructions in the charged lepton and neutrino sectors, we can easily get all possible models based on Γ 7 modular symmetry. Focusing on the cases with lower  [49]. The horizontal grey band denote the current experimental bound m ββ < (120 − 250) meV from KamLAND-Zen [51]. The vertical grey exclusion band denotes the bound on the lightest neutrino mass extracted from i mi < 120 meV by the Planck collaboration [52].
weight modular forms and free parameters as few as possible, we find four different types of models named as M 1,2,3,4 . The assignments of the weights and representations for the leptonic fields are listed in table 4. From the superpotential of C 1 and C 2 , we see that the phases of the parameters α, β and γ 1 can be absorbed into the right-handed lepton fields. Then the coupling constants α, β and γ 1 in the charged lepton mass matrix can be taken to be real and positive and γ 2 is a complex parameters for both C 1 and C 2 . As a consequence, the charged lepton sector only contain four real dimensionless parameters β/α, γ 1 /α, |γ 2 /α|, arg (γ 2 /α) and an overall scale αv d . In the neutrino sector, it is easy to check that light neutrino mass matrices for all the four cases W 1 , W 2 , S 1 and S 2 depend on two positive dimensionless parameters |g 2 /g 1 |, arg (g 2 /g 1 ) and an overall neutrino mass . We perform a numerical analysis for each model, the complex modulus τ is restricted to lie in the fundamental domain {τ | τ > 0, | τ | ≤ 1 2 , |τ | ≥ 1} of the modular group. The other dimensionless parameters randomly vary in the following regions arg (γ 2 /α), arg (g 2 /g 1 ) ∈ [0, 2π) , β/α, γ/α, γ 1 /α, |γ 2 /α|, |g 2 /g 1 | ∈ [0, 10 4 ] .
The overall parameters of charged lepton and neutrino mass matrices are fixed by the electron mass and the solar neutrino mass squared difference ∆m 2 21 . Then one can obtained all the predictions for six lepton masses and three lepton mixing angles as well as three CP violating phases. We find that good agreement with experimental data can be achieved for certain values of input parameters for all these four models in both NO and IO cases. Since NO is slightly preferred by present data, we only present the numerical results of the four models for NO. The predictions for lepton mixing parameters and neutrino masses are listed in table 5. The charged lepton mass matrix depends on four parameters α, β, γ 1 and γ 2 , the measured values of the three charged lepton masses can be reproduced exactly, hence we do not show the results of charged lepton masses in table 5. We find that the global best fit values of neutrino mixing angles and mass squared differences ∆m 2 21 , ∆m 2 31 from NuFIT v4.1 can be obtained, as shown in table 5. Moreover, the most stringent neutrino mass bound i m i < 120 meV [52] is saturated in models M 2 and M 3 , and the less stringent bound Best fit values for NO

Conclusion
We have considered the finite modular group Γ 7 P SL(2, Z 7 ) in the framework of the modular invariance approach to lepton flavour. Γ 7 is a quotient group of the infinite modular group Γ which is achieved by imposing the generator condition T 7 = 1. An additional condition (ST 3 ) 4 = 1 should also be satisfied, which is essential to make Γ 7 finite.
One crucial ingredient of modular-invariant theories is the introduction of modular forms, which involves the modulus field τ . Given level 7 and an even weight k, there are 14k − 2 linearly independent modular forms, which can all be decomposed into irreducible representations (irreps) of Γ 7 . At weight k = 2, we constructed all 26 modular forms with the help of the SageMath algebra system [46]. They are decomposed into a triplet 3, a septet 7 and two octets 8 of Γ 7 . A full list of linearly independent modular forms of weight up to 8 is provided in Appendix D. We have also considered two alternative ways to derive modular forms of level 7, the Dedekind eta function method proposed in [1] and the theta function method proposed in [20], discussed in Appendices B and C, respectively. Results from these methods are consistent with the former one, but incomplete: the Dedekind eta function method gives only 7 modular forms of weight 2, forming the septet of Γ 7 ; and the theta function method gives 23 modular forms of weight 2, decomposed as 7 + 8 + 8.
We also considered flavon-free lepton flavour models, constructed by assigning couplings and matter superfields transforming in irreps of Γ 7 . In the considered models, charged leptons gain masses via renormalisable Yukawa couplings, neutrinos gain masses via either the Weinberg operator or the type-I seesaw mechanism. Flavour textures arise after the modulus field τ gains a vacuum expectation value. Lists of charged lepton (C 1 , C 2 ) and neutrino mass matrices (W 1 , W 2 and S 1 , ..., S 4 ) involving τ are given in Tables 2 and 3, respectively.
In the numerical studies, we have considered two scenarios: the modular symmetry acting on only the neutrino sector, or both charged lepton and neutrino sectors. We have performed a χ 2 analysis with experimental data of neutrino oscillation parameters in the 1σ range taken into account. In the first scenario, we found that, in the normal neutrino mass ordering, only models W 1 with lepton doublet weight k L = 3 and S 4 with right-handed neutrino weight k N = 3, give results that agree with the experimental data. Similarly, in the inverted mass ordering, only models W 2 with k L = 3 and S 3 with k N = 3 are allowed by data. The effective neutrino masses m β and m ββ , which control the beta decay and neutrinoless double beta decay rates, respectively, are predicted to be compatible with current data. However, the prediction of the sum of neutrino masses is disfavoured by the current cosmological constraint i m i < 120 meV [52] although the less stringent bound i m i < 600 meV [52] is satisfied. In the second scenario, eight benchmark models (M 1 , ..., M 4 with two sets of modular weights), listed in Table 4, have been studied. While all of them are compatible with oscillation data, the models M 2 and M 3 are also favoured by the cosmological constraint on the sum of neutrino masses i m i < 120 meV, while the prediction of m ββ is within the sensitivity of the next generation of neutrinoless double beta decay experiments.
A Group Theory of Γ 7 ∼ = P SL(2, Z 7 ) The group Γ 7 ∼ = P SL(2, Z 7 ) is a non-Abelian finite subgroup of SU (3) of order 168. Γ 7 group can be generated by two generators S and T which satisfy the multiplication rules: The 168 elements of Γ 7 group are divided into 6 conjugacy classes:

Conjugacy Classes
where nC k denotes a class with n elements which is of order k. The character table of Γ 7 group is given in Table 6. Following the convention of Ref. [45], we find that Γ 7 group has ninety-two abelian subgroups in total: twenty-one Z 2 subgroups, twenty-eight Z 3 subgroups, fourteen K 4 subgroups, twenty-one Z 4 subgroups and eight Z 7 subgroups. In terms of the generators S and T , these abelian subgroups are given as follows: All the above twenty-one Z 2 subgroups are conjugate to each other.
All the fourteen K 4 subgroups are conjugate as well.

D Higher weight modular forms and constraints
Through the tensor products of the modular forms Y 8a and Y (2) 8b , one can find, at weight 4, the following modular multiplets: (D.32) Notice that not all of the above modular multiplets are linearly independent. From the q-expansions of Y i (τ ) given in Eq. (23), we find the following 297 constraints between the different weight 4 modular multiplets Y 1 Y