Landscape of Modular Symmetric Flavor Models

We study the moduli stabilization from the viewpoint of modular flavor symmetries. We systematically analyze stabilized moduli values in possible configurations of flux compactifications, investigating probabilities of moduli values and showing which moduli values are favorable from our moduli stabilization. Then, we examine their implications on modular symmetric flavor models. It is found that distributions of complex structure modulus $\tau$ determining the flavor structure are clustered at a fixed point with the residual $\mathbb{Z}_3$ symmetry in the $SL(2,\mathbb{Z})$ fundamental region. Also, they are clustered at other specific points such as intersecting points between $|\tau|^2=k/2$ and ${\rm Re}\,\tau=0,\pm 1/4, \pm1/2$, although their probabilities are less than the $\mathbb{Z}_3$ fixed point. In general, CP-breaking vacua in the complex structure modulus are statistically disfavored in the string landscape. Among CP-breaking vacua, the values ${\rm Re}\,\tau=\pm 1/4$ are most favorable in particular when the axio-dilaton $S$ is stabilized at ${\rm Re}\,S=\pm 1/4$. That shows a strong correlation between CP phases originated from string moduli.


Introduction
In particle physics, it is one of important issues to understand the origin of the flavor structure in the quark and lepton sectors, that is, the hierarchies of quark and lepton masses, their mixing angles and CP phases. Indeed, various scenarios have been proposed and studied. Among them, the approach by non-Abelian discrete flavor symmetries is one of interesting approaches, and many studies have been carried out by assuming various non-Abelian discrete flavor symmetries such as A N , S N , ∆(6N 2 ) [1][2][3][4][5]. For example, the A 4 flavor symmetry is one of most extensively studied models. It is a symmetry of tetrahedron. Thus, geometrical symmetries can be origins of these non-Abelian discrete flavor symmetries in higher-dimensional theories such as superstring theory [6,7]. Symmetries are important tools to connect physics between low-energy and high-energy scales.
The torus and orbifold compactifications also have another geometrical symmetry called the modular symmetry, which corresponds to a basis change of cycles on the torus. The twodimensional (2D) tours has the modular symmetry SL(2, Z), and higher-dimensional tori have larger modular symmetries, although the factorizable torus, e.g. T 2 × T 2 × T 2 , has the product symmetry of SL(2, Z). Furthermore, the symplectic modular symmetry appears in resolutions of toroidal orbifolds such as Calabi-Yau compactifications [8,9].
The modular group also transforms non-trivially zero-modes corresponding to quarks and leptons in four-dimensional (4D) effective field theory. That is the flavor symmetry. Yukawa couplings as well as higher-order couplings are functions of moduli. Thus, these couplings also transform non-trivially under the modular flavor symmetry. So far, the modular flavor symmetry are discussed in the context of Narain lattice in the heterotic string theory [10][11][12][13] and magnetized D-brane models in Type IIB string theory [14][15][16][17] from the top-down approach. Furthermore, unification of the modular symmetry and the 4D CP is also developed in the heterotic string theory and in Type IIB string theory, where the 4D CP is identified with an outer automorphism of the SL(2, Z) modular group on toroidal background [18,19] and Sp(2n, Z) modular group on Calabi-Yau threefolds [20]. That can be an origin of generalized CP symmetry [19,[21][22][23].
Interestingly, the modular group has finite subgroups such as S 3 , A 4 , S 4 , A 5 , ∆(96), ∆(384) [24]. These discrete groups have been used to construct phenomenologically interesting flavor models in the bottom-up approach [1][2][3][4][5]. Motivated by the above aspects, a modular symmetric flavor model was proposed in Ref. [25] as a phenomenological model building in the bottomup approach, and recently modular flavor symmetries are studied extensively. (See for early works [26,27].) Also such studies are extended to covering groups of Γ N [28]. In contrast with well-studied models with discrete flavor symmetries, these modular symmetric flavor models are controlled by the modular symmetry, and Yukawa couplings as well as higher-order couplings are written by modular forms, which are functions of the modulus τ . A vacuum expectation value of the modulus τ induces the spontaneously symmetry breaking of the modular symmetry, and consequently determines the mass hierarchy and mixing angles of quarks and leptons. Also, the axionic part of τ determines CP phases. Hence, the value of τ is the keypoint to realize realistic fermion mass matrices in these modular flavor models. In the bottom-up approach, the whole fundamental region of SL(2, Z) is often searched in order to compare the predictions of modular flavor models with the observational data. There is no guideline to pick up particular points in the moduli space. Furthermore, it is unclear whether the phenomenologically favored moduli values discussed in the bottom-up approach can be realized by some mechanism. In this respect, the ultra-violet theory would provide us with a hint on this problem, that is, the moduli stabilization problem.
In both bottom-up and top-down approaches, the moduli stabilization is important and necessary to determine the flavor structure. Unless these moduli fields are stabilized at a high scale larger enough than the current observational bound, the predictability of the modular flavor models is lost. For example, in Refs. [29,30] the possibilities of favorable moduli stabilization were studied by assuming modular symmetric non-perturbative effects.
One of powerful approaches to realize the moduli stabilization is flux compactifications, taking into account background field values of higher-dimensional tensor fields. For instance, in Type IIB flux compactifications, discrete modular symmetry is classified on toroidal background [31], and spontaneous CP-violation is discussed in toroidal [32] and Calabi-Yau back-grounds [20]. In this paper, we study how to fix the modulus value τ in the light of the moduli stabilization due to flux compactifications.
We deal with possible configurations of three-form fluxes in Type IIB string theory. Then, we analyze systematically stabilized values of moduli, investigating their probabilities in the SL(2, Z) fundamental region and showing which values are favorable from the moduli stabilization. That is the so-called string landscape. Then, we examine implications of our moduli stabilization from the viewpoint of modular flavor models. For illustrative purposes, we consider classes of modular A 4 models extensively studied in Ref. [33]. Then, we compare the predictions of these phenomenological models with our results in the string landscape. So far, the structure of Type IIB string landscape has been statistically studied in Refs. [34][35][36] and in particular, the enhanced symmetries appear in particular points in moduli spaces of the complex structure moduli and the axio-dilaton [37]. Such distributions of the moduli fields are useful to understand not only the nature of string landscape itself, but also phenomenological aspects of the 4D effective action. Indeed, the complex structure moduli determine the flavor structure of quarks and leptons on magnetized D-brane setup. It is expected that the probability distributions of complex structure moduli shed new light on the phenomenological approach. Our findings about the distribution of the complex structure modulus τ and the axio-dilaton S are summarized as follows: • Distribution of the complex structure modulus τ determining the flavor structure of quarks and leptons is mostly clustered at smaller vacuum expectation values rather than the scattered distribution. Remarkably, the Z 3 fixed point on the fundamental domain of the SL(2, Z) moduli space is statistically favored by the weak string coupling and the tadpole cancellation condition. Also, distributions are clustered at other specific points such as intersecting points between |τ | 2 = k/2 and Re τ = 0, ±1/4, ±1/2, although their probabilities are less than the Z 3 fixed point.
• Distribution of the complex structure modulus τ depends on values of the axio-dilaton S, although the Z 3 fixed point is universally favored in the landscape. For instance, the Z 2 fixed point in the moduli space of the complex structure modulus τ is favored on the same Z 2 fixed point in that of the axio-dilaton S.
• The CP-violating phase is determined by the vacuum expectation value of the complex structure modulus τ . We find that an axionic field of the complex structure modulus Re τ is mostly stabilized at the CP-conserving vacua Re τ = 0, ±1/2 or the CP-breaking vacua Re τ = ±1/4, although the CP-breaking vacua is statistically disfavored in the landscape. However, the CP-breaking vacua Re τ = ±1/4 are statically favored only if the axion Re S has certain values. Thus, there is a strong correlation among CP phases in Yukawa couplings and the CP phase due to the axion Re S.
• When we apply the distribution of the complex structure modulus to the modular A 4 models in Ref. [33], one can predict the theoretically favorable vacuum expectation values of the moduli and its probability in the string landscape.
This paper is organized as follows. First, we review the Type IIB supergravity action on T 6 /(Z 2 × Z 2 ) orientifold. The modular symmetries associated with the complex structure of tori and the axio-dilaton, and the stabilization of the moduli fields are reviewed in Section 2. We next analyze the distributions of complex structure moduli in Section 3, in which the phenomenologically applicable results are summarized. In Section 4, we compare the phenomenological predictions of modular A 4 models with the distributions of the complex structure moduli. Finally, we conclude this paper in Section 5. The detail of the numerical search is summarized in Appendix A.

Theoretical setup
In this section, we briefly review the phenomenologically attractive modular symmetry appearing in the string-derived effective action, with an emphasis on Type IIB supergravity on the factorizable T 6 torus and its orientifold. The dynamics of the moduli fields is discussed in the context of flux compactifications as reviewed in Section 2.2 on a simple T 6 /(Z 2 × Z 2 ) orientifold background. The reader who are interested in the distributions of the moduli fields may skip to Section 3.

Modular symmetry
Let us consider the factorizable 6D torus T 6 = (T 2 ) 1 ×(T 2 ) 2 ×(T 2 ) 3 , each of which is defined by a complex plane divided by the lattice Λ i . The geometrical symmetry of each torus is given by the SL(2, Z) i modular symmetry. When we represent the basis of the lattice Λ i by (e x i , e y i ), the basis of the lattice transforms as where are the element of SL(2, Z) i satisfying p i t i − q i s i = 1. Such a passive transformation of the basis induces the modular transformation of the complex structure τ i as well as the coordinate z i with reference to e x i as [17] where u i denote the coordinates of the complex plane. The generators of each SL(2, Z) i are given by satisfying the following algebraic relations These S-and T -transformations bring τ i to the fundamental domain of the SL(2, Z) i moduli space: To see the modular symmetry in the effective action of complex structure moduli, let us introduce the three-form basis in H 3 (T 6 , Z), 1 which satisfy the orientation with I, J = 0, 1, 2, 3. Then, the 4D kinetic terms of the complex structure moduli are obtained by the dimensional reduction of 10D Einstein-Hilbert action. It results in the following Kähler potential: where we define the holomorphic three-form on T 6 as The modular symmetries of the complex structure moduli are checked to see the Kählerinvariant quantity, which consists of the Kähler potential and the superpotential characterizing the 4D effective potential. Since the Kähler potential of the complex structure moduli transforms under the modular symmetry as the effective potential is modular invariant up to the Kähler transformation, only if the superpotential has the modular weight 1 under each SL(2, Z) i modular symmetry, that is, So far, there exists the 4D N = 8 supersymmetry originating from the 10D N = 2 supersymmetry in the case of Type II string theory. To reduce N = 8 supersymmetry to N = 1 supersymmetry, we impose the orientifold and orbifold projections on T 6 background. Especially, we focus on M = T 6 /(Z 2 × Z 2 ) orientifold throughout this paper. The Z 2 × Z 2 orbifold projections are chosen as and the orientifold projection is specified by We denote the world-sheet parity projection by Ω p and the left-moving fermion number operator by F L , respectively. It is known that under the latter two operators Ω p (−1) F L , the metric and the axio-dilaton are even and other higher-dimensional form fields such Kalb-Ramond field B and Ramond-Ramond (RR) fields C 2p behave as Note that all three-form basis in Eq. The kinetic terms of the complex structure moduli are also described in Eq. (9) in the same way as the factorizable T 6 background.

Flux compactifications on
The moduli stabilization is inevitable to discuss the flavor structure of quarks and leptons. In particular, the flux compactification is useful to determine the vacuum expectation values of the moduli fields in a controlled way. In this paper, we focus on T 6 /(Z 2 × Z 2 ) orientifold as explained in the previous section.
On the T 6 /(Z 2 × Z 2 ) background, the effective Kähler potential of the moduli fields is described by where we include the axio-dilaton S and the volume of the torus V in the Einstein-frame measured in units of the string length l s = 2π √ α .
Here and in what follows, we adopt the reduced Planck mass unit M Pl = 1. The effective action has three-types of SL(2, Z) symmetry, namely SL(2, Z) S for the axio-dilaton, SL(2, Z) T for the volume moduli and SL(2, Z) τ for the complex structure moduli. 2 In Type IIB string setup with D3/D7-branes, the complex structure moduli is of particular interest for the flavor structure of quarks and leptons localized on D-branes. These complex structure moduli and the axio-dilaton can be stabilized in the context of flux compactifications.
In the context of Type IIB string theory, there exists the kinetic term of three-form G 3 consisting of Ramond-Ramond (RR) F 3 and Neveu-Schwarz (NS) three-forms H 3 , namely G 3 = F 3 − SH 3 . After the dimensional reduction of this kinetic term on T 6 /(Z 2 × Z 2 ) background and taking into account background fluxes of these three-forms, one can obtain the modulidependent superpotential in the four-dimensional effective action [39] The RR and NS three-forms are expanded on the basis of Eq. (8); where we denote the integral flux quanta by {a 0,1,2,3 , b 0,1,2,3 , c 0,1,2,3 , d 0,1,2,3 }. On T 6 /(Z 2 × Z 2 ) geometry, these flux quanta are restricted to be in multiples of 4 or 8 depending on the existence of discrete torsion [40][41][42]. Throughout this paper, we analyze the latter case, since one can focus on the dynamics of 3 untwisted complex structure moduli, denoted by τ 1,2,3 in the effective action. Then, the explicit form of the superpotential is given by which generates the potential of the axio-dilaton and complex structure moduli. Indeed, by combining the superpotential (20) and the Kähler potential (17), the 4D scalar potential is constructed in the framework of 4D N = 1 supergravity. Since the superpotential is independent of the Kähler moduli, the 4D scalar potential can be simplified as with I, J = S, τ 1 , τ 2 , τ 3 . Here, we denote the inverse of the Kähler metric K IJ by K IJ as computed from the Kähler potential K IJ = ∂ I ∂J K, and the covariant derivative of the superpotential is given by Note that the negative term −3|W | 2 in the scalar potential is cancelled by the no-scale structure of the Kähler moduli. Before going into the detail to analyze the structure of the scalar potential, we discuss the symmetry possessed in the flux-induced effective action. To simplify our analysis, we concentrate on the isotropic torus, where the complex structure moduli are constrained on the following locus: and correspondingly, flux quanta are redefined as meaning that we have totally 8 independent flux quanta; {a 0 , a, b, b 0 } for RR flux and {c 0 , c, d, d 0 } for NS flux. This system has three types of symmetries as shown below. First, as mentioned in Section 2.1, the effective action is invariant under the SL(2, has the modular weight 3, which is different from 1 owing to the identification of three complex structure moduli. Under S-and T -transformations of SL(2, Z) τ modular group, the effective action is indeed modular invariant when the flux quanta transform as [46], Second, the axio-dilaton enjoys the SL(2, Z) S modular symmetry. Under the SL(2, Z) S transformation of the axio-dilaton with being the element of SL(2, Z) S satisfying pt − qs = 1, the effective action is modular invariant under the following transformation of the RR and NS fluxes, Specifically, the flux quanta transform under the S-and T -transformations of SL(2, Z) S as Third, the effective action is invariant under flipping the overall sign of flux quanta, namely which changes the sign of the superpotential, W → −W , but the effective action is invariant, taking into account the Kähler transformation. This symmetry can be used to reduce SL(2, Z) τ,S to P SL(2, Z) τ,S . Then, these three symmetries are employed to count the physically-distinct vacua as analyzed in Section 3. Finally, we comment on constraints for the flux quanta. These flux quanta induce the D3-brane charge through the Bianchi identity of the RR field, which should be cancelled by the contributions of D3/D7-branes and O3/O7-planes 3 . It was known that N flux are bounded as 0 ≤ N flux ≤ O(10) for some explicit setups on the toroidal T 6 /Z M and T 6 /(Z M × Z N ) orientifolds [43]. The positivity of N flux is a consequence of the supersymmetry. When we uplift the Type IIB orientifolds to its strong coupling regime, namely F-theory, O3-plane contributions are encoded in the Euler number of Calabi-Yau fourfolds, whose largest value is known as 1820448 [44,45]. In this case, the flux-induced D3-brane charge is bounded as Similar to the analysis in Ref. [46], we adopt this bound (35) in the following analysis, where we examine the distributions of moduli fields by changing the maximum value of the flux-induced D3-brane charge, N max flux . Since the flux quanta are in multiple of 8 on T 6 /(Z 2 × Z 2 ) geometry, N flux should be in multiple of 192, namely 4

Supersymmetric minima
In this section, we consider the supersymmetric minimum solutions of the moduli fields on the isotropic torus.
The supersymmetric minimum conditions of the moduli fields are given by where the third condition is due to the supersymmetric condition for Kähler moduli 5 .
First of all, we separate the superpotential into RR and NS parts as The minimum conditions ∂ S W = 0 and W = 0 lead Each of the equations, W RR = 0 and W NS = 0, is the cubic equation with respect to τ and the coefficients are real. The stabilized value Im τ must be positive, i.e. Im τ > 0; otherwise, the volume of T 6 is vanishing. To find a solution τ with its non-zero imaginary part, the cubic equation must have two common complex solutions which are conjugate with each other and a real one which need not be common. Taking into account them, W RR and W NS should be factorizable as 6 for a quadratic (integer-coefficient) polynomial P (τ ). It indicates that the complex τ is obtained from P (τ ) = 0. Next, let us consider the condition ∂ τ W = 0 whose explicit form is given by For a complex solution τ given by P (τ ) = 0, it is required to satisfy ∂ τ P (τ ) = 0 or Indeed, ∂ τ P (τ ) = 0 leads a real τ so that the vacuum expectation value of S is given as above. 7 To obtain an explicit form of τ , we set for m 2 − 4ln < 0, and then we arrive at the following vacuum expectation value of τ : As we reviewed in Section 2.2, there are two symmetries for S, τ , namely SL(2, Z) S , SL(2, Z) τ which leave the effective action invariant. Hence, there exists an equivalence relation between two given flux vacua, which can be found by checking the invariance of supersymmetric conditions under the modular transformations. For the detail about the equivalence, see, Ref. [46]. However, to be self-contained, we briefly summarize such an equivalence in Appendix A. We have to be careful not to make double-counting there.

Distributions of the moduli fields
In this section, we systematically analyze the distributions of complex structure modulus τ as well as the axio-dilaton S by the moduli stabilization through flux compactifications as explained in the previous section. We set the isotropic torus τ = τ 1 = τ 2 = τ 3 . On the basis of the setup in the previous section, we show our results on distributions of the complex structure modulus τ and the axio-dilaton S. In Section 3.1, we first investigate the distributions of stable vacua on the plane of complex structure modulus τ , taking into account the upper bound for the flux-induced D3-brane charge N max flux in Eq. (35), which fixes the maximal size of fluxes, and the lower bound for the string coupling g s = Im S −1 . These results indicate the correlation between the distributions of the complex structure modulus and the axio-dilaton. We next deal with the distributions of stable vacua on the plane of the axio-dilaton in Section 3.2, especially we pick up statistically favored vacuum expectation values of the axio-dilaton and discuss its consequence for the distribution of the complex structure modulus. Section 3.3 is devoted to analyze the CP-breaking vacua.

Distribution of complex structure modulus without fixing the axio-dilaton
We analyze the structure of the supersymmetric vacua as analytically demonstrated in Section 2.3. We search all the possible flux configurations (a 0 , a, b, b 0 , c 0 , c, d, d 0 ), which satisfy the tadpole cancellation condition (35) for each 8 N max flux = 10, 100, 250, 500, 1000.
By use of fluxes, we examine stabilized moduli values, satisfying the minimum conditions (37), i.e., Eqs. (44) and P (τ ) = 0. Since the effective action has three types of symmetries as mentioned in Section 2.2, it is required not to double count the number of stable vacua. As a result, we find the number of physically-distinct stable vacua as summarized in Table 1.
In this section, we first focus on the distribution of the complex structure modulus τ without imposing any bound for the axio-dilaton S. The distribution of the complex structure modulus τ in the fundamental region of the SL(2, Z) τ moduli space (6) is drawn in Figure 1 for the case with N max flux = 10 in the left panel and N max flux = 1000 in the right panel, respectively. We find that the stable vacua are clustered in the dark color region in both panels and especially, most of the stable vacua are clustered at smaller vacuum expectation values of Im τ rather than the scattered distribution. That is a consequence of the condition for the flux quanta through the tadpole cancellation condition (35).
Furthermore, the Z 3 fixed point of the SL(2, Z) modular group is favored in the landscape. We recall that the fixed points τ fix in the fundamental region of the SL(2, Z) moduli space are defined as where γ 0 is the element of SL(2, Z) transformation called as a stabilizer of τ fix . It is known that there exist three fixed points in the fundamental region such as , where the first two points and the latter one correspond to the Z 2 and Z 3 fixed points, respectively. (Re τ, Im τ )  Table 2: Probabilities of the stable vacua as functions of (Re τ, Im τ ) in the descending order of the probability for N max flux = 10.  Table 3: Probabilities of the stable vacua as functions of (Re τ, Im τ ) in the descending order of the probability for N max flux = 1000.
Indeed, we calculate the probability of stable vacua as shown in Tables 2 and 3. 10 These 10 Here and in what follows, we add ReX = 1/2 (X = τ, S) lines and the right side of unit circles to make results indicate that the stable vacua are clustered at the Z 3 fixed point, when N max flux is small. Note that the smallness of N max flux is more theoretically controllable, because we require the extension of Type IIB string setup to the F-theory extension in N max flux ∼ O(10 3 ) regime. In addition to the Z 3 fixed point, there exist other favorable modulus values such as the Z 2 fixed point (Re τ, Im τ ) = (0, 1) and (Re τ, Im τ ) = (0, √ k) with k = 2, 3, 6. All of these are on boundary of the GL(2, Z) fundamental region. Furthermore, all of them satisfy where k is positive integer. That is, all of the points in Tables 2 and 3 4 ) are intersecting points between |τ | 2 = k/2 and Re τ = 0, −1/2. 11 In addition to the points in Tables 2 and 3, the following points: have O(0.1-1)% of probability. These points are also on the circles with the radii k/2 including (Re τ, Im τ ) = (±1/4, √ 15/4). These points are important from the viewpoint of CP violation, as will be studied in Section 3.3.
In general, from the analytic expression of τ in Eqs. (46) and (47), one can see that i.e. all the vacuum expectation values of τ in the upper half-plane are on the semicircles and each radius of them is given by Eq. (54). Note that this does not mean that each semicircle becomes a solution to Eq. (37) in its entirety. Indeed, the solutions are discrete in the first place. If we assume that the vacua are distributed in those circles, it is needed to fold back into the region by using T -transformations since these circles protrude from the fundamental region. This phenomenon can be seen from the pattern of Fig. 1, and it implies that the overlap points including ( 1 4 , ) and (0, √ k) listed above are statistically favored. Our results are in agreement with the statistical approach developed in Refs. [34][35][36], in which the discrete flux quanta are taken to be continuous one. As demonstrated on the toroidal background [37], the number of W = 0 supersymmetric vacua with fixed complex structure modulus behaves as the figures symmetric. However, to define probabilities correctly, we would not do the same thing on the tables including the probabilities. In addition, we divided 2D moduli space into domains whose size are (0.01 × 0.01) to draw the figures, but this configuration still captures fine structures such as voids. 11 In Ref. [32], it was shown that the rescaled modulus φ satisfies |φ| 2 = 1 for CP-symmetric superpotential.
where t(x) = x for x ≡ 0 (mod 3) and otherwise t(x) = 3x, and {l, m, n} determine the vacuum expectation value of the complex structure modulus τ by Eq. (46). Note that the analysis in Ref. [37] adopts the parametrization gcd(l, m, n) = 1 in the equation P (τ ) = 0 determining the vacuum expectation value of the modulus τ , but the ratio between the number of vacua would be the same with our results when the continuous approximation of the flux quanta is valid, that is the N max flux 1 regime. Hence, they would capture the phenomena of our results. The approximated probabilities of Z 3 , Z 2 and other vacua shown in Tables 2 and 3 are calculated as in Table 4, in which we take the values of (l, m, n) appearing in P (τ ) = 0 such as (l, m, n) = (1, 0, 1) for Z 2 fixed point and (l, m, n) = (1, −1, 1) for Z 3 fixed point, respectively. The number of total vacua was approximately calculated as 1.59×10 −2 (N max flux ) 2 [37]. It turns out that Z 3 fixed point is statistically favored in the string landscape. Furthermore, the denominator of Eq. (55) determines the magnitude of Im τ , which means that the smaller value of Im τ is also favored in the string landscape. These phenomena are consistent with our results as seen in the right panel of Figure 1.  Table 4: Approximated probabilities of the stable vacua as functions of (Re τ, Im τ ) shown in Table 3, predicted by the statistical approach.
It is remarkable that the Z 3 fixed point is favorable from the moduli stabilization. At τ = ω, the T 2 /Z 2 orbifold has a high symmetry. The four fixed points on T 2 /Z 2 are equally spaced like a tetrahedron. The same stabilized value τ = ω is realized by another scenario of the moduli stabilization [47].
From the phenomenological point of view, these results support the studies of modular flavor models around specific modulus points. Indeed, phenomenological aspects of modular flavor models are studied around the fixed points of the SL(2, Z) modular group [27,[48][49][50]. The Z 3 and Z 2 residual symmetries remain on these fixed points. Such symmetries control mass matrices of quarks and leptons. Then, phenomenologically interesting results are obtained. It is important that such fixed points are statically favored in the string landscape. In addition to these fixed points, results in Tables 2 and 3 suggest us to study other non-trivial moduli values with definite probabilities. It is interesting to explore the modular flavor models around these modulus values in addition to the fixed points.
So far, we have not imposed any bound for the axio-dilaton S. Let us next impose the bound for the axio-dilaton Im S, focusing on the weak coupling regime. From the viewpoint of the string theory, the axio-dilaton determines the size of the string coupling, g s = (Im S) −1 . Hence, the weak coupling regime is theoretically controllable against perturbative and non-perturbative corrections with respect to the string coupling. We calculate the probability of the stable vacua as functions of the complex structure modulus τ and the axio-dilaton S as summarized in Tables 5 and 6 for N max flux = 10 and N max flux = 1000, respectively. Here we list probabilities in the descending order of highest probabilities of τ . It turns out that the distribution of the stable vacua for N max flux = 1000 are uniformly distributed with respect to the axio-dilaton, but for the weak coupling as well as the small flux regime: Im S > 1, N max flux = 10, all the stable vacua are clustered at the Z 3 fixed point. In this respect, the theoretically controllable weak coupling region also prefers the Z 3 fixed point in the fundamental domain of SL(2, Z) moduli space.
These results motivate us to study the contribution of the axio-dilaton to the distribution of the complex structure modulus. That will be discussed in the next section.
(Re τ , Im τ ) Im S All  Table 5: Probabilities (%) of the stable vacua as functions of (Re τ, Im τ ) and Im S in the descending order of the probability with respect to (Re τ, Im τ ). Here, we list the probability for N max flux = 10 and the column in "All" corresponds to the probability of Table 2 .
(Re τ , Im τ ) Im S All  Table 6: Probabilities (%) of the stable vacua as functions of (Re τ, Im τ ) and Im S in the descending order of the probability with respect to (Re τ, Im τ ). Here, we list the probability for N max flux = 1000 and the column in "All" corresponds to the probability of Table 3 .

Distribution of complex structure modulus with the fixed axiodilaton
In this section, we further study the distribution of the complex structure modulus τ , taking into account the distribution of the axio-dilaton S. As shown in Eq. (44), the vacuum expectation value of the axio-dilaton is correlated with that of the complex structure modulus τ .
Let us work with the distribution of the axio-dilaton without imposing any bound for the complex structure modulus. We project out the distribution of the stable vacua on the plane of the axio-dilaton S in the fundamental region of the SL(2, Z) S moduli space (6), as drawn in Figure 2, in which the stable vacua are clustered in the dark color region. There is a difference between the numbers of stable vacua for N max flux = 10 and N max flux = 1000 as shown in Table 1, but in both cases, the distributions of the stable vacua are not limited on the fixed points in the fundamental domain of SL(2, Z) S moduli space. Indeed, the probabilities of stable vacua in its descending order are summarized in Tables 7 and 8, from which the axiodilaton is distributed along the strong coupling regime Im S ∼ 1 and there are no particular points on the plane of the axio-dilaton. From this point of view, we pick up 10 points in the descending order of the probability of the stable vacua on the moduli space of the axio-dilaton. Then, we discuss their consequences to the distributions of the complex structure modulus as summarized in Tables 9 and 10. It turns out that the residual Z 3 fixed point is favoured in the landscape, but the existence of other points is sensitive to the value of the axio-dilaton. For example, Z 2 fixed point τ = i can be realized on specific points in the moduli space of the axio-dilaton such as S = i, 2i, at which the Z 2 fixed point is only allowed in the small N max flux regime. Also, τ = √ 2i, √ 3i vacua are statistically favored on specific values of the axio-dilaton S.  Table 7: Probabilities of the stable vacua as functions of (Re S, Im S) in the descending order of the probability for N max flux = 10.
(Re S, Im S)  Table 8: Probabilities of the stable vacua as functions of (Re S, Im S) in the descending order of the probability for N max flux = 1000.
(Re τ , Im τ ) (Re S, Im S) All (− 1 2 ,  Table 9: Probabilities of the stable vacua as functions of (Re τ, Im τ ) and (Re S, Im S) in the descending order of the probability with respect to (Re τ, Im τ ). Here, we list the probability for N max flux = 10 and the column in "All" corresponds to the probability of Table 2 .
(Re τ , Im τ ) (Re S, Im S) All (− 1 2 ,  Table 10: Probabilities of the stable vacua as functions of (Re τ, Im τ ) and (Re S, Im S) in the descending order of the probability with respect to (Re τ, Im τ ). Here, we list the probability for N max flux = 1000 and the column in "All" corresponds to the probability of Table 3 .

Distributions of CP-breaking vacua
Finally, we study the distribution of CP-breaking vacua. It is known that the 4D CP symmetry is embedded into the 10D proper Lorentz symmetry in terms of the 4D parity and the 6D orientation reversing transformations [51]. Here, we assume that the 10D spacetime is given by the product of 4D spacetime and 6D Calabi-Yau manifolds including the torus. Since the 6D orientation reversing transformation changes the sign of the 6D volume form dV , it corresponds to the anti-holomorphic involution of the holomorphic three-form of 6D space. Given the holomorphic three-form Ω = dz 1 ∧ dz 2 ∧ dz 3 in the local coordinates of the 6D space {z 1 , z 2 , z 3 }, the 6D orientation reversing transformation corresponds to due to the property iΩ ∧Ω ∝ dV . (See for more details in Ref. [20].) In the factorizable torus T 6 and its orbifolds, the 6D orientation reversing transformation is provided by the antiholomoprhic involution of the complex structure, Hence, the 4D CP symmetry is unbroken at Re τ i = 0. Furthermore, there exist the additional CP-conserving region in the fundamental region of the moduli space in models with modular symmetry [19] 12 .
As the simplest example having the CP symmetry, Re τ i = ±1/2 are the CP-conserving lines, since they are invariant under Eq. (57) up to the T -transformation, namely (58) 12 In Ref. [30], it was shown that the CP phase can be rephased concretely in modular flavor models.
In addition to it, the locus |τ i | = 1 is also CP-invariant up to the S-transformation. On this locus, the S-transformation is identified with the CP-transformation; Thus, the CP symmetry remains on Re τ i = 0 and the boundary of the SL(2, Z) fundamental domain, that is, the boundary of the GL(2, Z) fundamental region. This argument is applicable to the moduli space of the axio-dilaton. So far, we have considered the factorizable 6-torus with different complex structure moduli, but the above property also holds for the isotropic torus τ = τ 1 = τ 2 = τ 3 . The breakdown of the CP symmetry is of particular importance to determine the size of Cabbibo-Kobayashi-Maskawa and Maki-Nakagawa-Sakata phases through the moduli-dependent Yukawa couplings. Indeed, in both bottom-up and top-down approaches (such as magnetized D-branes in Type IIB string setup), Yukawa couplings depend on the complex structure moduli τ . We examine the distributions of CP-breaking supersymmetric vacua, allowing for general RR and NS fluxes which generically break the CP invariance of the action. Hence, it is not a scenario of spontaneous CP violation, but the CP symmetry would be spontaneously broken by particular choices of fluxes. Given the analytical solution satisfying the tadpole cancellation condition (35), we find that the axionic fields (Re τ, Re S) are mostly stabilized at the line Re τ = 0 or the boundary of the SL(2, Z) moduli spaces, especially in the small N max flux case as shown in Figure 3. The probabilities of the stable vacua as functions of (Re τ, Re S) are summarized in Table 11 for N max flux = 10 and Table 12 for N max flux = 1000, respectively. Hence, we conclude that CP-breaking vacua are disfavored in the complex structure moduli sector. This result is in agreement with the statistical approach. From the number of vacua as function of the vacuum expectation value of τ in Eq. (55), m = 0 vacua are statistically favored due to the suppression of the denominator. Since the m = 0 vacua correspond to Re τ = 0 in Eqs. (46) and (47), the CP-conversing vacua are favored in the context of the statistical approach. Furthermore, Re S seems to be more likely to lead CP-breaking vacua than Re τ , indicating from Table 13.
On the other hand, the CP-breaking vacua in Eq. (53) have O(1)% of probabilities. In particular, when Re S = ±1/4, the CP-breaking vacua for the complex structure modulus, are statistically favored in the landscape. This correlation may be important from the viewpoint of CP-violation phenomena. Also, Figure 3 shows the linear correlations, i.e., Re τ = ±Re S. These results suggest that different CP-violating phases such as CP phases in mass matrices and the strong CP phase may have a common origin 13 . We comment on the reason why the vacua in Eq. (60) is favored in the complex structure modulus from the viewpoint of the statistical approach. We come back to the denominator in Eq. (55), which determines the statistically favored moduli values. This denominator is rewritten in terms of the vacuum expectation value of the modulus field (46).
which indicates that the small l enhances the number of vacua, in particular for l = ±1. In that l = ±1 case, together with the relation the integer m also requires Hence, the existence of CP-breaking vacua Re τ = ±1/4 are also favored in the statistical approach. (Re τ, Re S)  Table 11: Probabilities of the stable vacua as functions of (Re τ, Re S) in the descending order of the probability for N max flux = 10.
(Re τ, Re S)  Table 12: Probabilities of the stable vacua as functions of (Re τ, Re S) in the descending order of the probability for N max flux = 1000.

Modular symmetric flavor models
In this section, to illustrate implications of our results, we compare the predictions in modular A 4 models with the distributions of complex structure modulus τ in flux vacua as shown in Section 3. In particular, we focus on 8 classes of models in Ref. [33] as summarized in Section 4.1. Section 4.2 is devoted to the distributions of the complex structure modulus τ without fixing the axio-dilaton S. Furthermore, these distributions predict favorable vacuum expectation values of the moduli fields consistent with the observational data in 8 classes of modular A 4 models.

Modular A 4 models
We briefly review modular A 4 models explaining the flavor structure. Before going to the detail of these models, we recall that there exist finite non-Abelian groups in the SL(2, Z) modular group. The principal congruence subgroup Γ(N ) of Γ = SL(2, Z) and its quotient groups are summarize as follows: Since the complex structure is invariant under the transformation S 2 , the complex structure moduli space is governed byΓ ≡ Γ/{±I}. Its quotient groups are also defined as and Γ 2,3,4,5 are isomorphic to the phenomenologically attractive S 3 , A 4 , S 4 , A 5 non-Abelian discrete groups, respectively. To discuss the Yukawa couplings of quarks and leptons, modular form f (τ ) of weight k and level N is of particular importance to be defined to satisfy where p, q, s, t denote the element of Γ(N ). In particular, one can assign the standard model particles to be irreducible representations of the modular group Γ N . When we denote F (τ ) = (f 1 (τ ), f 2 (τ ), · · · ) T represent multiplets of modular forms transforming in the irreducible representations of Γ N , it obeys for γ ∈Γ (1), where ρ(τ ) denotes a certain unitary representation matrix.
In this paper, we focus on the Γ 3 A 4 modular group which corresponds to a minimal subgroup admitting a triplet representation as well as three singlets. The singlet representations have the following representations under S-and T -transformations: whereas the triplet representation is assigned to be transformed under S-and T -transformations: Since the modular forms of the weight 2 are spanned on three-dimensional modular space, these modular forms of the weight 2 are assigned to the triplet of A 4 . Explicit forms of modular forms of the weight 2 are constructed in Ref. [25], and higher weight modular forms such as weights 4 and 6 are also constructed by tensor products of the modular forms of weight 2. We follow the supersymmetric modular A 4 models proposed in Ref. [33], in which comprehensive analysis of the lepton sector is performed under the smallest number of free parameters. In this analysis, we assign the charge assignment of SU (2), A 4 and the modular weight for lefthanded lepton doublets L = (L 1 , ) and the Higgs doublets H u , H d as shown in Table 14.
The modular weights of those fields, −k I , are free parameters, but they are constrained to obtain the modular-invariant superpotential. Since the Yukawa couplings transform as the unitary representation ρ Y of the modular group Γ N , the superpotential is modular invariant only if the modular weight of the coupling Y i 1 ···in , k Y , satisfy k Y = n k in (in global supersymmetric theory) with k in being the modular weight of fields Φ in . Furthermore, the product of representation ρ Y (γ) of the Yukawa coupling and that of fields should include the singlet representation, namely ρ Y ⊗ n ρ in 1.
To realize the tiny neutrino masses, the authors of Ref. [33] discussed two scenarios. First one is to introduce the Weinberg operator W = 1 Λ (H u H u LLY ) 1 with cutoff scale Λ, leading to ten models, A i (i = 1, · · · , 10) and the other is to consider the type-I seesaw mechanism, leading to 30 models B i , C i , D i (i = 1, · · · , 10). From the totally 40 classes of these models, it was found that 14 models are fitted with the observational data for both normal ordering and inverted ordering neutrino mass spectra. We compare these phenomenological results with the distribution of the complex structure modulus τ predicted from the string landscape in Section 3.

Distribution of complex structure modulus in modular A 4 models
We are now ready to examine predictions of modular A 4 models in light of our results on the moduli stabilization. As for the modular A 4 models, we focus on the successful 8 models labeled by B 9 , B 10 , D 5,6,7,8,9,10 , in which the type-I seesaw mechanism are employed to realize the tiny neutrino masses with the normal ordering. In these models, 6 (real) free parameters and two overall mass scales successfully lead to explain the flavor structures of the lepton sector, that is, the mass hierarchy of charged lepton masses, differences of neutrino masses squared and mixing angles.
The allowed regions for (Re τ, Im τ ) in 8 models are summarized in Table 15, which provide the constrained observables such as the mixing angles of the lepton sector, relatively large neutrino masses and maximal Dirac CP phase. Note that the analysis in Ref. [33] focuses on the following moduli space: That means that the Z 2 fixed point (Re τ, Im τ ) = (0, 1), the   Table 15: Phenomenologically allowed regions of (Re τ, Im τ ) [33].
Ref. [33], although these points lead to high probabilities in our results from the string landscape. The allowed regions of 4 ) points, respectively. In Figure 4, we append these phenomenologically predicted moduli values to the distributions of the complex structure modulus discussed in Section 3.1. 14 Since the string landscape predicts discrete vacuum expectation values of the modulus field at the finite number of vacua, the phenomenologically allowed region is tightly constrained in the string landscape. Indeed, in the small N max flux = 10 case, the boarder of the phenomenologically acceptable region is only allowed for B 9 , B 10 , D 5,6,7,8 models, and the other regions are rejected for all models. When we increase N max flux , other phenomenologically acceptable regions are allowed. Taking into account the results of Figure 4, the vacuum expectation values of the modulus field with the highest probability in the string landscape for 8 classes of modular A 4 models are shown in Table 16 for N max flux = 10 and Table 17 for N max flux = 1000, respectively. These moduli vacuum expectation values correspond to the theoretically favored values in the string landscape. The probabilities of these favored moduli vacuum expectation values in (i) each model and (ii) whole fundamental domain are also displayed in Tables 16 and 17, respectively. Furthermore, (iii) the probability of the phenomenologically allowed region for (Re τ, Im τ ) inside the string landscape is also shown in Tables 16 and 17. These tables indicate that the realistic modulus values τ in B 9 , B 10 , D 5,6,7 models are likely to be realized in the string landscape, in comparison with D 8,9,10 models.
We remark about our analysis in the modular A 4 models. First, we set the size of each bin as (0.01, 0.01) in Figure 4 just for visibility, but in Tables 16 and 17, we show moduli vacuum expectation values up to three decimal places in the calculation of various probabilities to compare with the results of Ref. [33]. Next, there are some strongly favored points on the unit circle as seen from Tables 16, 17. However, we excluded such points to compare with the results of [33], but in general if an allowed region contains such point on the unit circle |τ | = 1, it would become highly favored region in the string landscape. Indeed, the allowed regions of D 8 , (B 10 , D 6 ) and (B 10 , D 5 ) models allowing values around |τ | = 1 are favored around these region as mentioned before. In particular, the D 8 model predicts observed values close to the Z 3 fixed point. (See for the phenomenological studies around these points, Refs. [27,[48][49][50].) From the viewpoint of the CP violation, the value Re τ = 1/4 is important. Thus, most of models, B 9,10 and D 5,6,7 include such points. These models may be interesting from CP phenomenology.
-0. 50  None --0.000 Model D 10 None --0.000 Table 16: The theoretically favored vacuum expectation values of (Re τ, Im τ ) within each region for the N max flux = 10 case. We also show the probabilities of (i) the favorable value in each model, (ii) the favorable value in whole fundamental domain of SL(2, Z) moduli space and (iii) the allowed moduli region in the string landscape.  Table 17: The theoretically favored vacuum expectation values of (Re τ, Im τ ) within each region for the N max flux = 1000 case. We also show the probabilities of (i) the favorable value in each model, (ii) the favorable value in whole fundamental domain of SL(2, Z) moduli space and (iii) the allowed moduli region in the string landscape.

Conclusion
We have studied the moduli stabilization by flux compactifications on T 6 /(Z 2 × Z 2 ) orientifold from the viewpoint of modular flavor symmetries. We have systematically analyzed stabilized moduli values in possible configurations of flux compactifications, investigating probabilities of moduli values and showing which moduli values are favorable from our moduli stabilization. Then, we have examined their implications on modular symmetric flavor models.
In the bottom-up approach, the modular and CP symmetries lead to the construction of phenomenologically attractive models. Since Yukawa couplings as well as the higher-order couplings are controlled by the modular form of finite discrete subgroups of SL(2, Z), the stabilization of complex structure moduli appearing in the Yukawa couplings is of particular importance to compare the predictions of modular flavor models with the observational data. These moduli values determine not only the flavor structure of quarks and leptons but also the CP-violating phases through the Yukawa couplings written by the modular forms. So far, the random search in the fundamental domain of SL(2, Z) have often been carried out in the bottom-up approach. Our approach is to predict these moduli values and their probability distributions in the string landscape.
It turns out that values of the complex structure modulus τ are clustered at the Z 3 fixed point τ = ω in the fundamental region of SL(2, Z) τ associated with the geometrical symmetry of the torus. Especially, the Z 3 fixed point is statistically favored in the weak string coupling and small flux regime, although the distribution of the complex structure modulus is correlated with that of the axio-dilaton S. It is possible to realize the other Z 2 fixed point by focusing on specific vacuum expectation values of the axio-dilaton such as the Z 2 fixed point in the fundamental region of the SL(2, Z) S moduli space associated with the axio-dilaton. In general, the modulus values with definite probabilities are intersecting points between |τ | 2 = k/2 and Re τ = 0, ±1/4, ±1/2.
From both the top-down and bottom-up approaches, it is quite interesting that the Z 3 fixed point τ = ω as well as the Z 2 fixed point is favored. For example, when τ = ω, the four fixed points on the T 2 /Z 2 orbifold are equally spaced like a tetorahedron. That is quite symmetric. That affects the couplings among localized modes on fixed points, e.g. twisted strings in heterotic orbifold models, leading to the same strengths of couplings among different modes. That is, remnant of such a geometrical symmetry appears in string-derived effective field theory. Phenomenologically the residual Z 3 and Z 2 symmetries can appear in mass matrices of the bottom-up approach model building, constraining their forms. In this respect, our string landscape supports analyses in modular flavor models near these fixed points. In addition, it is interesting to extend such phenomenological analyses to other intersecting points between |τ | 2 = k/2 and Re τ = 0, ±1/4, ±1/2.
It is found that CP-breaking stable vacua are statistically disfavored in the complex structure modulus sector, meaning that the complex structure modulus is mostly stabilized at Re τ = 0 or the boundary of the SL(2, Z) τ moduli space. Among CP-breaking vacua, the values Re τ = ±1/4 are most favorable. Although the probability of CP-breaking vacua is small compared e.g., with the Z 3 fixed point τ = ω, its probability increases for certain values of the axiodilaton, e.g., Re S = 1/4. That means a strong correlation between CP phases due to Re τ and Re S. That suggests the possibility that different CP phases in the 4D effective field theory, e.g., CP phases in fermion mass matrices and the strong CP phase, may be related with each other in superstring theory.
As for the phenomenological implication, we examine the modular A 4 models in Ref. [33], where the lepton sector has the irreducible representations of A 4 realized as a principle congruence subgroup of SL(2, Z) τ . For the 8 phenomenologically attractive models with the normal ordering neutrino mass spectra, we compare these phenomenological predictions with the distribution of the complex structure modulus. We find that these phenomenological results are tightly constrained in our results of the string landscape in which the complex structure modulus is distributed at the discrete vacua as shown in Fig. 4. In this paper, we focused on modular A 4 models but our results in Section 3 are applicable to other modular flavor models. It is interesting to apply our method to other modular flavor models including concrete flavor models in string-derived effective field theories 15 . mum conditions taking into account the above perspective, let us redefine RR and NS fluxes {a 0 , a, b, b 0 , c 0 , c, d, d 0 } into the following integer variables {m, l, n, u, v, s, r} as provided in Section 3 16 . The relation between two sets is expressed as : rl = a 0 , rm + sl = −3a, rn + sm = −3b, sn = −b 0 , ul = c 0 , um + vl = −3c, un + vm = −3d, vn = −d 0 .
In fact, the relation is not a one-to-one correspondence and this is problematic for counting vacua correctly. One need to exclude sets {m, l, n, u, v, s, r} which gives same {a 0 , a, b, b 0 , c 0 , c, d, d 0 } carefully as discussed later.
The D3-brane charge is now represented as As a result of the flux quantization, this implies that N flux is a multiple of 64. Furthermore, the left-hand side becomes a multiple of 9 × 64 [38], and then N flux becomes a multiple of 192, which is not clear at this time. These relatively-large D3-brane charge would cause the backreactions to the geometry. However, we do not consider these backreactions in this paper.
• T q ∈ SL(2, Z) τ As discussed in Section 2, we can focus on only the l = 0 case. Hence, there is the equivalence relation: i.e., there are only 2l independent choice of m, and we can freely set m as m = −l, −l + 1, . . . , l − 1.
• T q ∈ SL(2, Z) S In this case, v = 0 is possible but u = 0 = v leads unstabilized S. Hence we have to consider these two cases: r ∼ r + qu (v = 0).
From now, we focus on the l, n > 0 case. The l, n < 0 case is completely analogous, so that we simply omit it. However, we included both results in the actual calculation. Recall that, since we require τ given in Eq. (46) to be complex, m 2 − 4ln must be negative in the l, n > 0 case. Furthermore, we can impose more severer condition on it: Since the inside of the unit circle and its outside are connected by S-transformations and Ttransformations have been already fixed, it is enough to search for Im τ ≥ √ 3 2 region. Let us come back to Eq. (74). Since rv − su must be integer, holds. Taking into consideration that N flux ≥ 0 holds for the supersymmetric solutions, we conclude l = 1, . . . , N flux , To constrain other fluxes {r, s, u, v}, we look at the expression of S in Eq. (44). As with the τ case, it is enough to require ImS ≥ By expanding |uτ + v| 2 = (uRe τ + v) 2 + u 2 (Im τ ) 2 and reducing the above condition for u, v, we find that with Re τ = − m 2l . As discussed before, we have to divide in two cases, namely (i) v = 0 and (ii) v = 0.
Then one can obtain r from the expression of N flux (74) as • (ii) v = 0 case Similarly, one can obtain s from N flux as and we adopt r ∼ r + qu so that r = 1, 2, . . . , |u| − 1.
Note that it is not clear that whether r, s are integer or not at this stage, therefore one has to check the integrability constraints in the actual calculation. Lastly, we summarize the protocol to obtain all the physically-distinct solutions: