Anomaly free Froggatt-Nielsen models of flavor

We introduce two anomaly free versions of Froggatt-Nielsen (FN) models, based on either $G_{\rm FN}=U(1)^3$ or $G_{\rm FN}=U(1)$ horizontal symmetries, that generate the SM quark and lepton flavor structures. The structure of these"inverted"FN models is motivated by the clockwork mechanism: the chiral fields, singlets under $G_{\rm FN}$, are supplemented by chains of vector-like fermions charged under $G_{\rm FN}$. Unlike the traditional FN models the hierarchy of quark and lepton masses is obtained as an expansion in $M/\langle \phi\rangle$, where $M$ is the typical vector-like fermion mass, and $\langle \phi\rangle$ the flavon vacuum expectation value. The models can be searched for through deviations in flavor observables such as $K-\bar K$ mixing, $\mu\to e$ conversion, etc., where the present bounds restrict the masses of vector-like fermions to be above ${\mathcal O}(10^7~{\rm GeV})$. If $G_{\rm FN}$ is gauged, the models can also be probed by searching for the flavorful $Z'$ gauge bosons. In principle, the $Z'$s can be very light, and can be searched for using precision flavor, astrophysics, and beam dump experiments.


Introduction
The masses of the Standard Model (SM) quarks and leptons exhibit large hierarchies [1]. The mass of the lightest quark, the up quark, is five orders of magnitude smaller than the mass of the heaviest quark, the top quark. The hierarchies in the leptonic sector are even more striking; the heaviest charged lepton, the tau, has a mass that is at least ten orders of magnitude larger than the mass of the heaviest neutrino. Furthermore, in the quark sector the measurable misalignment between the gauge and mass-eigenstates is also hierarchical, as encoded in the CKM matrix. For leptons, on the other hand, the three mixing angles in the PMNS matrix are all large.
This very peculiar flavor structure calls for a dynamical explanation. A very attractive solution to the above SM flavor puzzle is provided by the Froggatt-Nielsen (FN) models of flavor [2] (see also [3,4]). In FN models the SM mass hierarchy is due to a horizontal U (1) FN , spontaneously broken by the vacuum expectation value (vev) of the flavon field, φ . The SM fermions carry horizontal U (1) FN charges. The up and down quark Yukawa matrices are then given by 1 . For generic horizontal charges, and without additional field content, such a horizontal U (1) FN is anomalous and thus needs to be a global symmetry. It is also possible to gauge the FN horizontal symmetry for special choices of horizontal charges for which the symmetry is not anomalous. For example, Ref. [5] builds neutrino mass models based on an anomaly free U (1) 2 FN FN symmetry, while the full list of anomaly free charge assignments for U (1) FN can be found in [6,7] (for related work, see also [8][9][10][11][12][13][14]). The gauged U (1) FN can also be left naively anomalous, if the cancellation of U (1) FN gauge anomalies occurs through the Green-Schwarz mechanism [15][16][17][18][19][20][21][22][23][24][25][26].
In this paper we explore a twist to the original FN proposal. In our modified set of FN models the horizontal U (1) FN is anomaly free for any horizontal charge assignments, [q i ], [u i ], [d i ], making it much easier to obtain agreement with the measured fermion masses and mixings. In this case up and down quark Yukawa matrices are given by . (1.2) In these "inverted FN models" the SM mass hierarchy is obtained for flavon vev that is parametrically larger than the vector-like masses, e.g., for φ /M ∼ 5. Such a realization of the FN flavor model was discussed in Ref. [27], where it was given as a possible UV realization of the clockwork mechanism [28]. We extend the results of Ref. [27] in several ways. First of all, two choices for the horizontal symmetry group, G FN = U (1) 3 FN and G FN = U (1) FN , are explored. The choice G FN = U (1) 3 FN most closely resembles the analysis in Ref. [27], except that we allow O(1) variations for all the Yukawa couplings. We also assume that G FN is gauged, and explore the phenomenology of the associated Z gauge bosons.
The typical spectra of inverted FN models for the two choices of G FN are shown in Fig. 1. The flavor constraints require the vector-like fermions and the flavon to be heavy, for instance, the vector-like fermions need to be heavier than about O(10 7 GeV). The mass of Z , however, can span many orders of magnitude. For gauge couplings g ∼ O(1) the Z is heavy, while for small values of g it can be light and accessible at experiments. A better part of the present manuscript is devoted to analyzing the experimental possibilities for the light Z searches.
The paper is organized as follows. In Section 2 we present the inverted FN model based on U (1) 3 FN symmetry, starting with the mass suppression mechanism for a single quark generation, which is then readily extended to the case of three generations. Here we also explain that the scaling (1.2) can be understood using a spurion analysis for an approximate U (1) app that emerges at low energies, and is broken by M rather than φ . In subsection 2.2 we then perform a numerical scan over the parameters of the U (1) 3 FN model and compare the results with the measured values of quark masses and mixings. In Section 3 we present the U (1) FN model and the respective numerical scan. In Section 4 we extend both models to the lepton sector. In Section 5 we explore the phenomenological implications of the two models, covering the bounds from a number of flavor conserving and violating processes induced by the exchanges of Z i gauge bosons, with implications for precision flavor, astrophysics and beam dump experiments. Section 6 contains our conclusions, while Appendix A gives further details on K −K bounds for flavorful light Z . The representative benchmarks for the two models are given in Appendix B.

The set-up of the model
We start by assuming that there is a separate horizontal U (1) FN for each generation of the SM fermions, i.e., that the horizontal gauge group is G FN = U (1) 1 × U (1) 2 × U (1) 3 . The field content is that of the SM, supplemented by a number of vector-like fermions whose FN charges differ by one unit and are organized in chains of vector-like fermions. We first illustrate the set-up for the case of a single SM generation of down quarks, setting also G FN = U (1) FN for now. The Lagrangian is  In general, there could also be terms of the form φ * d L,n−1 d R,n and φ * q R,n−1 q L,n , which were dropped in (2.2) and (2.3). The reason for this choice is that we expect the above model to be supersymmetrized at some higher scale. Such terms are then absent because the superpotential is holomorphic.
The flavon obtains a vev, φ , which spontaneously breaks the horizontal U (1) FN . Each of the chains of vector-like fermions, q and d, has one massless zero mode, and N q and N d massive vector-like fermion states with mass O(M, φ ). For φ M , with Y d n , Y q n ∼ O(1) the two zero modes are mostly localized towards the ends of the respective chains. More precisely, the q L,0 and d R,0 massless mass eigenstates are given by (2.4) The notation is borrowed from [27], with V q L and V q R the N q × N q and (N q + 1) × (N q + 1) unitary matrices that diagonalize the N q × (N q + 1) mass matrix for the q fields through q L,N q q L,N q −1 q L,1 q L,0 d L,1 d L,  a bi-unitary transformation (and similarly for d L,R fields with obvious change of notation). For Y d n = Y q n = 1, M d n = M q n = M, the model realizes the clockworking mechanism, for which the zero mode profiles are known analytically [27], The zero mode obtains a mass once the electroweak symmetry is broken through the Higgs vev. The mass is suppressed by the zero mode wave functions on the zero node, where the overlaps with the zero-node are given by and are thus exponentially suppressed for large N ψ and φ M . The suppression can be understood from a spurion analysis. In the limit of a heavy flavon the fermion mass matrix has an emergent approximate U (1) app global symmetry after the flavon obtains a vev. This is easy to see once the fermion fields are re-arranged in a modified chain of nodes as shown in Fig. 3, while the heavy flavon is integrated out and need not be considered here. The fields on a particular node have U (1) FN charges that differ by one unit, but have the same U (1) app charges. The U (1) app charge reduces by one unit when hopping to the right by one node in the chain. That is, the q L,Nq has the U (1) app charge N q , the q L,Nq−1 andq L,Nq the U (1) app charge N q − 1, and so on, all the way to d R,N d which carries a charge −N d (the two greyed out nodes in Fig. 3 should be treated as a single node, with all the corresponding fields uncharged under U (1) app ). The mass terms proportional to M explicitly break U (1) app . Treating M as a spurion with U (1) app charge q L,N q q L,N q −1 Figure 3. The anomaly free inverted Froggatt-Nielsen model for a single generation of down quarks, but with rearranged fields on the nodes, see text for details.
of +1 the fermion mass matrix is formally invariant under U (1) app . The zero modes, i.e., the SM fermions are mostly q L,Nq and d R,N d , so that the zero mode mass term is given by which is U (1) app invariant. This is exactly the suppression found in (2.6). Furthermore, note that the effective FN charge of the zero mode is approximately the same as its U (1) app charge. The zero modes are localized toward the ends of the vector-like chains, see Fig. 5.
To leading order the q L,Nq and d R,N d zero modes thus carry effective FN charges N q and −N d , respectively.
The above discussion is readily extended to the case of three generations, including both down and up quarks. The horizontal gauge group is taken to be G FN = U (1) 1 × U (1) 2 × U (1) 3 , with a separate U (1) i for each of the three generations of fermions. The fermion Lagrangian is then (2.9) with i, j = 1, 2, 3, the generation indices, andH = iσ 2 H * . The kinetic Lagrangians are given in (2.2) and (2.3), but with extra generational superscript on fermions, q n , N q → N q(i) , etc., as well as φ → φ i , since each generation is gauged under a separate U (1) i . This results in three uncoupled chains, one per each generation, for q fields, for d fields, and for u fields. All of these chains attach to the central node with the Higgs, see Fig. 4. The Yukawa matrices Y d 0 , Y u 0 on the zero node are 3 × 3 complex matrices, while the Yukawa couplings on each of the chains, (2.2) and (2.3), are just complex numbers. 2 The Yukawa interactions with the Higgs are thus the only terms mixing different generations, see the last term in (2.9). 2 The flavon Yukawa couplings Y q(i) 1 , connecting the zero node with the first node of the q−chain, are in general arbitrary complex three-vectors. To simplify the notation in (2.9) we take them to be orthogonal, and then use the 3 × 3 unitary field redefinitions ofq  FN . The chiral fermions are on the zero node, denoted by the red ellipse, which also contains the Higgs. Up to higher corrections in v/M the SM fermions are equal to the zero modes. The zero mode wave functions peak toward the ends of the chains, cf. Eq. (2.5) and Fig. 5. This means that the SM fermion, q (i) L , has the largest component that is due to q (i) L,Nq , the left-handed part of the vector-like fermion that is the furthest away from the Higgs. For u R,Nu , and similarly for d R,N d . The SM fermions thus carry large U (1) i charges, and couplings to the U (1) i gauge bosons, Z i s, roughly given by the length of the corresponding vector-like chain, [ψ] i ≈ N ψ(i) . However, for the phenomenology of Z i also the subleading couplings are important, because they lead to flavor violating transitions. We explore the implications of these in Section 5.
The zero modes, identified with the SM fermions, obtain a nonzero mass after electroweak symmetry breaking from the overlaps with the Higgs. The Higgs is on the zero node, so that (2.10) The structure of the zero node overlaps, f a , explains the hierarchy of SM fermion masses. For M φ the zero node overlaps, f a , are more suppressed the longer the corresponding  chain of the vector-like fermions. This can be seen analytically, if all the Yukawa couplings are equal to one, see Eq. (2.7), but is also true in general.
Below we perform a scan over the inputs of the model in order to explore how well the eigenvalues of the quark mass matrices (2.10) resemble the observed values. To simplify the discussion, we set where M 1 , M 2 , M 3 are the typical values of vector-like masses for each of the three generations. This choice will also make it easier to compare with the single U (1) FN case, to be covered in Section 3. The scalings for the mass matrices in (2.10) are then, for (2.12)
The experimental values of quark and charged lepton masses (in GeV) at µ = 10 7 GeV, obtained from NNLO QCD RG evolution, and used for comparison with the numerical scan. Despite the large allowed experimental range for m u and m d we show only the central values, which coincide with the choice made in Ref. [27] at µ = 2 TeV.
form r a e iϕa . The Yukawa couplings between the FN fermions and the flavon are taken to be of the form Y f n φ = r a e iϕa q, so that the ratios |Y f φ /M f | are on average equal to q. The magnitudes and phases are taken to be uniformly distributed over r a ∈ [0.3, 0.9] and ϕ a ∈ [0, 2π), with r a , ϕ a uncorrelated between different couplings. The range of r a was chosen such that the top Yukawa,m t /(v/ √ 2) 0.55 is close to the median, while the ratio of the boundaries, r max /r min , is smaller than q. The predictions are not very sensitive to the precise ranges for r a since these cancel on average in the ratios M f n /(Y f n φ ), which control the hierarchies of the SM quark masses. The effect of changing the average value of r a is thus mostly due to a different average values of the Y d,u 0 matrix elements. For instance, using r a ∈ [0.6, 1.8] changes the distributions for the quark masses in Fig. 2.2 by an overall factor of 2. Relatively larger r a ranges, on the other hand, lead to wider distributions of quark masses and mixings.
The lengths of the chains with the vector-like quarks are taken to be almost the same as in Ref. [27]  In each run of the scan we generate and diagonalize the chains, then calculate the quark masses and the CKM matrix elements. The results of the scan with 5000 runs are shown in Fig. 6, with the CKM matrix elements shown in the first row, and the quark masses in the second and third row. For ease of comparison we also denote the measured central values with vertical solid lines, while the one sigma experimental bands are delineated with dashed lines. The quark masses are the MS masses at µ = 10 4 TeV, which we anticipate to be the rough lower bound on the vector-like fermion masses from flavor constraints, see Section 5 for details. The values of SM fermion masses at µ = 10 4 TeV, obtained through NNLO RG evolution, are listed in Table 1 (at µ = 2 TeV the values of the quark masses coincide with Ref. [27]). For the CKM matrix elements we checked using results in [29] that the effects of RG evolution are negligible, so that we use the unevolved values from the PDG [1]. 3 The difference arises partially from taking significantly higher µ and different median values of Y d,u 0 matrix elements. We take N d(3) = 3, N u(2) = 2 while in [27] these were set to N d(3) = 2, N u(2) = 1. We thank A. Kagan for suggesting the new charge assignments that improve the agreement with bottom and charm quark masses. Higgs u (1) u (3) q (3) d (3) d (2) d (1) q (1) q (2) u (2) and similarly for L u with the d → u replacement. For L q , on the other hand, there is one more left-handed field for each generation than the right-handed fields, so that The summation is over the nodes and the generations on each node. We label the fermions such that the i−th generation fermions have a vector-like chain of length N q(i) . The first generation has the longest fermion chain, so that the summation n = 1, . . . , N q(1) sums over all of the nodes. On node n there areN q | n generations, which are coupled to theN q | n−1 generations of fermions on the (n − 1)-th node through anN q | n−1 ×N q | n complex Yukawa matrix Y q n . This produces a set of coupled FN chains, where the mixing between generations is not only through 3 × 3 Yukawa matrices on the zero node, coupling to the Higgs, but also through the flavon Yukawa couplings, Y q n ji , and vector-like masses, M q n ji , cf. Fig. 7. The zero modes are obtained through unitary transformations and similarly for the up quark zero modes. As in the decoupled case, the Higgs provides the zero modes with nonzero masses after electroweak symmetry breaking from their overlaps with the zero node. However, in the coupled case the zero mode overlaps with the zero node are described by 3 × 3 matrices, V d R ,u R 0(j),0(i) , due to the inter-generational mixing in the chains. Therefore, the quark mass matrices become For φ M the zero-modes are, also in this case, mostly localized towards the ends of the respective FN chains, with small overlaps with the Higgs on the zero node. Parametrically, with M q,u,d , Y q,u,d denoting the typical values of the corresponding matrix elements. The hierarchy in the overlaps then translates into the hierarchy of the quark masses, cf. Eq. (2.10). However, due to additional flavour mixing on each node the parametric relations are now even more approximate compared to the case of decoupled of FN chains. We show this by performing a numerical scan. Since the expressions for SM quark masses effectively involve multiplications of a number of random matrices, Eq. (3.6), part of the hierarchy comes from the properties of random matrix multiplications [30] (see also Section 4).

Numerical scan
We choose the chain configurations with the following number of fermion generations on each nodeN Here n = 0, . . . , 3, respectively, for each of the three cases. This implies that the lengths of the vector-like chains are given by,   where N a + 1 is the length of the corresponding chain, cf. Eq. (3.2). The above chain configuration differs from the one for the decoupled FN chains in (2.13), in that the all three d R generations now have the same lengths of chains, and that the chain for the first u R generation is shorter.
We use a similar procedure as in Section 2.2 to produce the numerical scan for the coupled case. We fix q = 5, φ = 10 7 GeV, and vary the matrices . The vector-like masses are taken to be random complex numbers of the form r a e iϕa in the units of 10 7 GeV, while the elements of Yukawa matrices are taken to be of the form (Y f n ) ij φ = r a e iϕa q. The r a are varied in the range [0.3, 0.9] and phases in [0, 2π). For each scan run we generate and diagonalize the coupled chains and extract the quark masses and the CKM matrix elements. The results of the scan with 5000 runs are shown in Fig. 8, using the same layout as in Fig. 6.
Even though the d R chains are of the same length for all generations, we still have hierarchical masses for down quarks due to different lengths of q L chains. Part of the required hierarchy also comes from the fact that products of random matrices lead to matrices with hierarchical eigenvalues [30]. The resulting mass distributions are much broader than in the U (1) 3 FN case. This is particularly apparent in the distributions of CKM matrix elements, where the distributions for the V cd and V cb matrix elements almost • "|Vus|" • "|Vub|" • "|Vcd|" • "|Vcs|" • "|Vcb|" • "|Vtd|" • "|Vts|" • "|Vtb|"  Figure 9. The CKM elements that would be obtained from coupled FN chains for the case of a very large q = 1000. This shows the hierarchy between the would-be CKM matrix elements more clearly.
completely overlap. This is a result of relatively small expansion parameter 1/q 0.2, required to fit the observed values realized in Nature. In Fig. 9 we show that a hierarchical structure does appear also for the CKM matrix elements once a much larger value of q = 10 3 is taken.

Extension to leptons
It is straightforward to extend the above framework to leptons. To shorten the discussion we assume that the neutrinos have Majorana masses. For the G FN = U (1) 3 FN case the Lagrangian for the leptons is then given by replacing Q → L, d → e in Eq. (2.9), and ignore terms that involve u. Performing the same replacements in Eq. (3.1) gives the Lagrangian for the G FN = U (1) FN case.
For the neutrino masses we assume that they come from the dimension 5 Weinberg operator. On the zero node we thus add the mass term with c ij ∼ O(1). In the φ M L n limit the neutrino mass matrix takes the form after the vector-like fermions are integrated out. The FN charge assignments for charged leptons depend crucially on the assumed flavor structure of the neutrino mass matrix, see, e.g., [31]. We focus on the case where the neutrino masses are completely anarchic, with all the PMNS mixing angles taken to be O(1). This happens if all the left-handed leptons have the vector-like fermion chains of the same length, FN the hierarchy among charged lepton masses is then due to different lengths of FN chains for the right-handed leptons, giving .

(4.4)
The observed hierarchy between m e : m µ : m τ is obtain for Even with this identification, there is still significant freedom in phenomenologically viable charge assignment, since one can compensate for a particular choice of N L(i) by adjusting globally the N e(i) . In Fig. 10 (left) we show the charged lepton masses that are obtained by setting Fig . 11 shows the resulting values of PMNS matrix elements and of the neutrino masses, setting Λ LN = 10 11 GeV in order to approximately saturate the bound on the sum of neutrino masses from cosmology i m ν i 0.15 [32]. These are compared with the measured values with solid (dashed) lines denoting the central values (1σ bands) [1]. The scan is peformed in the same way as for the quarks in the case of decoupled FN chains, Section 2.2. The vector-like masses are taken to be random complex numbers with magnitudes in the ranges r a ∈ [0.3, 0.9] in units of 10 7 GeV, with arbitrary phases, while the Yukawa couplings Y f n φ are equal to q up to a randomized complex prefactor with magnitude in the range [0.3, 0.9] and a random phase, where φ = 10 7 GeV. This preferentially leads to normal hierarchy, see Fig. 11, and to PMNS phase and Majorana phases that are completely random.
For the case of coupled FN chains, G FN = U (1) FN , we exploit the fact that products of random matrices have hierarchical eigenvalues. We find that a completely anarchic charge  [1,32,33], with dashed lines denoting 1σ errors. For i m νi the upper bound is shown [32], while for ∆m 2 31 two lines are shown: blue (orange) for normal (inverted) ordering.
describes well the hierarchy among the charged leptons, see Fig. 10 (right). In the φ M limit the charged lepton mass matrix is then For (M L,e n ) ij = (c L,e n ) ij q φ , with the Y L,e n and c L,e n matrix elements randomly distributed, such that their average values vanish, one has [30] (det • |Uτ1| • |Uτ2| Tr where for simplicity we have assumed that (c L,e n ) −1 and Y L,e n all follow the same distribution with variance σ. Here N f = 3 is the number of families, while 2(N e + N L ) + 1 is the number of random matrices that get multiplied in Eq. (4.8). From (4.9) and (4.10) we get If the lengths of the FN chains grow the eigenvalues become more hierarchical following (4.11), which we also checked numerically. The distributions of electron, muon and tau masses that follow from the anarchic charge assignments for the coupled FN chains in Eq. (4.7) are shown in Fig. 10 (right). We see that despite the anarchic charges the hierarchy among the eigenvalues indeed reproduces well the measured hierarchy of charged lepton masses. In the scan we used the same approach as for quarks in Section 3.2; the vector-like masses are random complex variables with a magnitude in the range [0.3, 0.9] and a random phase, while the elements of Yukawa matrices have entries equal to q up to a similar random complex prefactor with magnitude in the range [0.3, 0.9].
The resulting PMNS matrix elements and neutrino masses for Λ LN = 10 13 GeV are shown in Fig. 12. We observe that there is a tendency for U e,2 , U µ,3 , U τ,2 , U τ,3 to be below the observed values, and for U e,1 , U µ,2 , U τ,3 to be above. Still, the agreement with the values realized in Nature remains reasonable. For the neutrino masses the normal ordering is heavily favored as can be seen from Fig. 12 (middle and bottom). The Majorana and PMNS phases are randomly distributed.

The phenomenology of flavorful Z bosons
How can one uncover experimentally whether any of the above anomaly free FN models is realized in Nature? The immediate answer is to search for new contributions to Flavor Changing Neutral Currents (FCNCs), e.g., B −B, K −K, mixing, µ → eγ, etc. The new contributions are due to the exchanges of flavons, heavy vector-like fermions, and, if G FN is gauged, also flavorful Z's. Below we will see that the FCNC constraints on possible tree level Z exchanges bound φ O(10 7 GeV). This in turn means that the vector-like fermions are heavier than about 10 7 GeV, since for O(1) Yukawa couplings their masses are comparable to the flavon vev, φ .
We will derive the FCNC constraints on inverted FN models assuming the Z contributions dominate. That is, we will assume that both the vector-like fermions and the flavon have masses m i ∼ φ . As a result, their contributions to FCNC transitions are subleading compared to the Z contributions. The vector-like fermions contribute to FCNCs only at one-loop, while the FCNC flavon couplings to the SM fermions In the construction of inverted FN models we made two choices for the anomaly free horizontal group G FN . If G FN is gauged, there is one extra gauge boson, Z , for the G FN = U (1) FN case, and three new gauge bosons, Z i , i = 1, 2, 3 for the G FN = U (1) 3 FN case. The rest of this section is devoted to the phenomenology of these flavorful Z s. If the G FN gauge couplings are small the Z can be light. This means that beside the indirect searches using FCNCs, the Z can also be searched for in on-shell production, e.g., in beam dumps or in astrophysical environments.
We first focus on the G FN = U (1) FN case. After the U (1) FN is broken the relevant terms in the Lagrangian are where B µν and Z µν are the field strength tensors of U Y (1) and U (1) FN , respectively, with B µ , Z µ the corresponding gauge bosons, while W a µ are the SU (2) L gauge bosons. Since the FN fermions are charged both under the SM gauge group and the U (1) FN , the fermionic kinetic terms give the couplings of the B µ , W a µ and Z µ to the chains of FN fermions. These result in the last three terms in (5.1), where where g is the SU (2) L coupling constant, and T a the corresponding generators (now acting inside the quark currents). The ellipses denote the couplings to leptonic, up-and down-quark FN chains. These are obtained from the terms explicitly shown in (5.2), (5.3) by making replacements q L → L L , e R , u R , d R and q R → L R , e L , u L , d L , respectively. Here g Y and g are the U Y (1) and U (1) FN gauge couplings, while Y f and nδ f are, respectively, the hypercharge and the horizontal U (1) FN quantum numbers of the fermion f The couplings to the SM fermions are obtained by performing the unitary transformations in (3.4), and keeping only the zero modes. This gives where, as before, the ellipses denote the couplings to right-handed quarks and to leptons, obtained through trivial replacements. The two electroweak currents, J µ Y and J W a , are flavor diagonal, since the FN chains carry the same SM charges as the corresponding SM fermions. In contrast, the FN current, J µ FN , has flavor violating couplings, and similarly for the other c f ij matrices. The off-diagonal entries arise because the horizontal charges of fermions on different nodes differ. Note that the zero modes are localized toward the ends of vector-like chains, and thus we expect c q L ij to have eigenvalues that are O(N q(i) ), see also the discussion below (2.8).
We can get rid of the kinetic mixing between B µ and Z µ , Eq. (5.1), by performing a field redifinition After the field redefinition the covariant derivative acting on the Higgs also contains Z , After FN symmetry and electroweak symmetry , the scalar kinetic terms, . This mixing can be thought of as occurring in two steps. In the → 0 limit the electroweak breaking mixes B µ , W 3 µ into a massless photon, A µ , and the massive Z µ . For nonzero the Z µ and Z µ further mix into two mass eigenstates, with s W = sin θ W the sine of the weak mixing angle. The mass ofẐ is the same as for the SM Z, m Z , up to O( 2 ) corrections, while theẐ has the mass Note that for m Z g φ we have θ → −s W and m Z → √ 2g φ . We can finally write down the couplings of Z to the SM fermions, where the sum runs over all the SM fermions, f = u, d, , ν, with c ij ν R = 0, and i, j = 1, . . . , 3, the generation indices. The couplings receive two contributions, the flavor diagonal one from J µ Y and J µ W 3 , while the contribution from the horizontal current, J µ FN , also contains the flavor violating couplings, 11) and similarly for d L,R with u → d replacements, for L,R with q → L, u → replacements, and for ν L with q → L, u L → ν L replacements in the above expressions. For later convenience we also introduce the vector and axial couplings The flavor diagonalˆ f term is of O( ), where T f 3 is the weak isospin for fermion f . Note that the second term vanishes in the m Z → ∞ limit, cf. Eq. (5.8), while the first term is the contribution from the kinetic mixing between Z and the photon. The unitary matrices V f in (5.11) diagonalize the corresponding SM fermion Yukawa matrices (or in the case of neutrinos the Weinberg operator mass term).
There are two distinct regimes for the couplings of Z to the SM fermions. If the c f are dominated byˆ f , then the phenomenology of Z is the same as for the dark photon. In the opposite limit, whenˆ f is negligible, the couplings of the Z are governed by the U (1) FN charges, giving both flavor diagonal and off-diagonal couplings of comparable strength. In the numerical examples below we set → 0, see Table 2 and Appendix B. In this limit the Z mostly couples through axial vector couplings, since left-handed and right-handed zero modes carry opposite effective U (1) FN charges, cf. Appendix B.
The above derivations generalize straightforwardly to the case of where J µ i,F N are the currents corresponding to the horizontal U (1) i symmetries for the i-th generation, with g i the gauge coupling of U (1) i . The last term in Eq. (5.14) contains kinetic mixings between different Z i . These enter the interactions with the SM fermions only at order where the coefficients have the same general form as in Eq. (5.11), but now for each Z i separately, i.e., The V u L and V u R diagonalize the up quark mass matrix, while the c q L k matrix has only one nonzero entry, and similarly for the other c f k matrices. Since the zero modes are localized toward the ends of vector-like chains we expect c q L k kk ∼ N q(k) . Theˆ f,k coefficients are of O( i ) and vanish in the limit i → 0. Since this is the limit we will work we do not display them explicitly. For later convenience we also define the vector and axial couplings as -2.7 1.9 -2.9 2.9 -1.9 3 -1.9 3 -0.001 0.026 -0.002 2.5 -1.9 2.1 -1.9 The bounds on flavorful Z s come from a variety of experimental observables. They can be grouped into four broad categories: the bounds that come from Z couplings to quarks or from Z couplings to leptons, in each case either due to flavor diagonal or from flavor violating couplings. In the rest of this section we work out the relevant bounds for two representative benchmarks, one for the coupled and one for the uncoupled FN chains. The numerical inputs as well as the resulting Z couplings for the two benchmarks are given in Appendix B. In both benchmarks we take the → 0 limit, therefore the Z couplings are completely dictated by J µ FN , Eq. (5.4). The various experimental bounds in the (m Z , g ) plane are compiled in Figs. [13][14][15][16][17][18]. Note that the couplings of Z are mostly axial, cf. Table  2. For other phenomenological analysis of light axial vectors, see [34,35].

The bounds on flavorful
We first derive the experimental bounds on the G FN = U (1) FN benchmark, see Appendix B.2. As we will show below, the K −K mixing bounds the flavon vev to be very large, φ 10 7 GeV. This means that the FN fermions are very heavy, with masses m F ∼ O( φ ). They only give suppressed contributions to flavor observables at one loop level and can be safely ignored in our analysis. To simplify the discussion, we also assume that the flavon is very massive m φ ∼ O( φ ), giving ∼ m d i m d j /m 2 φ suppressed contributions to the flavor observables compared to the Z and can thus be ignored. This assumption can be relaxed in the future, since m φ is a free parameter. To make the analysis tractable, we also limit the mass of the Z to be above, m Z > 10 MeV, an assumption that could also be relaxed in future studies. The results obtained for G FN = U (1) FN are straightforward to extend to the U (1) 3 FN model, which we do in Section 5.2.

Flavor diagonal couplings to quarks and/or leptons
We start by deriving constraints on Z couplings from flavor conserving processes. Representative Feynman diagrams are shown in Fig. 15 . recast in Ref. [36] for a generic Z with flavor diagonal vector couplings to quarks and leptons, or to new invisible states, with the results available in the form of a public code, Darkcast.
We use Darkcast to obtain the limits on Z from direct production, adapting it to the case in hand. The U (1) FN Z has predominantly axial vector couplings, cf. Table 2, unlike dark photon, which only has flavor diagonal vector couplings. For e + e − → Z γ searches at BaBar [37,38], KLOE [39] and LEP [40,41] in Z → e + e − , µ + µ − ,inv channels, as well as for the Z bremsstrahlung [42][43][44] in electron beam dump searches at A1 [45], APEX [46], E137 [47], E141 [48], E774 [49], Orsay [50], KEK [51], NA64 [52], and in proton beam dump search at ν-CAL I [53], with Z → e + e − , the Darkcast recast of bounds applies to our Z benchmark with the replacement eQ f → ((c 11 f V ) 2 + (c 11 f A ) 2 ) 1/2 for f = e, u, d, taking into account the change in the branching ratios due to possible decays to neutrinos and heavier charged fermions. The induced uncertainties due to this identification are of roughly the same size as the uncertainties due to the approximations done in the original recast of the experiments by Ref. [36]. In the same way, the LHCb searches for dark photon [54] can be recast for m Z > 1 GeV region, where the production is dominated by Drell-Yan production, q i q i → Z . We can also safely neglect off-diagonal couplings in direct Z production, which only lead to highly suppressed corrections. The resulting bounds are shown as red (from beam dumps) and orange (from e + e − colliders) excluded regions in Fig. 13 and, assuming only couplings to leptons, in the first three panels in Fig. 14. In addition, there are a number of searches for light new particles that are harder to recast for axial Z ; by LHCb [54] for m Z < 1 GeV, by NA60 [55], CHARM [56], ν-CAL I [57] and KLOE [58]. For instance, for dark photon with mass below 1 GeV the production in pp collision is dominated by production from π 0 → γ d γ, η → γ d γ, ω → γ d π 0 decays, which can be well estimated using vector meson dominance. For axial vector it is not clear what is the dominant production channel, and would require a dedicated phenomenological analysis (for π 0 → Z γ see [34]). Using NDA we expect that the exclusions from these remaining experiments are likely to fall within or close to the red and orange exclusion regions in Fig. 13, the same as they do for the dark photon. is the average vector coupling of Z to the nucleon. In the last equality we evaluated the average for the Cs nucleus, which has Z = 55, A = 133, N = A − Z = 78, since the most stringent bounds on NP contributions to APV come from measurements of the 6s − 7s transition in Cs [59]. Translating the results of Ref.
[60] to our notation gives g 2 c 11 A c V N > 3.9 · 10 −8 m 2 Z /GeV 2 , which is not very stringent and is comparable to the nEDM bound in Fig. 14 (bottom right).
Constraints from SN1987a: The Z bosons can be copiously produced in a core of a Supernova (SN) if they are light enough. If g is small enough the Z bosons can escape the SN core and contribute to the cooling of the proto-neutron star. Demanding that this cooling mechanism does not lead to an instantaneous energy flux that is bigger than the one from neutrinos when the SN core reaches peak density, gives the bound shown as the gray region in Fig. 13. For smaller g the Z bosons are not produced efficiently, while for larger g the Z are efficiently trapped inside the SN. We restrict the analysis to m Z ≥ 10 MeV, in which case two simplifications occur that make it easy to rescale reliably the results for dark photon from Ref. [61] to our case of a Z with predominantly axial vector couplings. First of all, for m Z ≥ 10 MeV we can neglect the corrections due to the coupling of Z with the electron plasma, an important effect for lighter gauge bosons. Furthermore, for m Z ≥ 10 MeV the main mechanism of production and absorption of Z in the SN core is bremsstrahlung in neutron-proton scattering. The bounds can then be obtained by simply rescaling the results from [61]. We use that the pn → pnZ cross sections, σ A,V , induced by the axial or vector couplings of Z to proton, respectively, satisfy the numerical relation σ A /σ V ∼ 3c 2 pA /c 2 pV , where c pA,V = 2c 11 uA,V + c 11 dA,V . We can thus translate the bounds in [61] by replacing with √ 3g c pA , giving the gray regions in Figs. 13 and 14.
Neutrino trident production: Neutrino scattering on nucleus, A, can produce lepton pairs through electroweak interactions. Such a trident process, ν i A → ν j + k − l A, can receive a contribution from a tree level exchange of an extra Z , see Fig. 15 (middle). This would result in a deviation of the total cross section from the SM prediction, signaling new physics [62]. For small Z couplings the main correction to the trident cross section comes from the interference with the leading SM contribution, which is due to the Z exchange. This means that the results of Ref. [63] for the L µ − L τ gauge boson directly translate to our case (see also [64]). The bound on g from Fig. 8 of [63] needs only to be re-interpreted as the bound on g c 22 of the weak mixing angle. The resulting bounds from the CCFR experiment [65] and the future projections for DUNE [63] are shown in Fig. 14 (top right) as solid black and dashed black lines, respectively.
Electron-neutrino scattering experiments: Measurements of ν i − e − andν i − e − scattering cross sections can bound the couplings c 11 A c ii νL and c 11 V c ii νL , where i is the flavor of the neutrino in the beam. The most stringent constraints are due to CHARM-II [66], which used O(10 GeV) ν µ andν µ beams, and from TEXONO [67], which used O(1 MeV) reactor ν e beam. We use the CHARM-II measured ratio σ(ν µ e)/σ(ν µ e) = 0.910 ± 0.113 [66], where σ(ν µ e) and σ(ν µ e) are the total ν µ e andν µ e scattering cross sections, respectively. Using this result we set the bound shown as a cyan line in Fig. 14 (top right). The TEXONO measured totalν e scattering rate, normalized to the SM prediction, R exp /R SM = 1.08±0.36 [67] results in the bound on g shown as a green solid line in Fig. 14 (top left). To derive these constraints we used the analytical results from Ref. [68].
Møller scattering: Measurements of parity violating contributions to the Møller scattering, e − e − → e − e − , bound the product c 11 V c 11 A . The most precise measurements were performed by E158 at SLAC at q 2 (0.16 GeV) 2 [69]. Comparison with the SM gives for m Z 100 MeV the bound g 2 |c 11 V c 11 A | 10 −8 , and for m Z 100 MeV the bound g 2 |c 11 V c 11 A | 3 · 10 −8 · (m Z /200 MeV) 2 (see also [34]). In our benchmark c 11 V c 11 A 0.30, giving relatively weak bounds g 2 · 10 −4 and g 2 · 10 −4 (m Z /100 MeV), respectively, which we thus do not plot in Fig. 14.
Isotope shift spectroscopy: The isotope shift spectroscopy constrains vector couplings of Z , giving g 2 |c 11 V c 11 f V | 10 −7 · m 2 Z /(10 MeV) 2 for f = u, d, when charge radius determination from Lamb shift in muonic atoms is used [70] (see also [71][72][73]). Since the Z has suppressed vector couplings, these bounds are not very constraining, and we do not consider them further.
White dwarf cooling: The tree level Z exchange contributes to the e + e − → νν process which can increase the cooling rate of the white dwarf (WD) core [74]. The cooling rate due to this additional cooling mechanism should not exceed the SM cooling rate due to the plasmon decaying into neutrinos. For our Z benchmarks the contribution to the star cooling is described by the effective Lagrangian since the Z is much heavier than the WD internal temperature of a few keV, and can be integrated out. Translating the limits from [74] to our notation gives (see also [75]), In our benchmark, c eff νL 3.4, and c eff e 2.5, giving the exclusion shown in Fig. 14 (top left) as a purple region.
Anomalous magnetic moments: At one loop the Z exchange contributes to the lepton anomalous magnetic moment, (g − 2) , cf. Fig. 17 (with i = j ). For our Z benchmark we only need to keep the contribution to (g − 2) µ from the diagonal couplings, with µ running in the loop. A similar diagram with a τ running in the loop does get a chirality flip enhancement of m τ /m µ ∼ 10, but is also suppressed by two off-diagonal couplings, For the (g − 2) e the two corresponding factors are m τ /m e ∼ 3 · 10 3 and |(c 13 L c 13 R )/(c 11 L ) 2 | ∼ 6 · 10 −5 , so that again the diagram with diagonal couplings dominates. Using the results of Ref. [76] with a trivial change of notation gives (see also [35]), with i = 1, 2 for = e, µ, respectively. The normalization of the loop functions, is such that for heavy Z , lim u→∞ F (u) = 1, lim u→∞ H(u) = 5/3, whereas for m Z = m , F (1) 0.31, H(1) 1.31, and for light Z , i.e., for u 1, F (u) 3u 2 /2, H(u) 1. In Fig. 14 (top right) we show in gray the 1σ band in the g , m Z parameter space that gives (g − 2) µ in agreement with experiment. The required value of g is excluded by a number of other measurements, among others also by the limit on the allowed NP contribution to (g − 2) e , denoted with a blue line in Fig. 14 (top left).
Note that here we consider only the contributions to (g − 2) f from Z running in the loop. The contributions from heavy vector-like fermions running in the loop are expected to give contributions that are parametrically of the same order. We do not attempt to include these contributions given that the (g − 2) f measurements put only weak bounds on our inverted FN benchmark. However, because of this approximation one should view the (g − 2) f bounds shown in Fig. 14 only as indicative. Similar comments apply to electron EDM, which we discuss next.
Electric dipole moments: The complex off-diagonal couplings of Z generate at one loop the electric dipole moments (EDMs). The contribution to the EDM of fermion f i from the fermion f j running in the loop is given by [77]  where C + 1 , C 0 , are the three-point one loop integrals arising from the evaluation of the Feynman diagram in Fig. 17 (right), the analytic expression for which can be found in [78]. The chromo-EDMs,d f i , are obtained from (5.24) by replacing g → g s . The experimental bound on electron EDM, |d e | < 1.1 · 10 −29 e cm [79], results in a bound shown as a black line in Fig. 14 (top left). The bound on neutron EDM, d n , on the other hand, translates to a bound denoted with a blue line in Fig. 14 (bottom right). The neutron EDM receives contributions both from quark EDMs, d q , and quark chromo-EDMs,d q , so that d n = q=u,d,s β q d q +β qdq . For matrix elements of quark EDM operators, β q , a recent lattice QCD calculation obtained β u = 0.784 (30)

Flavor violating couplings to quarks
The exchange of Z induces Flavor Changing Neutral Currents (FCNCs) already at tree level, as shown in Fig. 16. Strong constraints are obtained from the meson mixings, which we derive in detail in the next paragraphs.
K 0 −K 0 mixing: To estimate the contributions of tree level Z exchanges to K −K mixing we distinguish two regimes. For light Z , m Z m K , we can use ChPT with the Z acting as an external field. Using ChPT expansion we have (sγ µ γ 5 d) → −f K ∂ µ K 0 + · · · , where the ellipses denote higher orders in momentum expansion. Using this in Eq. (5.10) gives for the tree level Z exchange contribution to K −K mixing (cf. Eq. (5.12)) where f K 156 MeV is the kaon decay constant, and we use the same phase conventions for M 12 as in [82,83]. Note that this contribution is proportional to 1/m 2 Z even though the momentum flowing in the Z propagator is O(m K ), so that M Z 12 ∝ 1/ φ 2 . This is most easily understood in the Feynman gauge, where the dominant contribution to M Z 12 due to the exchange of the arg(φ) Goldstone boson, whose couplings to SM fermions are ∝ 1/ φ and do not depend on g . For m Z m K we can integrate out the Z at the scale µ m Z and match onto the effective weak Hamiltonian, H eff = a C a Q a , where the sum runs over the operators in Ref. [84]. The nonzero NP contributions to the Wilson coefficients are, To compare with the experimental results we run from µ m Z down to 2 GeV using the results of [84,85].
The comparison of the SM predictions and the experimental measurements for indirect CP violation parameter K is made through a ratio Estimating the errors due to charm loop contributions in the SM prediction of K is nontrivial and leads to differences between CKMfitter and UTFit collaborations [86,87] (see also the introduction in [88]). For consistency we use the UTFit collaboration extraction 0.87 < C K < 1.39 at 95% C.L.
[87] as well as their prediction for SM K using NLO charm contribution [89]. We denote the resulting bounds in Fig. 13 and 14 (bottom right) with black lines, where we use the heavy Z solution (5.26) for m Z > 3m K , the light Z solution (5.25) for m Z < m K /3, and connect the two with naive linear interpolation to guide the eye. The Z contribution to the K −K mixing is ∝ 1/ φ 2 . For heavy Z , m Z m K , this then leads to a bound φ = m Z /g > 1.2 · 10 7 GeV for our benchmark.
B 0 q −B 0 q mixing: We distinguish two limits, the heavy and light Z . For heavy Z , m Z m Bq , the Z is integrated out at µ m Z . The matching onto H eff = a C a Q a , gives, analogously to (5.26), where q = d, s, with i = 1, 2, respectively. We RG evolve the above Wilson coefficients from µ m Z down to µ = 4.2 GeV m b using the results of [84,85], and compare with the allowed deviations in the mixing matrix elements, The latest fit results from UTFit collaboration give 0.942 < C Bs < 1.288, −1.35 • < φ Bs < 2.21 • , and 0.83 Center : decay to three leptons, i → j j¯ j . Right: radiative decay, i → j γ. by with m b 4.2 GeV the b quark mass, while the remaining nonzero Wilson coefficients, C qb , are parametrically suppressed. We denote the resulting bound in Fig. 14

Flavor violating couplings to leptons
In this subsection we derive the experimental bounds from lepton flavor violating transitions. The relevant diagrams are summarized in Fig. 17.
µ → e conversion: Stopped muons can undergo coherent µ → e conversion in the field of the nucleus. The most stringent limits to date were obtained by the SINDRUM II experiment using 197 Au [91] Here Γ(µ − Au → e − Au) is the conversion rate, and Γ capt (µ − Au) the muon capture rate. State of the art estimates for both of these can be found in Ref. [92] for a whole range of nuclei, including gold. The tree level Z exchange induces four-fermion effective interaction for µ → e conversion, see Fig. 17, which in the notation of Refs. [92,93] takes the form 32) and similarly for g RV (f ) , but with L → R replacement, while g LA(f ) and g RA(f ) trivially follow from V → A replacement.
The g RV (f ) and g LV (f ) operators induce coherent µ → e conversion for which we use the results of Ref. [92]. The g RA(f ) and g LA(f ) operators induce spin dependent µ → e conversion, for which we use the results of [94], where the Wilson coefficients in the notation of [94] are C f f A,L = g LA(f ) /4, C f f A,L = g LA(f ) /4, while the axial structure factors S A we take from [95]. Note that the momentum exchange in µ → e conversion is q 2 −m 2 µ . This is small enough that it is easily absorbed by the nucleus without changing its nuclear structure. On the other hand, q 2 −m 2 µ is also large enough that even for very light Z , m Z m µ , the interaction can still be viewed as point-like in the calculation of conversion rate, at least for the heavy nuclei. The neutron and proton density distributions of heavy nuclei have radial extent which is parametrically bigger than the muon Compton wavelengths, for gold roughly r max ∼ 4/m µ . This means that we can still use the calculations of overlap integrals from Ref. [92] even for light Z , up to corrections of about O(30%).
Spin-dependent conversion starts to be dominant at very low masses. While the spinindependent conversion rate is coherently enhanced and at low masses goes as ∼ A 2 m µ , where A is the atomic number, the spin-dependent rate increases as m 3 µ /m 2 Z . For heavy nuclei, such as 197 Au used by the SINDRUM II experiment, the spin-dependent contribution starts to dominate around m Z ∼ 0.5 MeV, while for lighter nuclei, such as 26 Al used by Mu2e, it begins at m Z ∼ 10 MeV. In both cases, the spin-dependent rate is important for masses outside our range of interest, so we only consider the spin-independent rate in obtaining the bounds on the Z parameter space.
The constraint from Eq. (5.31) is denoted in Fig. 13 with a solid green line, while in Fig. 14 (top right) we also show the comparison with the other flavor violating processes involving muons. The Mu2e collaboration plans to achieve the sensitivity Br(µ → e) < 8 · 10 −17 for µ → e conversion on Al nuclei [96], which gives the projected sensitivity on g denoted with a dashed green line in Figs. 13 and 14 (top right).
Depending on the value of m Z there are three different regimes for Γ( i → 3 j ). For m Z m i , the Z can be integrated out, leading to effective four fermion interaction. In this limit, the i → 3 j decay is a genuine three body decay, with the decay width [97] (5.33) where we neglected the terms proportional to m j .
In the intermediate mass regime, 2m j < m Z < m i − m j , the Z can be produced on shell. The i → 3 j decay is thus a cascade decay, i → Z j followed by Z → + j − j , so that where we neglected the mass of j , while Br(Z → j j ) is the branching ratio for the Z → + j − j decay. In our benchmark Br(Z → µµ) ∼ 0.16, while Br(Z → ee) ∼ 0.23 for m µ < m Z < m τ and Br(Z → ee) ∼ 0.28 for m Z < m µ . Note that Γ( i → j Z ) ∝ g 2 /m 2 Z = 1/ φ 2 and thus does not diverge for m Z → 0. Finally, if m Z < m j the Z is off-shell and we again have a genuine three body decay. For m Z m j the axial current contributions are 1/m 4 Z enhanced, giving [97] Γ( i → 3 j ) The contributions from flavor diagonal vector couplings of Z are subleading and we can safely neglect them.

Rare meson decays (RMD):
The tree level Z exchanges induce rare FCNC meson decays such as is the initial (final) state meson. An example of such an exotic decay is K + → π + µ + e − which arises at one loop in the SM but with a tremendous suppression, with a matrix element proportional to the neutrino masses, and are thus vanishingly small in practice. The tree level transition matrix element due to a Z exchange instead has no suppression due to neutrino masses. For heavy Z the off-shell contribution is proportional to a product of two off-diagonal couplings, one for the quark flavor changing current and one for the lepton current. The resulting bounds are therefore expected to be weaker than the bounds that follow from the lepton FCNC transitions or from meson mixings. In contrast, when the Z is light enough to be produced on-shell, the M 1 → M 2 Z decays depend only on quark current off-diagonal couplings. In this case the bounds from rare meson decays in general become the most stringent ones.
For the M 1 → M 2 transitions the strongest experimental bounds follow from the M 1 → M 2 Z decays with M 1,2 both pseudoscalars, for which the matrix element is with J µ defined in (5.10), and p 1,2 the momenta of mesons M 1,2 , such that p 2 1,2 = m 2 1,2 , with m 1,2 the corresponding meson masses. The quark flavor-violating coupling c kl qV encodes the strength of the q l → q k transition, e.g., the c 23 dV coupling is obtained for the B → K and c 13 dV for the B → π matrix elements, respectively. The form factors f +,− (q 2 ) depend on the momentum exchange squared, q 2 , with q µ = (p 1 − p 2 ) µ . We follow the common practice and trade the form factor f − (q 2 ) for We use the z-expansion based parametrization of the form factors [102], For light Z , with a mass below the lepton pair production threshold, m Z ≤ m i + m j , the M 1 → M 2 − i + j transition is a genuine 3-body decay, mediated by an off-shell Z . The longitudinal Z contribution dominates, leading to a decay width proportional to ∝ g 4 /m 4 Z 1/(4 φ 4 ). The expression for the differential width is well approximated by taking the m Z → 0 limit, giving where m 2 i j = (p i + p j ) 2 is the invariant mass squared of the leptonic pair, and similarly m 2 For heavy Z , with a mass m Z m 1 , the Z can be integrated out. This also gives a decay width proportional to ∝ g 4 /m 4 Z 1/(4 φ 4 ). In the limit of massless final states, m 2 , m i , m j → 0, the differential rate is given by a rather compact expression, (5.42) In the numerical analysis we work with the full dependence on m 2 , m i , m j .
Phenomenologically, the most relevant is the intermediate mass regime, m i + m j < m Z < m 1 − m 2 . In this case the Z is produced on-shell in the 2-body decay, M 1 → M 2 Z , followed by a decay into two leptons, Z → i¯ j . The total decay width is given by The two-body decay width for M 1 → M 2 Z transition depends only on the quark current matrix element, and is given by where y i = m i /m 1 . The branching fractions for Z → i¯ j are easily calculated. Ignoring the kinematical factors for final state particles they are given by Br(Z → i¯ j ) = (|c ij . For our benchmark we have Br(Z → ee, µµ) 0.1 and Br(Z → α ν α ν α ) 0.4 for the flavor conserving channels, and Br(Z → eµ) ∼ 10 −5 for the flavor violating decay.
In the numerical analysis we use the PDG values for the measurements of the branching ratios and the upper bounds [103], while the numerical inputs for the form factors are taken from Ref. [104]. The envelope of constraints on g as a function of Z mass due to rare meson decays (RMD) is shown as a blue solid line in Figs. 13 and 14. We see that the RMD constraints are the most stringent ones up to several GeV Z masses, i.e., as long as the B → πZ decays are still kinematically allowed.
The kinematical thresholds from K, D, B decays are clearly visible in Figs. 13 and 14. At the lowest Z masses the K + → π + Z , Z → νν decay leads to the strongest constraints, a result of K → π + inv being well constrained experimentally, as well as the large Z → νν branching ratio. Conservatively, we use the integrated branching ratio Br(K + → π + νν) < 2.3 × 10 −10 to set the limit on g in Figs. 13 and 14. A more careful analysis of missing mass differential rates, taking care of spill-overs between different bins, would further improve these bounds.
Flavor violating radiative decays: At one loop the off diagonal couplings of Z to leptons generate the i → j γ transition with the decay width [105], see also Fig. 17, where x a/b = m a /m b , and with Q k the charge of the fermion k running in the loop, while c γ R can be obtained by making the replacement L ↔ R in the above expression. The loop functions are where (below we shorten x k/Z → x) The present bounds on Br(τ → µ(e)γ) < 4.4(3.3) · 10 −8 [106] and Br(µ → eγ) < 4.2 · 10 −13 [107] are denoted in Fig. 14 (bottom left) by a brown(black) solid line and in Fig. 14 (top right) by a blue solid line. The projected sensitivities at Belle II [101] and MEG-II [108] Br(τ → µ(e)γ) < 1(3) · 10 −9 , Br(µ → eγ) < 6 · 10 −14 , are shown with the corresponding dashed lines. Note that in our benchmark for all three transitions the largest contribution comes from a diagram with a τ running in the loop. Note also, that the contributions from vector-like fermions, which we do not take into account, are parametrically of the same order. The resulting bounds in Fig. 14 should thus be taken only as indicative. For all three Z i the most stringent constraints arise from rare meson decays, µ → e conversion, K −K mixing and SN constraints. This is the same as for the U (1) FN , Section 5.1. However, the relative importance of these constrains has shifted because of the changes in the flavor patterns of the Z i couplings, see Fig. 18.

Bounds on the
In Fig. 18 we also show the constraints from beam-dump (red regions) and from e + e − experiments (orange regions). Here some care is needed, since for Z 2 and Z 3 the 12 off-diagonal couplings can dominate over the diagonal 11 couplings, see, e.g., Eq. (B.37). This means that we need to adjust our treatment of bounds from direct production of Z i compared to what was done in Section 5.1. For m e + m µ < m Z < 2m µ the Z i → eµ decay is now the dominant decay channel. In the G FN = U (1) 3 FN benchmark we have Γ(Z 2 → ee)/Γ(Z 2 → eµ) 0.16 and Γ(Z 3 → ee)/Γ(Z 3 → eµ) 0.52, which is then used to rescale appropriately the bounds from Darkcast, cf. Section 5.1.1. The beam dump constraints are stronger for Z 1 than they are for Z 2 and Z 3 , because Z 2,3 have suppressed couplings to the first generation quarks. This is also the reason that the regions excluded . by SN constraints (gray regions in Fig. 18), are shifted to larger values of gauge coupling g for Z 2 and Z 3 . The expressions for these can be taken directly from Section 5.1, with trivial relabeling corresponding to Z → Z i replacements. For Z 1 the most stringent constraints come from the FCNCs involving the first two generations, rare kaon decays, µ → e conversion and K −K mixing, since the Z 1 has appreciable 12 offdiagonal couplings. The D −D mixing, even though less stringent, can still be an important constraint, especially for future projections, while the FCNCs involving third generations lead to weaker bounds. For illustration we show the most stringent one, from B d −B d mixing. The Z 2 and Z 3 have much smaller 12 off-diagonal couplings. Even so, rare meson decays, µ → e conversion and K −K mixing still lead to the most stringent constraints in our benchmark, but the FCNCs involving the third generation, B d −B d and B s −B s mixing, are relatively more important.
In discussing the above constraints it is important to keep in mind that we show bounds for a single benchmark. We expect other G FN = U (1) 3 FN benchmarks to lead to a similar pattern of constraints. However, the relative strengths of different constraints may easily differ by relative O(1) factors (for instance for some benchmarks we would expect µ → e conversion to be more stringent than K −K mixing even for heavy Z 1 , unlike what was found here). Fig. 18 should thus be taken only as an illustration of how important different probes are for G FN = U (1) 3 FN inverted FN models.

Conclusions
We introduced a class of anomaly-free Froggatt-Nielsen (FN) models that can simultaneously explain the fermion mass hierarchy as well as the mixing patterns. These inverted FN models differ from the traditional FN models in that the expansion parameter has the inverted form, M/ φ , where M is the vector-like mass parameter and φ the vev of the flavon field. The observed pattern of masses and mixings is obtained for φ M , while in the traditional FN models the opposite is required, M φ . Different realizations of the inverted FN models differ in the choice of the FN horizontal group G FN and the assignment of FN charges. The common feature, on the other hand, is that the fields that are chiral under the SM do not carry G FN charges. These are the fields that couple directly to the Higgs, and also to the chains of vector-like fermions that are charged under G FN . This set-up gives zero modes that are localized toward the ends of chains, with exponentially suppressed overlaps with the zero node where the Higgs resides. After electroweak symmetry breaking this results in the hierarchy of the SM fermion masses. The set-up also ensures that G FN is anomaly free irrespective of the choice of the FN charges, making it easy to extend our work to other charge assignments or gauge groups.
In the paper we explored in detail two choices for G FN : the "decoupled chains" model, where G FN = U (1) 3 FN , and the "coupled chains" model, with G FN = U (1) FN . For both cases we showed using a numerical scan of input parameters that one can obtain the observed hierarchy of SM quark masses and CKM matrix elements given the appropriate choice of FN charges. The scan was over random O(1) complex flavon and Higgs Yukawa couplings and over random vector-like mass parameters, M a , taken to be on average a factor 5 smaller than φ . The obtained solutions for the pattern of quark masses and mixings do not contain any tunings -the associated Barbieri-Giudice measure only reaches values of a few. Furthermore, the solutions simultaneously give the hierarchical quark masses and hierarchical mixing angles. Demanding that the quark masses are the measured ones, the probability for the scan to give all three mixing angles below their SM values is about 8% in both FN models (to be compared with the best case scenario of 0.5 3 = 12.5%, and much higher than 3 · 10 −6 for mixing angles that are distributed completely randomly). Similarly, we showed that the neutrino mass differences and PMNS matrix elements are well described by a completely anarchic form of a Weinberg dimension 5 operator. Since the Weinberg operator couples to chains of vector-like fermions the neutrino masses are still hierarchical, preferentially resulting in a normal ordering.
The inverted FN models do share a problem common to all FN models. Since one introduces a relatively large number of states charged under QCD, the theory is no longer asymptotically free in the UV, see, e.g., [27]. In the two models we presented there is a QCD Landau pole about four orders of magnitude above φ , i.e., for φ ∼ 10 7 GeV this would be around 10 11 GeV. One could potentially push this to even higher scales by adjusting the lengths of FN chains along with the size of φ /M expansion parameter. Alternatively, for φ ∼ 10 14 − 10 15 GeV, depending on the exact field content, the Landau pole is above the Planck scale. Another interesting possibility is that with an appropriate field content one may realize an asymptotically safe theory [109]. The possible Landau pole in the U (1) FN is less constraining; for g < 0.11 at µ 10 7 GeV the Landau pole in g is above the Planck scale.
In the second part of the paper we explored the phenomenological implications of the two inverted FN models. The structures we introduced leave imprints in FCNC transitions, which then leads to stringent bounds on the allowed values of φ and M a . We focused on the case where G FN is gauged, resulting in three (one) gauge bosons Z i (Z ) for G FN = U (1) There are several ways in which our study could be extended in the future. Most immediately, one could explore the bounds on inverted FN models, which are not gauged so that there are no Z tree level contributions to the FCNCs. It would also be interesting to explore whether or not the inverted FN models can aid in exploring the experimental anomalies in b → s + − transitions, potentially along the lines of Ref. [110].
In order to construct the appropriate ChPT Lagrangian in the presence of flavor violating Z we use the spurion analysis [111]. The QCD Lagrangian is now L QCD+Z = q(i / ∂ + g s / G a T a + χ V / Z + χ A / Z γ 5 )q −qM q q, where q = (u, d, s), M q = diag(m u , m d , m s ), and χ V,A are 3 × 3 Hermitian matrices of the form and that take the values v µ = χ V Z µ , a µ = χ A Z µ , s = M q , p = 0. The LO ChPT Lagrangian, including Z as the light degree of freedom, is therefore (we work in the unitary gauge for the Z gauge boson) where U = exp(iλ a π a /f ), with λ a the Gell-Mann matrices, f is related to the meson decay constant Since we consider radiative corrections due to the Z gauge bosons, the Z cannot be treated as merely an external field that enters into the spurions. In constructing the ChPT Lagrangian this means that there are two additional hermitian spurions, Q R,L = χ V ± χ A , that transform as Q R,L → g R,L (x)Q R,L g † R,L (x). This gives the following contributions to the LO ChPT Lagrangian, where we only show the terms that will be relevant for K −K mixing, where they will serve as counter-terms to one-loop corrections. For this we also need part of the O(p 4 ) Lagrangian, (A.4) 5 We use the normalization 0|q1γµq2|P (p) = ipµfP .

The ChPT Lagrangian is invariant under parity
In addition, since the terms containing Q R,L can only arise from loops of Z , the ChPT terms without external Z fields need to have even powers of g . This is insured by requiring that the ChPT Lagrangian is invariant under the Z 2 transformation Q L,R → −Q L,R . The chiral counting of the spurions is Q R,L ∼ O(p 0 ), v µ , a µ , ∼ O(p), s + ip ∼ O(p 2 ). In principle one could thus have an arbitrary number of Q R,L insertions in L (2),(4) Q R,L . However, we work only to O(g 2 ) and thus merely keep the terms that are ∼ Q 2 R,L . Note that the first term in (A.3), the analogue of the ChPT+QED Lagrangian [112,113], is in our counting O(p 0 ). However, its coefficient is C where we abbreviated m uds = m u + m d + m s , and m ds = m d + m s . Expanding the NLO Lagrangian (A.2) in meson fields gives (A.6) Above, we kept only the terms relevant for K 0 −K 0 mixing, defined φ 1

and replaced
√ 2f with the kaon decay constant, f K = 155.6 ± 0.4 MeV, to account for the SU(3) breaking. The first term in (A.5) is due to the axial vector coupling of the Z boson and leads to the tree level contribution to the K 0 −K 0 mixing amplitude, Z ′ Figure 19. The one loop diagrams contributing to K −K mixing.
. Eq. (5.25). The second term is due to the vector coupling and contributes at one loop, see left diagram in Fig. 19. The first and the last terms in (A.5) contribute to the middle diagram in Fig. 19, while the Z µ Z µ term in the second line of (A.5) gives the last diagram in Fig. 19.
The one-loop contributions to K −K mixing are with c π 0 = −1/ √ 2, c η = 3/2, and A 0 and B 0 the Passarino-Veltman loop functions, see, e.g., [114]. The counter-terms give (A.8) Since there are no one-loop contributions of the form ∼ c 12 d A + c 12 d V , this implies a relation between the counter-tems in the last line, m 2 d C 3 .

B Benchmarks
In this appendix we give the details of the two benchmarks, one for the case of the decoupled chains with U (1) 3 FN horizontal symmetry, and one for the coupled chains with the U (1) FN horizontal symmetry.  where again the entries correspond to values on different nodes. Note that the ratios |(Y f φ ) ij /(M f ) ij | are on average equal to q = 5, but have a distribution that allows for O(1) deviations (in relative terms) from this value. The Yukawa couplings on the zero node, which couple chiral fermions to the Higgs, are The absolute values, r a , for each of the above entries were taken to be in the interval r a ∈ [0.3, 0.9], when constructing the benchmark. Finally, the matrix of coefficients in the neutrino mass term, Eq. In the benchmark we set the kinetic mixing between different Z i to zero, as we do the kinetic mixing of Z i with hypercharge. The mass eigenstates, Z i , therefore correspond to the gauge bosons coupling to the i−th fermion generation in the flavor basis, that is, before the electroweak symmetry breaking. After fermion mass diagonalization that includes the electroweak symmetry breaking terms, the real parts of the Hermitian coupling matrices in Eq.