Anomaly free Froggatt-Nielsen models of flavor

TL;DR

This paper presents two new anomaly-free Froggatt-Nielsen models with horizontal symmetries, inspired by the clockwork mechanism, which generate flavor structures and can be tested through flavor observables and potential light gauge bosons.

Contribution

The paper introduces novel anomaly-free FN models based on $U(1)^3$ or $U(1)$ symmetries, utilizing a clockwork-inspired structure with vector-like fermion chains, expanding the mechanisms for flavor hierarchy generation.

Findings

Models generate quark and lepton flavor hierarchies via a new expansion mechanism.

Current flavor bounds imply vector-like fermion masses above ~10^7 GeV.

Gauged models predict light $Z'$ bosons detectable through flavor and astrophysical experiments.

Abstract

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.

Tables

Table 1. Table 1: The experimental values of quark and charged lepton masses (in GeV) at μ = 10 7 𝜇 superscript 10 7 \mu=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 subscript 𝑚 𝑢 m_{u} and m d subscript 𝑚 𝑑 m_{d} we show only the central values, which coincide with the choice made in Ref.  Alonso:2018bcg at μ = 2 𝜇 2 \mu=2 TeV.

${\overline{m}}_u$	${\overline{m}}_d$	${\overline{m}}_c$	${\overline{m}}_s$	${\overline{m}}_t$	${\overline{m}}_b$	${\overline{m}}_e$	${\overline{m}}_{\mu }$	${\overline{m}}_{\tau }$
0.0007	0.0015	0.33	0.03	95	1.6	0.0005	0.104	1.798

Table 2. Table 2: The flavor diagonal couplings of Z ′ superscript 𝑍 ′ Z^{\prime} to the SM fermions, Eq. ( 44 ), for the G FN = U ​ ( 1 ) FN subscript 𝐺 FN 𝑈 subscript 1 FN G_{\rm FN}=U(1)_{\rm FN} benchmark. For a complete list see Appendix B .

$i$	$c_{u_L}^{ii}$	$c_{u_R}^{ii}$	$c_{d_L}^{ii}$	$c_{d_R}^{ii}$	$c_{{\ell }_L}^{ii}$	$c_{{\ell }_R}^{ii}$	$c_{{\nu }_L}^{ii}$
$1$	-2.3	3	-2	3	-2	3	-2
$2$	-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

Equations

$$\big{(}Y_{u}\big{)}_{ij}\sim\biggr{(}\frac{\langle\phi\rangle}{M}\biggr{)}^{[q_{i}]-[u_{j}]},\qquad\big{(}Y_{d}\big{)}_{ij}\sim\biggr{(}\frac{\langle\phi\rangle}{M}\biggr{)}^{[q_{i}]-[d_{j}]},$$

$$\big{(}Y_{u}\big{)}_{ij}\sim\biggr{(}\frac{M}{\langle\phi\rangle}\biggr{)}^{[q_{i}]-[u_{j}]},\qquad\big{(}Y_{d}\big{)}_{ij}\sim\biggr{(}\frac{M}{\langle\phi\rangle}\biggr{)}^{[q_{i}]-[d_{j}]}.$$

$${\cal L}_{1}={\cal L}_{q}+{\cal L}_{d}+\big{(}Y^{d}_{0}\bar{q}_{L,0}d_{R,0}H+{\rm h.c.}\big{)},$$

$${\mathcal{L}}_q$$

$${\mathcal{L}}_d$$

$$q_{L,0}^{\prime }=\sum _{n=0}^{N_q}{V}_{n0}^{q_L}{q}_{L,n},d_{R,0}^{\prime }=\sum _{n=0}^{N_d}{V}_{n0}^{d_R}{d}_{R,n}.$$

$$V_{n0}^{q_{L}}={\cal N}_{0}^{q_{L}}\biggr{(}\frac{M}{\langle\phi\rangle}\biggr{)}^{N_{q}-n},\qquad V_{n0}^{d_{R}}={\cal N}_{0}^{d_{R}}\biggr{(}\frac{M}{\langle\phi\rangle}\biggr{)}^{N_{d}-n},$$

$$m_d\simeq f_{q_L}{f}_{d_R}{Y}_0^{d}v/\sqrt{2},$$

$$f_{\psi}\equiv V_{00}^{\psi}\simeq\biggr{(}\frac{M}{\langle\phi\rangle}\biggr{)}^{N_{\psi}},$$

$${\cal L}_{\rm SM}\sim Y_{0}^{d}\biggr{(}\frac{M}{\langle\phi\rangle}\biggr{)}^{N_{q}+N_{d}}\bar{q}_{L,N_{q}}d_{R,N_{d}}H,$$

$${\cal L}_{f}=\sum_{i}\Big{(}{\cal L}_{q{(i)}}+{\cal L}_{d{(i)}}+{\cal L}_{u{(i)}}\Big{)}+\sum_{ij}\Big{[}\big{(}Y^{d}_{0}\big{)}_{ij}\bar{q}_{L,0}^{(i)}d_{R,0}^{(j)}H+\big{(}Y^{u}_{0}\big{)}_{ij}\bar{q}_{L,0}^{(i)}u_{R,0}^{(j)}\tilde{H}+{\rm h.c.}\Big{]},$$

$$(m_{d})_{ij}\simeq f_{q_{L}^{(i)}}f_{d_{R}^{(j)}}\big{(}Y_{0}^{d}\big{)}_{ij}\frac{v}{\sqrt{2}},\qquad(m_{u})_{ij}\simeq f_{q_{L}^{(i)}}f_{u_{R}^{(j)}}\big{(}Y_{0}^{u}\big{)}_{ij}\frac{v}{\sqrt{2}}.$$

$$\frac{M_1}{⟨{\varphi }_1⟩}=\frac{M_2}{⟨{\varphi }_2⟩}=\frac{M_3}{⟨{\varphi }_3⟩}=\frac{1}{q}\simeq \lambda =0.2,$$

$$(m_d{)}_{ij}\simeq {\lambda }^{N_{q(i)}+N_{d(j)}}\frac{v}{\sqrt{2}},(m_u{)}_{ij}\simeq {\lambda }^{N_{q(i)}+N_{u(j)}}\frac{v}{\sqrt{2}}.$$

$$\begin{array}{rl}{N}_{q(1)}& =3,N_{q(2)}=2,N_{q(3)}=0,\\ N_{u(1)}& =4,N_{u(2)}=2,N_{u(3)}=0,\\ N_{d(1)}& =4,N_{d(2)}=3,N_{d(3)}=3,\end{array}$$

$${\cal L}_{1}={\cal L}_{q}+{\cal L}_{d}+{\cal L}_{u}+\sum_{ij}\Big{[}\big{(}Y^{d}_{0}\big{)}_{ij}\bar{q}_{L,0}^{(i)}d_{R,0}^{(j)}H+\big{(}Y^{u}_{0}\big{)}_{ij}\bar{q}_{L,0}^{(i)}u_{R,0}^{(j)}\tilde{H}+{\rm h.c.}\Big{]},$$

$$\begin{split}{\cal L}_{d}=&i\sum_{n=0}^{N_{d(i)}}\bar{d}_{R,n}^{(i)}\not{D}d_{R,n}^{(i)}+i\sum_{n=1}^{N_{d(i)}}\bar{d}_{L,n}^{(i)}\not{D}d_{L,n}^{(i)}\\ &-\sum_{n=1}^{N_{d(1)}}\Big{(}\sum_{i,j=1}^{\hat{N}_{d}|_{n}}\big{(}M_{n}^{d})_{ij}\bar{d}_{L,n}^{(i)}d_{R,n}^{(j)}-\sum_{i=1}^{\hat{N}_{d}|_{n}}\sum_{j=1}^{\hat{N}_{d}|_{n-1}}(Y^{d}_{n})_{ij}\phi\bar{d}_{L,n}^{(i)}d_{R,n-1}^{(j)}+{\rm h.c.}\Big{)},\end{split}$$

$$\begin{split}{\cal L}_{q}=&i\sum_{n=1}^{N_{q(i)}}\bar{q}_{R,n}^{(i)}\not{D}q_{R,n}^{(i)}+i\sum_{n=0}^{N_{q(i)}}\bar{q}_{L,n}^{(i)}\not{D}q_{L,n}^{(i)}\\ &-\sum_{n=1}^{N_{q(1)}}\Big{(}\sum_{i,j=1}^{\hat{N}_{q}|_{n}}\big{(}M_{n}^{q}\big{)}_{ji}\bar{q}_{L,n}^{(j)}q_{R,n}^{(i)}-\sum_{i=1}^{\hat{N}_{q}|_{n}}\sum_{j=1}^{\hat{N}_{q}|_{n-1}}\big{(}Y^{q}_{n}\big{)}_{ji}\phi\bar{q}_{L,n-1}^{(j)}q_{R,n}^{(i)}+{\rm h.c.}\big{)}.\end{split}$$

$$q^{\prime }{}_{L,0}^{(i)}=\sum _{n=0}^{N_{q(1)}}\sum _{j=1}^{{\hat{N}}_q{|}_n}{V}_{n(j),0(i)}^{q_L}{q}_{L,n}^{(j)},d^{\prime }{}_{R,0}^{(i)}=\sum _{n=0}^{N_{d(1)}}\sum _{j=1}^{{\hat{N}}_q{|}_n}{V}_{n(j),0(i)}^{d_R}{d}_{R,n}^{(j)},$$

$$(m_{d})_{kl}\simeq V^{q_{L}}_{0(i),0(k)}\big{(}Y_{0}^{d}\big{)}_{ij}V_{0(j),0(l)}^{d_{R}}\frac{v}{\sqrt{2}},\qquad(m_{d})_{kl}\simeq V^{q_{L}}_{0(i),0(k)}\big{(}Y_{0}^{u}\big{)}_{ij}V_{0(j),0(l)}^{u_{R}}\frac{v}{\sqrt{2}}.$$

$$f_{q_{L}^{(i)}}\sim\Big{(}\frac{M^{q}}{Y^{q}\langle\phi\rangle}\Big{)}^{N_{q(i)}},\qquad f_{u_{R}^{(i)}}\sim\Big{(}\frac{M^{u}}{Y^{u}\langle\phi\rangle}\Big{)}^{N_{u(i)}},\qquad f_{d_{R}^{(i)}}\sim\Big{(}\frac{M^{d}}{Y^{d}\langle\phi\rangle}\Big{)}^{N_{d(i)}},$$

$$\begin{array}{rl}{\hat{N}}_q{|}_n=& \{3,2,2,1\},\\ {\hat{N}}_u{|}_n=& \{3,2,2,1\},\\ {\hat{N}}_d{|}_n=& \{3,3,3,3\}.\end{array}$$

$$\begin{array}{rl}{N}_{q(1)}& =3,N_{q(2)}=2,N_{q(3)}=0,\\ N_{u(1)}& =3,N_{u(2)}=2,N_{u(3)}=0,\\ N_{d(1)}& =3,N_{d(2)}=3,N_{d(3)}=3,\end{array}$$

$${\mathcal{L}}_{\dim 5}\supset c_{ij}(L_{i}H)(L_{j}H)/{\Lambda }_{\mathrm{LN}},$$

$$m_{ij}^{\nu}\simeq c_{ij}\frac{v^{2}}{\Lambda_{\rm LN}}\Big{(}\frac{M^{L}}{\langle\phi\rangle}\Big{)}^{N_{L(i)}+N_{L(j)}},$$

$$N_{L(1)}=N_{L(2)}=N_{L(3)}.$$

$$m_{ij}^{e}\sim v\Big{(}\frac{M}{\langle\phi\rangle}\Big{)}^{N_{L(i)}+N_{e(j)}}.$$

$$N_{e(1)}=N_{e(2)}+3=N_{e(3)}+4.$$

$$N_{L(1)}=N_{L(2)}=N_{L(3)}=3,N_{e(1)}=4,N_{e(2)}=1,N_{e(3)}=0.$$

$$N_{L(i)}\equiv N_L=2,N_{e(i)}\equiv N_e=3,$$

$$m_{ij}^{e}\simeq\frac{v}{\sqrt{2}}M_{N_{L}}^{L}\big{(}Y_{N_{L}}^{L}\langle\phi\rangle\big{)}^{-1}\cdots\big{(}Y_{1}^{L}\langle\phi\rangle\big{)}^{-1}Y_{0}^{e}\big{(}Y_{1}^{e}\langle\phi\rangle\big{)}^{-1}M_{1}^{e}\cdots\big{(}Y_{N_{e}}^{e}\langle\phi\rangle\big{)}^{-1}M_{N_{e}}^{e}.$$

Full text

1 Introduction

The masses of the Standard Model (SM) quarks and leptons exhibit large hierarchies Tanabashi:2018oca . 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 Froggatt:1978nt (see also Leurer:1992wg ; Leurer:1993gy ). In FN models the SM mass hierarchy is due to a horizontal $U(1)_{\rm FN}$, spontaneously broken by the vacuum expectation value (vev) of the flavon field, $\langle\phi\rangle$. The SM fermions carry horizontal $U(1)_{\rm FN}$ charges. The up and down quark Yukawa matrices are then given by111To shorten the discussion we focus manly on quarks, but also show in Section 4 how the discussion naturally generalizes to leptons.

[TABLE]

with $M$ the typical mass of the new FN vector-like fermions, while $[q_{i}],[u_{i}],[d_{i}]$ are the horizontal charges of the SM left-handed quark doublets, the up quarks, and the down quarks, respectively, with $i,j=1,2,3$ the generational indices. For $\langle\phi\rangle/M\sim 0.2$ one can reproduce the parametric size of the CKM matrix and the quark masses with the appropriate choice of horizontal charges, $[q_{i}],[u_{i}],[d_{i}]$ (for details, see below). For generic horizontal charges, and without additional field content, such a horizontal $U(1)_{\rm 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. Froggatt:1998he builds neutrino mass models based on an anomaly free $U(1)_{\rm FN}^{2}$ FN symmetry, while the full list of anomaly free charge assignments for $U(1)_{\rm FN}$ can be found in Allanach:2018vjg ; Costa:2019zzy (for related work, see also Ellis:2017nrp ; DelleRose:2017xil ; Correia:2019woz ; Correia:2019pnn ; Berger:2000sc ; Chen:2006hn ; Rathsman:2019wyk ). The gauged $U(1)_{\rm FN}$ can also be left naively anomalous, if the cancellation of $U(1)_{\rm FN}$ gauge anomalies occurs through the Green-Schwarz mechanism Ibanez:1994ig ; Binetruy:1994ru ; Jain:1994hd ; Dudas:1995yu ; Dudas:1996fe ; Binetruy:1996xk ; Irges:1998ax ; King:1999mb ; King:1999cm ; Shafi:2000su ; Dreiner:2003yr ; Dreiner:2007vp .

In this paper we explore a twist to the original FN proposal. In our modified set of FN models the horizontal $U(1)_{\rm 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

[TABLE]

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 $\langle\phi\rangle/M\sim 5$. Such a realization of the FN flavor model was discussed in Ref. Alonso:2018bcg , where it was given as a possible UV realization of the clockwork mechanism Giudice:2016yja . We extend the results of Ref. Alonso:2018bcg in several ways. First of all, two choices for the horizontal symmetry group, $G_{\rm FN}=U(1)_{\rm FN}^{3}$ and $G_{\rm FN}=U(1)_{\rm FN}$, are explored. The choice $G_{\rm FN}=U(1)_{\rm FN}^{3}$ most closely resembles the analysis in Ref. Alonso:2018bcg , except that we allow ${\mathcal{O}}(1)$ variations for all the Yukawa couplings. We also assume that $G_{\rm FN}$ is gauged, and explore the phenomenology of the associated $Z^{\prime}$ gauge bosons.

The typical spectra of inverted FN models for the two choices of $G_{\rm 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 ${\mathcal{O}}(10^{7}~{}{\rm GeV})$. The mass of $Z^{\prime}$, however, can span many orders of magnitude. For gauge couplings $g^{\prime}\sim{\mathcal{O}}(1)$ the $Z^{\prime}$ is heavy, while for small values of $g^{\prime}$ 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^{\prime}$ searches.

The paper is organized as follows. In Section 2 we present the inverted FN model based on $U(1)_{\rm FN}^{3}$ 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 (2) can be understood using a spurion analysis for an approximate $U(1)_{\rm app}$ that emerges at low energies, and is broken by $M$ rather than $\langle\phi\rangle$. In subsection 2.2 we then perform a numerical scan over the parameters of the $U(1)_{\rm FN}^{3}$ model and compare the results with the measured values of quark masses and mixings. In Section 3 we present the $U(1)_{\rm 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^{\prime}_{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-\bar{K}$ bounds for flavorful light $Z^{\prime}$. The representative benchmarks for the two models are given in Appendix B.

2 Decoupled FN chains - the $U(1)_{\rm FN}^{3}$ model

2.1 The set-up of the model

We start by assuming that there is a separate horizontal $U(1)_{\rm FN}$ for each generation of the SM fermions, i.e., that the horizontal gauge group is $G_{\rm FN}=U(1)_{1}\times U(1)_{2}\times 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_{\rm FN}=U(1)_{\rm FN}$ for now. The Lagrangian is

[TABLE]

where

[TABLE]

The $q_{L/R,j}$ and $d_{L/R,j}$ have the SM quantum numbers of $q_{L}$ and $d_{R}$ SM quark fields, respectively. The flavon fields, $\phi$, carry the horizontal $U(1)_{\rm FN}$ charge of $[\phi]=+1$, while $d_{R/L,n}$ and $q_{R/L,n}$ have charges $[d_{R/L,n}]=+n$ and $[q_{R/L,n}]=-n$, respectively. The above Lagrangian can be represented as a chain of nodes that contain fields with the same $U(1)_{\rm FN}$ charges, connected by flavon fields, see Fig. 2. For simplicity we set all the vector-like masses to be equal, $M$. This can be trivially relaxed, without affecting the main properties of the model as long as they remain ${\mathcal{O}}(M)$ (without loss of generality they can all be made real positive through phase redefinitions). More important is the assumption that the flavon couples on each node only to one field out of each pair of vector-like fermions, cf. Fig. 2. In general, there could also be terms of the form $\phi^{*}\bar{d}_{L,n-1}d_{R,n}$ and $\phi^{*}\bar{q}_{R,n-1}q_{L,n}$, which were dropped in (4) and (5). 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, $\langle\phi\rangle$, which spontaneously breaks the horizontal $U(1)_{\rm 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 ${\mathcal{O}}(M,\langle\phi\rangle)$. For $\langle\phi\rangle\gg M$, with $Y_{n}^{d},Y_{n}^{q}\sim{\mathcal{O}}(1)$ the two zero modes are mostly localized towards the ends of the respective chains. More precisely, the $q_{L,0}^{\prime}$ and $d_{R,0}^{\prime}$ massless mass eigenstates are given by

[TABLE]

The notation is borrowed from Alonso:2018bcg , with $V^{q_{L}}$ and $V^{q_{R}}$ the $N_{q}\times N_{q}$ and $(N_{q}+1)\times(N_{q}+1)$ unitary matrices that diagonalize the $N_{q}\times(N_{q}+1)$ mass matrix for the $q$ fields through a bi-unitary transformation (and similarly for $d_{L,R}$ fields with obvious change of notation). For $Y_{n}^{d}=Y_{n}^{q}=1$, $M_{n}^{d}=M_{n}^{q}=M,$ the model realizes the clockworking mechanism, for which the zero mode profiles are known analytically Alonso:2018bcg ,

[TABLE]

with the normalization constants {\cal N}_{0}^{\psi}=\sqrt{\big{(}q^{2}-1\big{)}/(q^{2}-1/q^{2N_{\psi}})}\simeq 1, where $q=\langle\phi\rangle/M\gg 1$.

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,

[TABLE]

where the overlaps with the zero-node are given by

[TABLE]

and are thus exponentially suppressed for large $N_{\psi}$ and $\langle\phi\rangle\gg 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)_{\rm 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)_{\rm FN}$ charges that differ by one unit, but have the same $U(1)_{\rm app}$ charges. The $U(1)_{\rm app}$ charge reduces by one unit when hopping to the right by one node in the chain. That is, the $q_{L,N_{q}}$ has the $U(1)_{\rm app}$ charge $N_{q}$, the $q_{L,N_{q}-1}$ and $\bar{q}_{L,N_{q}}$ the $U(1)_{\rm 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)_{\rm app}$). The mass terms proportional to $M$ explicitly break $U(1)_{\rm app}$. Treating $M$ as a spurion with $U(1)_{\rm app}$ charge of $+1$ the fermion mass matrix is formally invariant under $U(1)_{\rm app}$. The zero modes, i.e., the SM fermions are mostly $q_{L,N_{q}}$ and $d_{R,N_{d}}$, so that the zero mode mass term is given by

[TABLE]

which is $U(1)_{\rm app}$ invariant. This is exactly the suppression found in (8). Furthermore, note that the effective FN charge of the zero mode is approximately the same as its $U(1)_{\rm app}$ charge. The zero modes are localized toward the ends of the vector-like chains, see Fig. 5. To leading order the $q_{L,N_{q}}$ 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_{\rm FN}=U(1)_{1}\times U(1)_{2}\times U(1)_{3}$, with a separate $U(1)_{i}$ for each of the three generations of fermions. The fermion Lagrangian is then

[TABLE]

with $i,j=1,2,3,$ the generation indices, and $\tilde{H}=i\sigma_{2}H^{*}$. The kinetic Lagrangians are given in (4) and (5), but with extra generational superscript on fermions, $q_{a}^{(i)},u_{a}^{(i)},d_{a}^{(i)}$ (replacing $d\to u$ in ${\cal L}_{d{(i)}}$ gives ${\cal L}_{u{(i)}}$), and with $M_{n}^{q}\to M_{n}^{q(i)}$, $Y_{n}^{q}\to Y_{n}^{q(i)}$, $N_{q}\to N_{q(i)}$, etc., as well as $\phi\to\phi_{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_{0}^{d},Y_{0}^{u}$ on the zero node are $3\times 3$ complex matrices, while the Yukawa couplings on each of the chains, (4) and (5), are just complex numbers.222The flavon Yukawa couplings $Y_{1}^{q(i)}$, connecting the zero node with the first node of the $q-$chain, are in general arbitrary complex three-vectors. To simplify the notation in (11) we take them to be orthogonal, and then use the $3\times 3$ unitary field redefinitions of $\bar{q}_{L,0}^{(i)}$ to ensure the $Y_{1}^{q(i)}$ have nonzero components only in the $i$-th direction. We assume the same for $Y_{1}^{d(i)}$ and $Y_{1}^{u(i)}$. The Yukawa interactions with the Higgs are thus the only terms mixing different generations, see the last term in (11).

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. (7) and Fig. 5. This means that the SM fermion, $q_{L}^{(i)}$, has the largest component that is due to $q_{L,N_{q}}^{(i)}$, the left-handed part of the vector-like fermion that is the furthest away from the Higgs. For $u_{R}^{(i)}$ the largest component is $u_{R,N_{u}}^{(i)}$, and similarly for $d_{R}^{(i)}$ it is $d_{R,N_{d}}^{(i)}$. The SM fermions thus carry large $U(1)_{i}$ charges, and couplings to the $U(1)_{i}$ gauge bosons, $Z_{i}^{\prime}$s, roughly given by the length of the corresponding vector-like chain, $[\psi]_{i}\approx N_{\psi{(i)}}$. However, for the phenomenology of $Z_{i}^{\prime}$ 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

[TABLE]

The structure of the zero node overlaps, $f_{a}$, explains the hierarchy of SM fermion masses. For $M\ll\langle\phi\rangle$ 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. (9), 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 (12) resemble the observed values. To simplify the discussion, we set

[TABLE]

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)_{\rm FN}$ case, to be covered in Section 3. The scalings for the mass matrices in (12) are then, for $Y_{0}^{d},Y_{0}^{u}\sim{\mathcal{O}}(1)$,

[TABLE]

2.2 Numerical scan

For the numerical scan we set $\langle\phi\rangle=10^{7}$ GeV and $q=5$, cf. Eq. (13). The vector-like masses, $M_{n}^{f}$, $f=u,d,q$, Eqs. (4), (5), in the units of $10^{7}$ GeV, and the zero node Yukawa matrix elements, $(Y_{0}^{d})_{ij},(Y_{0}^{u})_{ij}$, Eq. (11), are taken to be random complex numbers of the form $r_{a}e^{i\varphi_{a}}$. The Yukawa couplings between the FN fermions and the flavon are taken to be of the form $Y_{n}^{f}\langle\phi\rangle=r_{a}e^{i\varphi_{a}}q$, so that the ratios $|Y^{f}\langle\phi\rangle/M^{f}|$ are on average equal to $q$. The magnitudes and phases are taken to be uniformly distributed over $r_{a}\in[0.3,0.9]$ and $\varphi_{a}\in[0,2\pi)$, with $r_{a},\varphi_{a}$ uncorrelated between different couplings. The range of $r_{a}$ was chosen such that the top Yukawa, $\bar{m}_{t}/(v/\sqrt{2})\simeq 0.55$ is close to the median, while the ratio of the boundaries, $r_{\rm max}/r_{\rm 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_{n}^{f}/(Y_{n}^{f}\langle\phi\rangle)$, 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_{0}^{d,u}$ matrix elements. For instance, using $r_{a}\in[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. Alonso:2018bcg ,333The difference arises partially from taking significantly higher $\mu$ and different median values of $Y_{0}^{d,u}$ matrix elements. We take $N_{d(3)}=3$, $N_{u(2)}=2$ while in Alonso:2018bcg 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.

[TABLE]

giving the configuration of the vector-like FN chains shown in Fig. 4.

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 $\overline{\text{MS}}$ masses at $\mu=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 $\mu=10^{4}$ TeV, obtained through NNLO RG evolution, are listed in Table 1 (at $\mu=2$ TeV the values of the quark masses coincide with Ref. Alonso:2018bcg ). For the CKM matrix elements we checked using results in Babu:1987im that the effects of RG evolution are negligible, so that we use the unevolved values from the PDG Tanabashi:2018oca .

3 Coupled FN chains - the $U(1)_{\rm FN}$ model

3.1 The set-up of the model

We discuss next the coupled FN chains, which are obtained in the case of a single horizontal $U(1)_{\rm FN}$, see Fig. 4. The Yukawa couplings between vector-like fermions and the flavon field are now complex matrices. The Lagrangian is thus

[TABLE]

where

[TABLE]

and similarly for ${\cal L}_{u}$ with the $d\to u$ replacement. For ${\cal L}_{q}$, on the other hand, there is one more left-handed field for each generation than the right-handed fields, so that

[TABLE]

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,\ldots,N_{q(1)}$ sums over all of the nodes. On node $n$ there are $\hat{N}_{q}|_{n}$ generations, which are coupled to the $\hat{N}_{q}|_{n-1}$ generations of fermions on the $(n-1)$-th node through an $\hat{N}_{q}|_{n-1}\times\hat{N}_{q}|_{n}$ complex Yukawa matrix $Y_{n}^{q}$.

This produces a set of coupled FN chains, where the mixing between generations is not only through $3\times 3$ Yukawa matrices on the zero node, coupling to the Higgs, but also through the flavon Yukawa couplings, \big{(}Y^{q}_{n}\big{)}_{ji}, and vector-like masses, \big{(}M_{n}^{q}\big{)}_{ji}, cf. Fig. 7. The zero modes are obtained through unitary transformations

[TABLE]

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\times 3$ matrices, $V_{0(j),0(i)}^{d_{R},u_{R}}$, due to the inter-generational mixing in the chains. Therefore, the quark mass matrices become

[TABLE]

For $\langle\phi\rangle\gg 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,

[TABLE]

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. (12). 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. (21), part of the hierarchy comes from the properties of random matrix multiplications vonGersdorff:2017iym (see also Section 4).

3.2 Numerical scan

We choose the chain configurations with the following number of fermion generations on each node

[TABLE]

Here $n=0,\ldots,3$, respectively, for each of the three cases. This implies that the lengths of the vector-like chains are given by,

[TABLE]

where $N_{a}+1$ is the length of the corresponding chain, cf. Eq. (17). The above chain configuration differs from the one for the decoupled FN chains in (15), 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$, $\langle\phi\rangle=10^{7}$ GeV, and vary the matrices $M_{ij}^{f},(Y_{n}^{f})_{ij}\langle\phi\rangle$, $f=u,d,q$, cf. Eqs. (17) (18), and $(Y_{0}^{d})_{ij},(Y_{0}^{u})_{ij}$, cf. Eq. (16). The vector-like masses are taken to be random complex numbers of the form $r_{a}e^{i\varphi_{a}}$ in the units of $10^{7}$ GeV, while the elements of Yukawa matrices are taken to be of the form $(Y_{n}^{f})_{ij}\langle\phi\rangle=r_{a}e^{i\varphi_{a}}q$. The $r_{a}$ are varied in the range $[0.3,0.9]$ and phases in $[0,2\pi)$. 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 vonGersdorff:2017iym . The resulting mass distributions are much broader than in the $U(1)_{\rm FN}^{3}$ 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 completely overlap. This is a result of relatively small expansion parameter $1/q\simeq 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.

4 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_{\rm FN}=U(1)_{\rm FN}^{3}$ case the Lagrangian for the leptons is then given by replacing $Q\to L$, $d\to e$ in Eq. (11), and ignore terms that involve $u$. Performing the same replacements in Eq. (16) gives the Lagrangian for the $G_{\rm FN}=U(1)_{\rm 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

[TABLE]

with $c_{ij}\sim{\mathcal{O}}(1)$. In the $\langle\phi\rangle\gg M_{n}^{L}$ limit the neutrino mass matrix takes the form

[TABLE]

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., CRZZ:2019 . We focus on the case where the neutrino masses are completely anarchic, with all the PMNS mixing angles taken to be ${\mathcal{O}}(1)$. This happens if all the left-handed leptons have the vector-like fermion chains of the same length,

[TABLE]

For $G_{\rm FN}=U(1)_{\rm FN}^{3}$ the hierarchy among charged lepton masses is then due to different lengths of FN chains for the right-handed leptons, giving

[TABLE]

The observed hierarchy between $m_{e}:m_{\mu}:m_{\tau}$ is obtain for

[TABLE]

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

[TABLE]

Fig. 11 shows the resulting values of PMNS matrix elements and of the neutrino masses, setting $\Lambda_{\rm LN}=10^{11}$ GeV in order to approximately saturate the bound on the sum of neutrino masses from cosmology $\sum_{i}m_{\nu_{i}}\lesssim 0.15$ Vagnozzi:2017ovm . These are compared with the measured values with solid (dashed) lines denoting the central values (1$\sigma$ bands) Tanabashi:2018oca . 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}\in[0.3,0.9]$ in units of $10^{7}$ GeV, with arbitrary phases, while the Yukawa couplings $Y_{n}^{f}\langle\phi\rangle$ are equal to $q$ up to a randomized complex prefactor with magnitude in the range $[0.3,0.9]$ and a random phase, where $\langle\phi\rangle=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_{\rm FN}=U(1)_{\rm FN}$, we exploit the fact that products of random matrices have hierarchical eigenvalues. We find that a completely anarchic charge assignment

[TABLE]

describes well the hierarchy among the charged leptons, see Fig. 10 (right). In the $\langle\phi\rangle\gg M$ limit the charged lepton mass matrix is then

[TABLE]

For $(M^{L,e}_{n})_{ij}=(c_{n}^{L,e})_{ij}q\langle\phi\rangle$, with the $Y^{L,e}_{n}$ and $c^{L,e}_{n}$ matrix elements randomly distributed, such that their average values vanish, one has vonGersdorff:2017iym

[TABLE]

where for simplicity we have assumed that $(c_{n}^{L,e})^{-1}$ and $Y_{n}^{L,e}$ all follow the same distribution with variance $\sigma$. 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. (31). From (32) and (33) we get

[TABLE]

If the lengths of the FN chains grow the eigenvalues become more hierarchical following (34), 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. (30) 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 $\Lambda_{\text{LN}}=10^{13}$ GeV are shown in Fig. 12. We observe that there is a tendency for $U_{e,2},U_{\mu,3},U_{\tau,2},U_{\tau,3}$ to be below the observed values, and for $U_{e,1},U_{\mu,2},U_{\tau,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.

5 The phenomenology of flavorful $Z^{\prime}$ 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-\bar{B}$, $K-\bar{K}$, mixing, $\mu\to e\gamma$, etc. The new contributions are due to the exchanges of flavons, heavy vector-like fermions, and, if $G_{\rm FN}$ is gauged, also flavorful $Z$’s. Below we will see that the FCNC constraints on possible tree level $Z^{\prime}$ exchanges bound $\langle\phi\rangle\gtrsim{\mathcal{O}}(10^{7}{\rm~{}GeV})$. This in turn means that the vector-like fermions are heavier than about $10^{7}$ GeV, since for ${\mathcal{O}}(1)$ Yukawa couplings their masses are comparable to the flavon vev, $\langle\phi\rangle$.

We will derive the FCNC constraints on inverted FN models assuming the $Z^{\prime}$ contributions dominate. That is, we will assume that both the vector-like fermions and the flavon have masses $m_{i}\sim\langle\phi\rangle$. As a result, their contributions to FCNC transitions are subleading compared to the $Z^{\prime}$ contributions. The vector-like fermions contribute to FCNCs only at one-loop, while the FCNC flavon couplings to the SM fermions $f_{i},f_{j}$ are suppressed by $\sqrt{m_{f_{i}}m_{f_{j}}}/\langle\phi\rangle$.

In the construction of inverted FN models we made two choices for the anomaly free horizontal group $G_{\rm FN}$. If $G_{\rm FN}$ is gauged, there is one extra gauge boson, $Z^{\prime}$, for the $G_{\rm FN}=U(1)_{\rm FN}$ case, and three new gauge bosons, $Z_{i}^{\prime}$, $i=1,2,3$ for the $G_{\rm FN}=U(1)_{\rm FN}^{3}$ case. The rest of this section is devoted to the phenomenology of these flavorful $Z^{\prime}$s. If the $G_{\rm FN}$ gauge couplings are small the $Z^{\prime}$ can be light. This means that beside the indirect searches using FCNCs, the $Z^{\prime}$ can also be searched for in on-shell production, e.g., in beam dumps or in astrophysical environments.

We first focus on the $G_{\rm FN}=U(1)_{\rm FN}$ case. After the $U(1)_{\rm FN}$ is broken the relevant terms in the Lagrangian are

[TABLE]

where $B_{\mu\nu}$ and $Z_{\mu\nu}^{\prime}$ are the field strength tensors of $U_{Y}(1)$ and $U(1)_{\rm FN}$, respectively, with $B_{\mu}$, $Z_{\mu}^{\prime}$ the corresponding gauge bosons, while $W_{\mu}^{a}$ are the $SU(2)_{L}$ gauge bosons. Since the FN fermions are charged both under the SM gauge group and the $U(1)_{\rm FN}$, the fermionic kinetic terms give the couplings of the $B_{\mu},W_{\mu}^{a}$ and $Z_{\mu}^{\prime}$ to the chains of FN fermions. These result in the last three terms in (35), where

[TABLE]

while $J_{W^{a}}^{\mu}$ is obtained from $J_{Y}^{\mu}$ by replacing $g_{Y}Y_{q}$ with $gT^{a}$, 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 (36), (37) by making replacements $q_{L}\to L_{L},e_{R},u_{R},d_{R}$ and $q_{R}\to L_{R},e_{L},u_{L},d_{L}$, respectively. Here $g_{Y}$ and $g^{\prime}$ are the $U_{Y}(1)$ and $U(1)_{\rm FN}^{\prime}$ gauge couplings, while $Y_{f}$ and $n\delta_{f}$ are, respectively, the hypercharge and the horizontal $U(1)_{\rm FN}$ quantum numbers of the fermion $f_{L/R,n}^{(i)}$, where $\delta_{q}=\delta_{L}=-1$ and $\delta_{u}=\delta_{d}=\delta_{e}=1$.

The couplings to the SM fermions are obtained by performing the unitary transformations in (19), and keeping only the zero modes. This gives

[TABLE]

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}^{\mu}$ 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_{\rm FN}^{\mu}$, has flavor violating couplings,

[TABLE]

and similarly for the other $c^{\prime}{}_{\negmedspace ij}^{f}$ 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^{\prime}{}_{\negmedspace ij}^{q_{L}}$ to have eigenvalues that are ${\mathcal{O}}(N_{q(i)})$, see also the discussion below (10).

We can get rid of the kinetic mixing between $B_{\mu}$ and $Z_{\mu}^{\prime}$, Eq. (35), by performing a field redifinition $B_{\mu}\to B_{\mu}-\epsilon Z^{\prime}_{\mu}$.444To have canonically normalized fields one also needs to rescale $Z^{\prime}$ at ${\mathcal{O}}(\epsilon^{2})$, which we can safely ignore since we work to ${\mathcal{O}}(\epsilon)$. This induces a coupling of $Z^{\prime}$ with $J_{Y}^{\mu}$ proportional to $\epsilon$,

[TABLE]

After the field redefinition the covariant derivative acting on the Higgs also contains $Z^{\prime}$, D_{\mu}H=\big{(}\partial_{\mu}+ig_{Y}(B_{\mu}-\epsilon Z_{\mu}^{\prime})/2+igT^{a}W_{\mu}^{a}\big{)}H. After FN symmetry and electroweak symmetry are broken, $\phi\to\langle\phi\rangle$, $H\to(0,v/\sqrt{2})$, the scalar kinetic terms,

[TABLE]

mix the $B_{\mu},W_{\mu}^{3}$ and $Z_{\mu}^{\prime}$ (here D_{\mu}\phi=\big{(}\partial_{\mu}+ig^{\prime}Z_{\mu}^{\prime}\big{)}\phi). This mixing can be thought of as occurring in two steps. In the $\epsilon\to 0$ limit the electroweak breaking mixes $B_{\mu},W_{\mu}^{3}$ into a massless photon, $A_{\mu}$, and the massive $Z_{\mu}$. For nonzero $\epsilon$ the $Z_{\mu}$ and $Z_{\mu}^{\prime}$ further mix into two mass eigenstates, $\hat{Z}_{\mu}^{\prime}=c_{\theta}Z_{\mu}^{\prime}+s_{\theta}Z_{\mu}$, $\hat{Z}_{\mu}=-s_{\theta}Z_{\mu}^{\prime}+c_{\theta}Z_{\mu}$, where $s_{\theta}=\sin\theta$, $c_{\theta}=\cos\theta$. The mixing angle is, up to ${\mathcal{O}}(\epsilon^{2})$ corrections,

[TABLE]

with $s_{W}=\sin\theta_{W}$ the sine of the weak mixing angle. The mass of $\hat{Z}$ is the same as for the SM $Z$, $m_{Z}$, up to ${\mathcal{O}}(\epsilon^{2})$ corrections, while the $\hat{Z}^{\prime}$ has the mass

[TABLE]

Note that for $m_{Z}\gg g^{\prime}\langle\phi\rangle$ we have $\theta\to-s_{W}\epsilon$ and $m_{Z^{\prime}}\to\sqrt{2}g^{\prime}\langle\phi\rangle$.

We can finally write down the couplings of $Z^{\prime}$ to the SM fermions,

[TABLE]

where the sum runs over all the SM fermions, $f=u,d,\ell,\nu$, with $c_{\nu_{R}}^{ij}=0$, and $i,j=1,\ldots,3,$ the generation indices. The couplings receive two contributions, the flavor diagonal one from $J_{Y}^{\mu}$ and $J_{W^{3}}^{\mu}$, while the contribution from the horizontal current, $J_{\rm FN}^{\mu}$, also contains the flavor violating couplings,

[TABLE]

and similarly for $d_{L,R}$ with $u\to d$ replacements, for $\ell_{L,R}$ with $q\to L$, $u\to\ell$ replacements, and for $\nu_{L}$ with $q\to L$, $u_{L}\to\nu_{L}$ replacements in the above expressions. For later convenience we also introduce the vector and axial couplings

[TABLE]

The flavor diagonal $\hat{\epsilon}_{f}$ term is of ${\mathcal{O}}(\epsilon)$,

[TABLE]

where $T_{3}^{f}$ is the weak isospin for fermion $f$. Note that the second term vanishes in the $m_{Z}\to\infty$ limit, cf. Eq. (42), while the first term is the contribution from the kinetic mixing between $Z^{\prime}$ and the photon. The unitary matrices $V_{f}$ in (45) 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^{\prime}$ to the SM fermions. If the $c_{f}$ are dominated by $\hat{\epsilon}_{f}$, then the phenomenology of $Z^{\prime}$ is the same as for the dark photon. In the opposite limit, when $\hat{\epsilon}_{f}$ is negligible, the couplings of the $Z^{\prime}$ are governed by the $U(1)_{\rm FN}$ charges, giving both flavor diagonal and off-diagonal couplings of comparable strength. In the numerical examples below we set $\epsilon\to 0$, see Table 2 and Appendix B. In this limit the $Z^{\prime}$ mostly couples through axial vector couplings, since left-handed and right-handed zero modes carry opposite effective $U(1)_{\rm FN}$ charges, cf. Appendix B.

The above derivations generalize straightforwardly to the case of $G_{\rm FN}=U(1)_{\rm FN}^{3}$, i.e., the decoupled FN chains. The $Z^{\prime}$ terms in the $G_{\rm FN}=U(1)_{\rm FN}$ Lagrangian, Eq. (35), are replaced by

[TABLE]

where $J^{\mu}_{i,FN}$ are the currents corresponding to the horizontal $U(1)_{i}$ symmetries for the $i$-th generation,

[TABLE]

with $g_{i}^{\prime}$ the gauge coupling of $U(1)_{i}$. The last term in Eq. (48) contains kinetic mixings between different $Z_{i}^{\prime}$. These enter the interactions with the SM fermions only at order ${\mathcal{O}}(\epsilon g^{\prime})$ or ${\mathcal{O}}(\epsilon^{2})$ and can be safely ignored. The field redefinition $B_{\mu}\to B_{\mu}-\sum_{i=1}^{3}\epsilon_{i}Z^{\prime}_{i,\mu}$ gets rid of the kinetic mixing between $B_{\mu}$ and $Z_{i}^{\prime}$ and trades it for mixing between $Z$ and $Z_{i}^{\prime}$ through the gauge boson mass terms, mirroring the discussion for the $G_{\rm FN}=U(1)_{\rm FN}$ case above. This results in ${\mathcal{O}}(\epsilon)$ flavor diagonal couplings of $Z_{i}^{\prime}$ to the SM fermions, while the $J^{\mu}_{i,{\rm FN}}$ result in both flavor diagonal and off-diagonal couplings. In general we can write

[TABLE]

where the coefficients have the same general form as in Eq. (45), but now for each $Z_{i}^{\prime}$ separately, i.e.,

[TABLE]

The $V_{u_{L}}$ and $V_{u_{R}}$ diagonalize the up quark mass matrix, while the $c^{\prime}{}_{\negmedspace k}^{q_{L}}$ matrix has only one nonzero entry,

[TABLE]

and similarly for the other $c^{\prime}{}_{\negmedspace k}^{f}$ matrices. Since the zero modes are localized toward the ends of vector-like chains we expect \big{(}c^{\prime}{}_{\negmedspace k}^{q_{L}}\big{)}_{kk}\sim N_{q(k)}. The $\hat{\epsilon}_{f,k}$ coefficients are of ${\mathcal{O}}(\epsilon_{i})$ and vanish in the limit $\epsilon_{i}\to 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

[TABLE]

The bounds on flavorful $Z^{\prime}$s come from a variety of experimental observables. They can be grouped into four broad categories: the bounds that come from $Z^{\prime}$ couplings to quarks or from $Z^{\prime}$ 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^{\prime}$ couplings for the two benchmarks are given in Appendix B. In both benchmarks we take the $\epsilon\to 0$ limit, therefore the $Z^{\prime}$ couplings are completely dictated by $J_{\rm FN}^{\mu}$, Eq. (38). The various experimental bounds in the $(m_{Z^{\prime}},g^{\prime})$ plane are compiled in Figs. 13-18. Note that the couplings of $Z^{\prime}$ are mostly axial, cf. Table 2. For other phenomenological analysis of light axial vectors, see Kahn:2016vjr ; Kozaczuk:2016nma .

5.1 The bounds on flavorful $Z^{\prime}$ for the $G_{\rm FN}=U(1)_{\rm FN}$ benchmark

We first derive the experimental bounds on the $G_{\rm FN}=U(1)_{\rm FN}$ benchmark, see Appendix B.2. As we will show below, the $K-\bar{K}$ mixing bounds the flavon vev to be very large, $\langle\phi\rangle\gtrsim 10^{7}$ GeV. This means that the FN fermions are very heavy, with masses $m_{F}\sim{\mathcal{O}}(\langle\phi\rangle)$. 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_{\phi}\sim{\mathcal{O}}(\langle\phi\rangle)$, giving $\sim m_{d_{i}}m_{d_{j}}/m_{\phi}^{2}$ suppressed contributions to the flavor observables compared to the $Z^{\prime}$ and can thus be ignored. This assumption can be relaxed in the future, since $m_{\phi}$ is a free parameter. To make the analysis tractable, we also limit the mass of the $Z^{\prime}$ to be above, $m_{Z^{\prime}}>10$ MeV, an assumption that could also be relaxed in future studies. The results obtained for $G_{\rm FN}=U(1)_{\rm FN}$ are straightforward to extend to the $U(1)_{\rm FN}^{3}$ model, which we do in Section 5.2.

5.1.1 Flavor diagonal couplings to quarks and/or leptons

We start by deriving constraints on $Z^{\prime}$ couplings from flavor conserving processes. Representative Feynman diagrams are shown in Fig. 15, while the numerical values of $Z^{\prime}$ couplings with fermions in our benchmark are listed in Appendix B.2. The $Z^{\prime}$ couplings are almost flavor diagonal, see Eqs. (137)-(140). Furthermore, the couplings to top and bottom are highly suppressed, except the coupling to $b_{R}$. This is easy to understand from the $U(1)_{\rm FN}$ charge assignments, since the $q_{L}^{(3)}$ and $u_{R}^{(3)}$ are not charged under $U(1)_{\rm FN}$, see Eq. (23) and Fig. 7.

Direct $Z^{\prime}$ production:

A number of experiments have searched for a dark photon, a heavy gauge boson that kinematically mixes with the hypercharge. These searches were recast in Ref. Ilten:2018crw for a generic $Z^{\prime}$ 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^{\prime}$ from direct production, adapting it to the case in hand. The $U(1)_{\rm FN}$ $Z^{\prime}$ has predominantly axial vector couplings, cf. Table 2, unlike dark photon, which only has flavor diagonal vector couplings. For $e^{+}e^{-}\to Z^{\prime}\gamma$ searches at BaBar Lees:2014xha ; Lees:2017lec , KLOE Anastasi:2015qla and LEP Abdallah:2003np ; Abdallah:2008aa in $Z^{\prime}\to e^{+}e^{-},\mu^{+}\mu^{-}$,inv channels, as well as for the $Z^{\prime}$ bremsstrahlung Bjorken:2009mm ; Andreas:2012mt ; Blumlein:2013cua in electron beam dump searches at A1 Merkel:2014avp , APEX Abrahamyan:2011gv , E137 Bjorken:1988as , E141 Riordan:1987aw , E774 Bross:1989mp , Orsay Davier:1989wz , KEK Konaka:1986cb , NA64 Banerjee:2018vgk , and in proton beam dump search at $\nu$-CAL I Blumlein:1991xh , with $Z^{\prime}\to e^{+}e^{-}$, the Darkcast recast of bounds applies to our $Z^{\prime}$ benchmark with the replacement $\epsilon eQ_{f}\to((c_{fV}^{11})^{2}+(c_{fA}^{11})^{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. Ilten:2018crw . In the same way, the LHCb searches for dark photon Aaij:2017rft can be recast for $m_{Z^{\prime}}>1$ GeV region, where the production is dominated by Drell-Yan production, $\bar{q}_{i}q_{i}\to Z^{\prime}$. We can also safely neglect off-diagonal couplings in direct $Z^{\prime}$ 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^{\prime}$; by LHCb Aaij:2017rft for $m_{Z^{\prime}}<1$ GeV, by NA60 Arnaldi:2016pzu , CHARM Bergsma:1985qz , $\nu$-CAL I Blumlein:1990ay and KLOE Archilli:2011zc . For instance, for dark photon with mass below $1$ GeV the production in $pp$ collision is dominated by production from $\pi^{0}\to\gamma_{d}\gamma,\eta\to\gamma_{d}\gamma,\omega\to\gamma_{d}\pi^{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 $\pi^{0}\to Z^{\prime}\gamma$ see Kahn:2016vjr ). 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.

Atomic Parity Violation (APV):

The APV measurements bound the parity violating combination of couplings $g^{\prime}{}^{2}|c_{\ell A}^{11}c_{N}^{V}|$ (see Fig. 15 left), as a function of $m_{Z^{\prime}}$, where $c_{\ell A}^{11}$ is the axial coupling to electrons, while $c_{N}^{V}=c_{uV}^{11}(2Z+N)/{A}+c_{dV}^{11}(2N+Z)/A\simeq 1.41c_{uV}^{11}+1.59c_{dV}^{11}$ is the average vector coupling of $Z^{\prime}$ 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 Wood:1997zq . Translating the results of Ref. Dzuba:2017puc to our notation gives g^{\prime}{}^{2}c_{\ell A}^{11}c_{N}^{V}>3.9\cdot 10^{-8}\big{(}m_{Z^{\prime}}^{2}/{\rm GeV}^{2}\big{)}, which is not very stringent and is comparable to the nEDM bound in Fig. 14 (bottom right).

Constraints from SN1987a:

The $Z^{\prime}$ bosons can be copiously produced in a core of a Supernova (SN) if they are light enough. If $g^{\prime}$ is small enough the $Z^{\prime}$ 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^{\prime}$ the $Z^{\prime}$ bosons are not produced efficiently, while for larger $g^{\prime}$ the $Z^{\prime}$ are efficiently trapped inside the SN. We restrict the analysis to $m_{Z^{\prime}}\geq 10$ MeV, in which case two simplifications occur that make it easy to rescale reliably the results for dark photon from Ref. Chang:2016ntp to our case of a $Z^{\prime}$ with predominantly axial vector couplings. First of all, for $m_{Z^{\prime}}\geq 10$ MeV we can neglect the corrections due to the coupling of $Z^{\prime}$ with the electron plasma, an important effect for lighter gauge bosons. Furthermore, for $m_{Z^{\prime}}\geq 10$ MeV the main mechanism of production and absorption of $Z^{\prime}$ in the SN core is bremsstrahlung in neutron–proton scattering. The bounds can then be obtained by simply rescaling the results from Chang:2016ntp . We use that the $pn\to pnZ^{\prime}$ cross sections, $\sigma_{A,V}$, induced by the axial or vector couplings of $Z^{\prime}$ to proton, respectively, satisfy the numerical relation $\sigma_{A}/\sigma_{V}\sim 3c_{pA}^{2}/c_{pV}^{2}$, where $c_{pA,V}=2c_{uA,V}^{11}+c_{dA,V}^{11}$. We can thus translate the bounds in Chang:2016ntp by replacing $\epsilon^{\prime}$ with $\sqrt{3}g^{\prime}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, $\nu_{i}A\to\nu_{j}\ell_{k}^{+}\ell_{l}^{-}A$, can receive a contribution from a tree level exchange of an extra $Z^{\prime}$, see Fig. 15 (middle). This would result in a deviation of the total cross section from the SM prediction, signaling new physics Altmannshofer:2014pba . For small $Z^{\prime}$ 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. Altmannshofer:2019zhy for the $L_{\mu}-L_{\tau}$ gauge boson directly translate to our case (see also Ballett:2019xoj ). The bound on $g^{\prime}$ from Fig. 8 of Altmannshofer:2019zhy needs only to be re-interpreted as the bound on g^{\prime}\big{[}c_{\nu_{L}}^{22}\big{(}(-\frac{1}{4}+s_{W}^{2})c_{\ell_{L}}^{22}+s_{W}^{2}c_{\ell_{R}}^{22}\big{)}/\big{(}-\frac{1}{4}+2s_{W}^{2}\big{)}\big{]}^{1/2}, with $s_{W}\equiv\sin\theta_{W}$ the sine of the weak mixing angle. The resulting bounds from the CCFR experiment Mishra:1991bv and the future projections for DUNE Altmannshofer:2019zhy are shown in Fig. 14 (top right) as solid black and dashed black lines, respectively.

Electron-neutrino scattering experiments:

Measurements of $\nu_{i}-e^{-}$ and $\bar{\nu}_{i}-e^{-}$ scattering cross sections can bound the couplings $c_{\ell A}^{11}c_{\nu L}^{ii}$ and $c_{\ell V}^{11}c_{\nu L}^{ii}$, where $i$ is the flavor of the neutrino in the beam. The most stringent constraints are due to CHARM-II Vilain:1994qy , which used ${\cal O}(10~{}{\rm GeV})$ $\nu_{\mu}$ and $\bar{\nu}_{\mu}$ beams, and from TEXONO Deniz:2009mu , which used ${\cal O}(1~{}{\rm MeV})$ reactor $\bar{\nu}_{e}$ beam. We use the CHARM-II measured ratio ${\sigma(\nu_{\mu}e)}/{\sigma(\bar{\nu}_{\mu}e)}=0.910\pm 0.113$ Vilain:1994qy , where $\sigma(\nu_{\mu}e)$ and $\sigma(\bar{\nu}_{\mu}e)$ are the total $\nu_{\mu}e$ and $\bar{\nu}_{\mu}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 $\bar{\nu}_{e}$ scattering rate, normalized to the SM prediction, ${R_{\rm exp}}/{R_{\rm SM}}=1.08\pm 0.36$ Deniz:2009mu results in the bound on $g^{\prime}$ shown as a green solid line in Fig. 14 (top left). To derive these constraints we used the analytical results from Ref. Lindner:2018kjo .

Møller scattering:

Measurements of parity violating contributions to the Møller scattering, $e^{-}e^{-}\to e^{-}e^{-},$ bound the product $c_{\ell V}^{11}c_{\ell A}^{11}$. The most precise measurements were performed by E158 at SLAC at $q^{2}\simeq(0.16~{}{\rm GeV})^{2}$ Anthony:2005pm . Comparison with the SM gives for $m_{Z^{\prime}}\lesssim 100$ MeV the bound $g^{\prime}{}^{2}|c_{\ell V}^{11}c_{\ell A}^{11}|\lesssim 10^{-8}$, and for $m_{Z^{\prime}}\gtrsim 100$ MeV the bound $g^{\prime}{}^{2}|c_{\ell V}^{11}c_{\ell A}^{11}|\lesssim 3\cdot 10^{-8}\cdot(m_{Z^{\prime}}/200~{}{\rm MeV})^{2}$ (see also Kahn:2016vjr ). In our benchmark $c_{\ell V}^{11}c_{\ell A}^{11}\simeq 0.30$, giving relatively weak bounds $g^{\prime}\lesssim 2\cdot 10^{-4}$ and $g^{\prime}\lesssim 2\cdot 10^{-4}(m_{Z^{\prime}}/100~{}{\rm MeV})$, respectively, which we thus do not plot in Fig. 14.

Isotope shift spectroscopy:

The isotope shift spectroscopy constrains vector couplings of $Z^{\prime}$, giving $g^{\prime}{}^{2}|c_{\ell V}^{11}c_{fV}^{11}|\lesssim 10^{-7}\cdot{m_{Z^{\prime}}^{2}}/{(10~{}{\rm MeV})^{2}}$ for $f=u,d$, when charge radius determination from Lamb shift in muonic atoms is used Delaunay:2017dku (see also Delaunay:2016brc ; Frugiuele:2016rii ; Berengut:2017zuo ). Since the $Z^{\prime}$ has suppressed vector couplings, these bounds are not very constraining, and we do not consider them further.

White dwarf cooling:

The tree level $Z^{\prime}$ exchange contributes to the $e^{+}e^{-}\to\nu\bar{\nu}$ process which can increase the cooling rate of the white dwarf (WD) core Dreiner:2013tja . 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^{\prime}$ benchmarks the contribution to the star cooling is described by the effective Lagrangian

[TABLE]

since the $Z^{\prime}$ is much heavier than the WD internal temperature of a few keV, and can be integrated out. Translating the limits from Dreiner:2013tja to our notation gives (see also Bauer:2018onh ),

[TABLE]

with $(c_{\nu_{L}}^{\rm eff})^{2}=\sum_{ij}|c_{\nu_{L}}^{ij}|^{2}$, and $(c_{e}^{\rm eff})^{2}=|c_{\ell A}^{11}|^{2}+|c_{\ell V}^{11}|^{2}$. In our benchmark, $c_{\nu L}^{\rm eff}\simeq 3.4$, and $c_{e}^{\rm eff}\simeq 2.5$, giving the exclusion shown in Fig. 14 (top left) as a purple region.

Anomalous magnetic moments:

At one loop the $Z^{\prime}$ exchange contributes to the lepton anomalous magnetic moment, $(g-2)_{\ell}$, cf. Fig. 17 (with $\ell_{i}=\ell_{j}$). For our $Z^{\prime}$ benchmark we only need to keep the contribution to $(g-2)_{\mu}$ from the diagonal couplings, with $\mu$ running in the loop. A similar diagram with a $\tau$ running in the loop does get a chirality flip enhancement of $m_{\tau}/m_{\mu}\sim 10$, but is also suppressed by two off-diagonal couplings, $|(c_{\ell_{L}}^{23}c_{\ell_{R}}^{23})/(c_{\ell_{L}}^{22})^{2}|\sim 10^{-3}$. For the $(g-2)_{e}$ the two corresponding factors are $m_{\tau}/m_{e}\sim 3\cdot 10^{3}$ and $|(c_{\ell_{L}}^{13}c_{\ell_{R}}^{13})/(c_{\ell_{L}}^{11})^{2}|\sim 6\cdot 10^{-5}$, so that again the diagram with diagonal couplings dominates. Using the results of Ref. Fayet:2007ua with a trivial change of notation gives (see also Kozaczuk:2016nma ),

[TABLE]

with $i=1,2$ for $\ell=e,\mu$, respectively. The normalization of the loop functions,

[TABLE]

is such that for heavy $Z^{\prime}$, $\lim_{u\to\infty}F(u)=1$, $\lim_{u\to\infty}H(u)=5/3$, whereas for $m_{Z^{\prime}}=m_{\ell}$, $F(1)\simeq 0.31,H(1)\simeq 1.31$, and for light $Z^{\prime}$, i.e., for $u\ll 1$, $F(u)\simeq 3u^{2}/2$, $H(u)\simeq 1$.

In Fig. 14 (top right) we show in gray the $1\sigma$ band in the $g^{\prime},m_{Z^{\prime}}$ parameter space that gives $(g-2)_{\mu}$ in agreement with experiment. The required value of $g^{\prime}$ 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^{\prime}$ 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^{\prime}$ 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 Hollik:1998wk

[TABLE]

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 Beenakker:1991ca . The chromo-EDMs, $\tilde{d}_{f_{i}}$, are obtained from (58) by replacing $g\to g_{s}$. The experimental bound on electron EDM, $|d_{e}|<1.1\cdot 10^{-29}e\,{\rm cm}$ Andreev:2018ayy , 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, $\tilde{d}_{q}$, so that d_{n}=\sum_{q=u,d,s}\big{(}\beta_{q}d_{q}+\tilde{\beta}_{q}\tilde{d}_{q}\big{)}. For matrix elements of quark EDM operators, $\beta_{q}$, a recent lattice QCD calculation obtained $\beta_{u}=0.784(30)$, $\beta_{d}=-0.204(15)$, $\beta_{s}=-0.0027(16)$ at $\mu=2$ GeV in $\overline{\rm MS}$ scheme Gupta:2018lvp . In contrast, the matrix elements of chromo-EDM operators are poorly known with the estimates ranging over an order of magnitude, $\tilde{\beta}_{u}=-(0.09-0.9)$, $\tilde{\beta}_{d}=-(0.2-1.8)$ Engel:2013lsa . In Fig. 14 (bottom right) we used the estimated best values, $\tilde{\beta}_{u}=-0.35$, $\tilde{\beta}_{d}=-0.7$ Engel:2013lsa .

5.1.2 Flavor violating couplings to quarks

The exchange of $Z^{\prime}$ 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}-\bar{K}^{0}$ mixing:

To estimate the contributions of tree level $Z^{\prime}$ exchanges to $K-\bar{K}$ mixing we distinguish two regimes. For light $Z^{\prime}$, $m_{Z^{\prime}}\lesssim m_{K}$, we can use ChPT with the $Z^{\prime}$ acting as an external field. Using ChPT expansion we have $(\bar{s}\gamma^{\mu}\gamma_{5}d)\to-f_{K}\partial^{\mu}K^{0}+\cdots$, where the ellipses denote higher orders in momentum expansion. Using this in Eq. (44) gives for the tree level $Z^{\prime}$ exchange contribution to $K-\bar{K}$ mixing (cf. Eq. (46))

[TABLE]

where $f_{K}\simeq 156$ MeV is the kaon decay constant, and we use the same phase conventions for $M_{12}$ as in Buras:1998raa ; Buras:2010pza . Note that this contribution is proportional to $1/m_{Z^{\prime}}^{2}$ even though the momentum flowing in the $Z^{\prime}$ propagator is ${\mathcal{O}}(m_{K})$, so that $M_{12}^{Z^{\prime}}\propto 1/\langle\phi\rangle^{2}$. This is most easily understood in the Feynman gauge, where the dominant contribution to $M_{12}^{Z^{\prime}}$ due to the exchange of the $\arg(\phi)$ Goldstone boson, whose couplings to SM fermions are $\propto 1/\langle\phi\rangle$ and do not depend on $g^{\prime}$. In the unitary gauge the dominant contribution arises from the longitudinal component of the $Z^{\prime}$ propagator. The vector couplings of $Z^{\prime}$ contribute only at higher order in ChPT and are relatively suppressed by ${\mathcal{O}}(m_{K}^{2}/(4\pi f_{K})^{2})\sim{\mathcal{O}}(0.2)$, on top of the suppression of vector couplings themselves. For completeness we calculate these contributions in Appendix A.

For $m_{Z^{\prime}}\gg m_{K}$ we can integrate out the $Z^{\prime}$ at the scale $\mu\simeq m_{Z^{\prime}}$ and match onto the effective weak Hamiltonian, $H_{\rm eff}=\sum_{a}C_{a}Q_{a}$, where the sum runs over the operators in Ref. Bona:2007vi . The nonzero NP contributions to the Wilson coefficients are,

[TABLE]

where $Q_{1}^{sd}=(\bar{d}\gamma^{\mu}s_{L})^{2}$, $\tilde{Q}_{1}^{sd}=(\bar{d}\gamma^{\mu}s_{R})^{2}$, $Q_{5}^{sd}=(\bar{d}^{\alpha}s_{R}^{\beta})(\bar{d}^{\beta}s_{L}^{\alpha})$. To compare with the experimental results we run from $\mu\simeq m_{Z^{\prime}}$ down to $2$ GeV using the results of Bona:2007vi ; Ciuchini:1998ix .

The comparison of the SM predictions and the experimental measurements for indirect CP violation parameter $\epsilon_{K}$ is made through a ratio

[TABLE]

Estimating the errors due to charm loop contributions in the SM prediction of $\epsilon_{K}$ is nontrivial and leads to differences between CKMfitter and UTFit collaborations Charles:2015gya ; UTFit (see also the introduction in Ligeti:2016qpi ). For consistency we use the UTFit collaboration extraction $0.87<C_{\epsilon_{K}}<1.39$ at 95% C.L. UTFit as well as their prediction for $\epsilon_{K}^{\rm SM}$ using NLO charm contribution Ciuchini:2000de . We denote the resulting bounds in Fig. 13 and 14 (bottom right) with black lines, where we use the heavy $Z^{\prime}$ solution (60) for $m_{Z^{\prime}}>3m_{K}$, the light $Z^{\prime}$ solution (59) for $m_{Z^{\prime}}<m_{K}/3$, and connect the two with naive linear interpolation to guide the eye. The $Z^{\prime}$ contribution to the $K-\bar{K}$ mixing is $\propto 1/\langle\phi\rangle^{2}$. For heavy $Z^{\prime}$, $m_{Z^{\prime}}\gg m_{K}$, this then leads to a bound $\langle\phi\rangle=m_{Z^{\prime}}/g^{\prime}>1.2\cdot 10^{7}$ GeV for our benchmark.

$B_{q}^{0}-\bar{B}_{q}^{0}$ mixing:

We distinguish two limits, the heavy and light $Z^{\prime}$. For heavy $Z^{\prime}$, $m_{Z^{\prime}}\gg m_{B_{q}}$, the $Z^{\prime}$ is integrated out at $\mu\simeq m_{Z^{\prime}}$. The matching onto $H_{\rm eff}=\sum_{a}C_{a}Q_{a}$, gives, analogously to (60),

[TABLE]

where $q=d,s$, with $i=1,2$, respectively. We RG evolve the above Wilson coefficients from $\mu\simeq m_{Z^{\prime}}$ down to $\mu=4.2{\rm~{}GeV}\simeq m_{b}$ using the results of Bona:2007vi ; Ciuchini:1998ix , and compare with the allowed deviations in the mixing matrix elements,

[TABLE]

The latest fit results from UTFit collaboration give $0.942<C_{B_{s}}<1.288$, $-1.35^{\circ}<\phi_{B_{s}}<2.21^{\circ}$, and $0.83<C_{B_{d}}<1.29$, $-6.0^{\circ}<\phi_{B_{d}}<1.5^{\circ}$ UTFit .

For light $Z^{\prime}$ we perform an operator product expansion, where, to leading order in $\Lambda_{\rm QCD}/m_{b}$, the dominant Wilson coefficients of the $B_{q}^{0}-\bar{B}_{q}^{0}$ operators are at $\mu\simeq m_{b}$ given by

[TABLE]

with $m_{b}\simeq 4.2$ GeV the $b$ quark mass, while the remaining nonzero Wilson coefficients, $C_{2}^{qb}=-C_{1}^{qb}m_{Z^{\prime}}^{2}/m_{b}^{2}$, $\tilde{C}_{2}^{qb}=-\tilde{C}_{1}^{qb}m_{Z^{\prime}}^{2}/m_{b}^{2}$, $C_{5}^{qb}=2C_{4}^{qb}m_{Z^{\prime}}^{2}/m_{b}^{2}$, are parametrically suppressed.

We denote the resulting bound in Fig. 14 (bottom right) with a green (brown) solid line for $B_{s}-\bar{B}_{s}$ ($B_{d}-\bar{B}_{d}$) mixing, while the corresponding dashed lines shows the projected bounds after Belle II and LHCb Upgrade II Cerri:2018ypt .

Charm mixing:

The $Z^{\prime}$ contributions to $D-\bar{D}$ mixing take the same form as for $B_{q}^{0}-\bar{B}_{q}^{0}$, after making the replacements $i\to 1$, $3\to 2$ and $d,q\to u,b\to c$ in Eqs. (62), (64). The expressions for the case of light $Z^{\prime}$ are only approximate, valid to the extent that one can use the operator product expansion despite the relatively light charm quark mass. The resulting present bound UTFit and projected sensitivity Cerri:2018ypt are shown in Fig. 14 (bottom right) as solid and dashed red lines, respectively.

5.1.3 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.

$\mu\to e$ conversion:

Stopped muons can undergo coherent $\mu\to e$ conversion in the field of the nucleus. The most stringent limits to date were obtained by the SINDRUM II experiment using 197Au Bertl:2006up

[TABLE]

Here $\Gamma(\mu^{-}\text{Au}\to e^{-}\text{Au})$ is the conversion rate, and $\Gamma_{\rm capt}(\mu^{-}\text{Au})$ the muon capture rate. State of the art estimates for both of these can be found in Ref. Kitano:2002mt for a whole range of nuclei, including gold. The tree level $Z^{\prime}$ exchange induces four-fermion effective interaction for $\mu\to e$ conversion, see Fig. 17, which in the notation of Refs. Kitano:2002mt ; Kuno:1999jp takes the form {\cal L}_{\rm int}=-\big{[}(G_{F}/\sqrt{2})\sum_{f=u,d}g_{LV(f)}(\bar{e}\gamma_{\mu}P_{L}\mu)(\bar{f}\gamma^{\mu}f)+g_{RV(f)}(\bar{e}\gamma_{\mu}P_{R}\mu)(\bar{f}\gamma^{\mu}f)\big{]}+\big{[}(\bar{f}\gamma^{\mu}f)\to(\bar{f}\gamma^{\mu}\gamma_{5}f),V\to A\big{]}, with

[TABLE]

and similarly for $g_{RV(f)}$, but with $L\to R$ replacement, while $g_{LA(f)}$ and $g_{RA(f)}$ trivially follow from $V\to A$ replacement.

The $g_{RV(f)}$ and $g_{LV(f)}$ operators induce coherent $\mu\to e$ conversion for which we use the results of Ref. Kitano:2002mt . The $g_{RA(f)}$ and $g_{LA(f)}$ operators induce spin dependent $\mu\to e$ conversion, for which we use the results of Cirigliano:2017azj , where the Wilson coefficients in the notation of Cirigliano:2017azj are $C_{A,L}^{ff}=g_{LA(f)}/4$, $C_{A,L}^{ff}=g_{LA(f)}/4$, while the axial structure factors $S_{A}$ we take from Engel:1995gw . Note that the momentum exchange in $\mu\to e$ conversion is $q^{2}\simeq-m_{\mu}^{2}$. This is small enough that it is easily absorbed by the nucleus without changing its nuclear structure. On the other hand, $q^{2}\simeq-m_{\mu}^{2}$ is also large enough that even for very light $Z^{\prime}$, $m_{Z^{\prime}}\ll m_{\mu}$, 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}\sim 4/m_{\mu}$. This means that we can still use the calculations of overlap integrals from Ref. Kitano:2002mt even for light $Z^{\prime}$, up to corrections of about ${\mathcal{O}}(30\%)$.

Spin-dependent conversion starts to be dominant at very low masses. While the spin-independent conversion rate is coherently enhanced and at low masses goes as $\sim A^{2}m_{\mu}$, where $A$ is the atomic number, the spin-dependent rate increases as $m_{\mu}^{3}/m_{Z^{\prime}}^{2}$. For heavy nuclei, such as 197Au used by the SINDRUM II experiment, the spin-dependent contribution starts to dominate around $m_{Z^{\prime}}\sim 0.5$ MeV, while for lighter nuclei, such as 26Al used by Mu2e, it begins at $m_{Z^{\prime}}\sim 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^{\prime}$ parameter space.

The constraint from Eq. (65) 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 ${\rm Br}(\mu\to e)<8\cdot 10^{-17}$ for $\mu\to e$ conversion on Al nuclei Bernstein:2019fyh , which gives the projected sensitivity on $g^{\prime}$ denoted with a dashed green line in Figs. 13 and 14 (top right).

Decays to three leptons:

The tree level $Z^{\prime}$ exchanges mediate the lepton flavor violating decays such at $\tau\to 3\mu,e2\mu,\ldots$, see Fig. 17. To shorten the discussion we focus on the decays of the form $\ell_{i}\to 3\ell_{j}$: $\tau\to 3\mu$, $\tau\to 3e$ and $\mu\to 3e$. Similar constraints are obtained also from the remaining modes, $\tau^{-}\to\mu^{-}e^{+}e^{-},e^{-}\mu^{+}e^{-},\mu^{-}\mu^{+}e^{-},\mu^{-}e^{+}\mu^{-}$, etc.

Depending on the value of $m_{Z^{\prime}}$ there are three different regimes for $\Gamma(\ell_{i}\to 3\ell_{j})$. For $m_{Z^{\prime}}\gg m_{\ell_{i}}$, the $Z^{\prime}$ can be integrated out, leading to effective four fermion interaction. In this limit, the $\ell_{i}\to 3\ell_{j}$ decay is a genuine three body decay, with the decay width Heeck:2016xkh

[TABLE]

where we neglected the terms proportional to $m_{\ell_{j}}$.

In the intermediate mass regime, $2m_{\ell_{j}}<m_{Z^{\prime}}<m_{\ell_{i}}-m_{\ell_{j}}$, the $Z^{\prime}$ can be produced on shell. The $\ell_{i}\to 3\ell_{j}$ decay is thus a cascade decay, $\ell_{i}\to Z^{\prime}\ell_{j}$ followed by $Z^{\prime}\to\ell_{j}^{+}\ell_{j}^{-}$, so that

[TABLE]

where

[TABLE]

where we neglected the mass of $\ell_{j}$, while ${\rm Br}(Z^{\prime}\to\ell_{j}\ell_{j})$ is the branching ratio for the $Z^{\prime}\to\ell_{j}^{+}\ell_{j}^{-}$ decay. In our benchmark ${\rm Br}(Z^{\prime}\to\mu\mu)\sim 0.16$, while ${\rm Br}(Z^{\prime}\to ee)\sim 0.23$ for $m_{\mu}<m_{Z^{\prime}}<m_{\tau}$ and ${\rm Br}(Z^{\prime}\to ee)\sim 0.28$ for $m_{Z^{\prime}}<m_{\mu}$. Note that $\Gamma(\ell_{i}\to\ell_{j}Z^{\prime})\propto g^{\prime}{}^{2}/m_{Z^{\prime}}^{2}=1/\langle\phi\rangle^{2}$ and thus does not diverge for $m_{Z^{\prime}}\to 0$.

Finally, if $m_{Z^{\prime}}<m_{\ell_{j}}$ the $Z^{\prime}$ is off-shell and we again have a genuine three body decay. For $m_{Z^{\prime}}\ll m_{\ell_{j}}$ the axial current contributions are $1/m_{Z^{\prime}}^{4}$ enhanced, giving Heeck:2016xkh

[TABLE]

The contributions from flavor diagonal vector couplings of $Z^{\prime}$ are subleading and we can safely neglect them.

The above expressions constrain $g^{\prime}$ as a function of $m_{Z^{\prime}}$ given present bounds on the $\tau\to 3\mu,3e$ and $\mu\to 3e$ decays. For our benchmark the most stringent is the $90\%$ C.L. bound ${\rm Br}(\mu\to 3e)<1.0\cdot 10^{-12}$ Bellgardt:1987du , shown as the purple solid line in Fig. 14 (top right), while the the purple dashed line corresponds to the projected reach by Mu3e of $5.2\cdot 10^{-15}$ Perrevoort:2018ttp . The present (projected Belle II) bounds ${\rm Br}(\tau\to 3\mu)<2.1(0.03)\cdot 10^{-8}$ and ${\rm Br}(\tau\to 3e)<2.7(0.04)\cdot 10^{-8}$ Hayasaka:2010np ; Kou:2018nap are shown as blue and purple solid (dashed) lines in Fig. 14 (bottom left). The different thresholds appear as the breaks in the slopes of the lines.

Rare meson decays (RMD):

The tree level $Z^{\prime}$ exchanges induce rare FCNC meson decays such as $M_{1}\to M_{2}\ell_{i}^{+}\ell_{j}^{-}$, where $M_{1}(M_{2})$ is the initial (final) state meson. An example of such an exotic decay is $K^{+}\to\pi^{+}\mu^{+}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^{\prime}$ exchange instead has no suppression due to neutrino masses. For heavy $Z^{\prime}$ 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^{\prime}$ is light enough to be produced on-shell, the $M_{1}\to M_{2}Z^{\prime}$ 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}\to M_{2}$ transitions the strongest experimental bounds follow from the $M_{1}\to M_{2}Z^{\prime}$ decays with $M_{1,2}$ both pseudoscalars, for which the matrix element is

[TABLE]

with $J^{\mu}$ defined in (44), and $p_{1,2}$ the momenta of mesons $M_{1,2}$, such that $p_{1,2}^{2}=m_{1,2}^{2}$, with $m_{1,2}$ the corresponding meson masses. The quark flavor-violating coupling $c_{qV}^{kl}$ encodes the strength of the $q_{l}\to q_{k}$ transition, e.g., the $c_{dV}^{23}$ coupling is obtained for the $B\to K$ and $c_{dV}^{13}$ for the $B\to\pi$ matrix elements, respectively. The form factors $f_{+,-}(q^{2})$ depend on the momentum exchange squared, $q^{2}$, with $q_{\mu}=(p_{1}-p_{2})_{\mu}$. We follow the common practice and trade the form factor $f_{-}(q^{2})$ for

[TABLE]

We use the $z$-expansion based parametrization of the form factors Bourrely:2008za ,

[TABLE]

where

[TABLE]

with $t_{+}=(m_{1}+m_{2})^{2}$ , $t_{0}=(m_{1}+m_{2})(\sqrt{m}_{1}-\sqrt{m}_{2})^{2}$.

For light $Z^{\prime}$, with a mass below the lepton pair production threshold, $m_{Z^{\prime}}\leq m_{\ell_{i}}+m_{\ell_{j}}$, the $M_{1}\to M_{2}\ell_{i}^{-}\ell_{j}^{+}$ transition is a genuine 3-body decay, mediated by an off-shell $Z^{\prime}$. The longitudinal $Z^{\prime}$ contribution dominates, leading to a decay width proportional to $\propto g^{\prime}{}^{4}/m_{Z^{\prime}}^{4}\simeq 1/(4\langle\phi\rangle^{4})$. The expression for the differential width is well approximated by taking the $m_{Z^{\prime}}\to 0$ limit, giving

[TABLE]

where $m_{\ell_{i}\ell_{j}}^{2}=(p_{\ell_{i}}+p_{\ell_{j}})^{2}$ is the invariant mass squared of the leptonic pair, and similarly $m_{M_{1}\ell_{i}}^{2}=(p_{1}+p_{\ell_{i}})^{2}$.

For heavy $Z^{\prime}$, with a mass $m_{Z^{\prime}}\gg m_{1}$, the $Z^{\prime}$ can be integrated out. This also gives a decay width proportional to $\propto g^{\prime}{}^{4}/m_{Z^{\prime}}^{4}\simeq 1/(4\langle\phi\rangle^{4})$. In the limit of massless final states, $m_{2},m_{\ell_{i}},m_{\ell_{j}}\to 0$, the differential rate is given by a rather compact expression,

[TABLE]

In the numerical analysis we work with the full dependence on $m_{2},m_{\ell_{i}},m_{\ell_{j}}$.

Phenomenologically, the most relevant is the intermediate mass regime, $m_{\ell_{i}}+m_{\ell_{j}}<m_{Z}^{\prime}<m_{1}-m_{2}$. In this case the $Z^{\prime}$ is produced on-shell in the 2-body decay, $M_{1}\to M_{2}Z^{\prime}$, followed by a decay into two leptons, $Z^{\prime}\to\ell_{i}\bar{\ell}_{j}$. The total decay width is given by

[TABLE]

The two-body decay width for $M_{1}\to M_{2}Z^{\prime}$ transition depends only on the quark current matrix element, and is given by

[TABLE]

where $y_{i}=m_{i}/m_{1}$. The branching fractions for $Z^{\prime}\to\ell_{i}\bar{\ell}_{j}$ are easily calculated. Ignoring the kinematical factors for final state particles they are given by ${\rm Br}(Z^{\prime}\to\ell_{i}\bar{\ell}_{j})=(|c_{\ell V}^{ij}|^{2}+|c_{\ell A}^{ij}|^{2})/(\sum_{ff^{\prime}}|c_{\ell V}^{ff^{\prime}}|^{2}+|c_{\ell A}^{ff^{\prime}}|^{2})$. For our benchmark we have ${\rm Br}(Z^{\prime}\to ee,\mu\mu)\simeq 0.1$ and ${\rm Br}(Z^{\prime}\to\sum_{\alpha}\nu_{\alpha}\nu_{\alpha})\simeq 0.4$ for the flavor conserving channels, and ${\rm Br}(Z^{\prime}\to e\mu)\sim 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 10.1093/ptep/ptaa104 , while the numerical inputs for the form factors are taken from Ref. Bondarenko_2018 . The envelope of constraints on $g^{\prime}$ as a function of $Z^{\prime}$ 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^{\prime}$ masses, i.e., as long as the $B\to\pi Z^{\prime}$ 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^{\prime}$ masses the $K^{+}\to\pi^{+}Z^{\prime},~{}Z^{\prime}\to\nu\nu$ decay leads to the strongest constraints, a result of $K\to\pi+{\rm inv}$ being well constrained experimentally, as well as the large $Z^{\prime}\to\nu\bar{\nu}$ branching ratio. Conservatively, we use the integrated branching ratio ${\rm Br}(K^{+}\to\pi^{+}\nu\nu)<2.3\times 10^{-10}\,$ to set the limit on $g^{\prime}$ 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.

The bounds on $K^{+}\to\pi^{+}\mu^{+}e^{-}$ decays, ${\rm Br}(K^{+}\to\pi^{+}\mu^{+}e^{-})<1.3\times 10^{-11}\,,{\rm Br}(K^{+}\to\pi^{+}\mu^{-}e^{+})<6.6\times 10^{-11}$ lead to the most stringent constraints on $g^{\prime}$ in the intermediate $Z^{\prime}$ mass region, despite the small ${\rm Br}(Z^{\prime}\to e\mu)$. The related decays, $K_{L}\to\mu e$, with the decay width given by

[TABLE]

are always less constraining. When $K\to\pi Z^{\prime}$ decays are kinematically forbidden the decays of heavier mesons become relevant. In Figs. 13 and 14. we show the limits that follow from the bounds on the decays, ${\rm Br}(D^{+}\to\pi^{+}\mu\mu)<7.3\times 10^{-8}\,,{\rm Br}(B^{+}\to\pi^{+}\mu\mu)<6.9\times 10^{-8}\,,$ and ${\rm Br}(B^{+}\to K^{+}\mu\mu)<4.8\times 10^{-7}\,$.

Flavor violating radiative decays:

At one loop the off diagonal couplings of $Z^{\prime}$ to leptons generate the $\ell_{i}\to\ell_{j}\gamma$ transition with the decay width Lavoura:2003xp , see also Fig. 17,

[TABLE]

where $x_{a/b}=m_{a}/m_{b}$, and

[TABLE]

with $Q_{k}$ the charge of the fermion $\ell_{k}$ running in the loop, while $c_{R}^{\gamma}$ can be obtained by making the replacement $L\leftrightarrow R$ in the above expression. The loop functions are

[TABLE]

where (below we shorten $x_{k/Z^{\prime}}\to x$)

[TABLE]

The present bounds on ${\rm Br}(\tau\to\mu(e)\gamma)<4.4(3.3)\cdot 10^{-8}$ Aubert:2009ag and ${\rm Br}(\mu\to e\gamma)<4.2\cdot 10^{-13}$ TheMEG:2016wtm 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 Kou:2018nap and MEG-II Baldini:2018nnn ${\rm Br}(\tau\to\mu(e)\gamma)<1(3)\cdot 10^{-9}\,,{\rm Br}(\mu\to e\gamma)<6\cdot 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 $\tau$ 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.

5.2 Bounds on the $Z^{\prime}_{i}$ in the $G_{\rm FN}=U(1)_{\rm FN}^{3}$ model

Next, we extend the above analysis to the $G_{\rm FN}=U(1)_{\rm FN}^{3}$ benchmark, with the couplings of the three gauge bosons, $Z_{i}^{\prime}$, $i=1,2,3$, listed in Appendix B.1. The $Z_{3}^{\prime}$ only couples to $d_{R}$ and $\ell_{L},\nu_{L}$, see (130)-(132). This follows immediately from the charge assignments in Eqs. (15) and (29), because $q_{L}^{(i)}$ and $u_{R}^{(i)}$ are not charged under $U(1)_{3}$. The large mixings in $d_{R}$ and lepton sectors result in an almost anarchic form of $Z_{3}^{\prime}$ couplings to $d_{R}$ and $\ell_{L},\nu_{L}$, with large couplings to the 2nd and 3rd generation and only slightly suppressed couplings to the 1st generation SM fermions. The $Z_{1}^{\prime}$ and $Z_{2}^{\prime}$, on the other hand, predominantly couple to the first and second generation, respectively, with the exception of $Z_{2}^{\prime}$ couplings to $d_{R}$ and $\ell_{L}$, $\nu_{L}$, where there are also large couplings to the third generation, again due to the large mixings, see Eqs. (122)-(129).

For all three $Z^{\prime}_{i}$ the most stringent constraints arise from rare meson decays, $\mu\to e$ conversion, $K-\bar{K}$ mixing and SN constraints. This is the same as for the $U(1)_{\rm FN}$, Section 5.1. However, the relative importance of these constrains has shifted because of the changes in the flavor patterns of the $Z^{\prime}_{i}$ couplings, see Fig. 18.

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^{\prime}_{2}$ and $Z^{\prime}_{3}$ the 12 off-diagonal couplings can dominate over the diagonal 11 couplings, see, e.g., Eq. (130). This means that we need to adjust our treatment of bounds from direct production of $Z^{\prime}_{i}$ compared to what was done in Section 5.1. For $m_{e}+m_{\mu}<m_{Z^{\prime}}<2m_{\mu}$ the $Z^{\prime}_{i}\to e\mu$ decay is now the dominant decay channel. In the $G_{\rm FN}=U(1)_{\rm FN}^{3}$ benchmark we have $\Gamma(Z^{\prime}_{2}\to ee)/\Gamma(Z^{\prime}_{2}\to e\mu)\simeq 0.16$ and $\Gamma(Z^{\prime}_{3}\to ee)/\Gamma(Z^{\prime}_{3}\to e\mu)\simeq 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}^{\prime}$ than they are for $Z_{2}^{\prime}$ and $Z_{3}^{\prime}$, because $Z_{2,3}^{\prime}$ 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^{\prime}$ for $Z_{2}^{\prime}$ and $Z_{3}^{\prime}$.

Fig. 18 also gives the FCNC constraints from rare meson decays, $D-\bar{D}$, $B_{d}-\bar{B}_{d}$ and $B_{s}-\bar{B}_{s}$ mixing (blue, red, brown, cyan lines). The expressions for these can be taken directly from Section 5.1, with trivial relabeling corresponding to $Z^{\prime}\to Z_{i}^{\prime}$ replacements. For $Z_{1}^{\prime}$ the most stringent constraints come from the FCNCs involving the first two generations, rare kaon decays, $\mu\to e$ conversion and $K-\bar{K}$ mixing, since the $Z_{1}^{\prime}$ has appreciable 12 off-diagonal couplings. The $D-\bar{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}-\bar{B}_{d}$ mixing. The $Z_{2}^{\prime}$ and $Z_{3}^{\prime}$ have much smaller 12 off-diagonal couplings. Even so, rare meson decays, $\mu\to e$ conversion and $K-\bar{K}$ mixing still lead to the most stringent constraints in our benchmark, but the FCNCs involving the third generation, $B_{d}-\bar{B}_{d}$ and $B_{s}-\bar{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_{\rm FN}=U(1)_{\rm FN}^{3}$ benchmarks to lead to a similar pattern of constraints. However, the relative strengths of different constraints may easily differ by relative ${\mathcal{O}}(1)$ factors (for instance for some benchmarks we would expect $\mu\to e$ conversion to be more stringent than $K-\bar{K}$ mixing even for heavy $Z_{1}^{\prime}$, unlike what was found here). Fig. 18 should thus be taken only as an illustration of how important different probes are for $G_{\rm FN}=U(1)_{\rm FN}^{3}$ inverted FN models.

6 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/\langle\phi\rangle$, where $M$ is the vector-like mass parameter and $\langle\phi\rangle$ the vev of the flavon field. The observed pattern of masses and mixings is obtained for $\langle\phi\rangle\gg M$, while in the traditional FN models the opposite is required, $M\gg\langle\phi\rangle$. Different realizations of the inverted FN models differ in the choice of the FN horizontal group $G_{\rm 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_{\rm 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_{\rm 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_{\rm 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_{\rm FN}$: the "decoupled chains" model, where $G_{\rm FN}=U(1)_{\rm FN}^{3}$, and the "coupled chains" model, with $G_{\rm FN}=U(1)_{\rm 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 ${\cal 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 $\langle\phi\rangle$. 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\cdot 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., Alonso:2018bcg . In the two models we presented there is a QCD Landau pole about four orders of magnitude above $\langle\phi\rangle$, i.e., for $\langle\phi\rangle\sim 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 $\langle\phi\rangle/M$ expansion parameter. Alternatively, for $\langle\phi\rangle\sim 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 Sannino:2019sch . The possible Landau pole in the $U(1)_{\rm FN}$ is less constraining; for $g^{\prime}<0.11$ at $\mu\simeq 10^{7}$ GeV the Landau pole in $g^{\prime}$ 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 $\langle\phi\rangle$ and $M_{a}$. We focused on the case where $G_{\rm FN}$ is gauged, resulting in three (one) gauge bosons $Z_{i}^{\prime}$ ($Z^{\prime}$) for $G_{\rm FN}=U(1)_{\rm FN}^{3}(U(1)_{\rm FN})$. Since the FN charges have opposite signs for the weak $SU(2)_{L}$ doublets and the weak $SU(2)_{L}$ singlets, the $Z_{i}^{\prime}$ have predominantly axial vector couplings to the SM fermions. The non-universal $G_{\rm FN}$ charges induce both flavor conserving and flavor violating couplings of $Z_{i}^{\prime}$ to the SM fermions. The most important FCNC contributions then come from tree level $Z_{i}^{\prime}$ exchanges $K-\bar{K}$ mixing and $\mu\to e$ conversion bounds limiting $\langle\phi\rangle\gtrsim 10^{7}$ GeV, and consequently also the masses of vector-like fermions to roughly the same scale. The $Z_{i}^{\prime}$, on the other hand, can be light, if $g^{\prime}$ is small. For light $Z^{\prime}$, with masses from several 10s MeV to several 100s MeV, the rare meson decays and beam dump searches are the most constraining, see Figs. 13, 14 and 18 for the bounds on two representative benchmarks.

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^{\prime}$ 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\to s\ell^{+}\ell^{-}$ transitions, potentially along the lines of Ref. Grinstein:2018fgb .

Acknowledgements.

We thank Ben Allanach, Stefan Antusch, Fady Bishara, Giancarlo d’Ambrosio, Gino Isidori, Jernej Kamenik, Nejc Košnik, and Uli Nierste for insightful discussions, and Joachim Brod, Alex Kagan and Sourov Roy for useful comments on the manuscript. MT and JZ acknowledge support in part by the DOE grant DE-SC0011784. JZ and AS thank American Slovenian Education Foundation for financial support. AS is supported by the Young Researchers Programme of the Slovenian Research Agency under the grant No. 50510. JZ thanks the theory group of the Jozef Stefan Institute for the hospitality and acknowledges the financial support from the Slovenian Research Agency (research core funding No. P1-0035). JZ acknowledges financial support from LAPTh, where part of the research was performed, and especially for their hospitality.

Appendix A Light $Z^{\prime}$ contributions to $K-\bar{K}$ mixing

In this appendix we give further details on the contributions from light $Z^{\prime}$ exchanges to the $K-\bar{K}$ mixing. For $m_{Z^{\prime}}\ll\Lambda_{\rm ChPT}\simeq 4\pi f_{\pi}\sim{\mathcal{O}}(1{\rm~{}GeV})$ we can use Chiral Perturbation Theory (ChPT) to describe these interactions. We estimate the leading contributions from both axial and vector couplings of $Z^{\prime}$. The former start at tree level, the latter at one loop.

In order to construct the appropriate ChPT Lagrangian in the presence of flavor violating $Z^{\prime}$ we use the spurion analysis Gasser:1984gg . The QCD Lagrangian is now ${\cal L}_{{\rm QCD}+Z^{\prime}}=\bar{q}(i\not{\partial}+g_{s}\not{G}^{a}T^{a}+\chi_{V}\not{Z}^{\prime}+\chi_{A}\not{Z}^{\prime}\gamma_{5})q-\bar{q}{\cal M}_{q}q$, where $q=(u,d,s)$, ${\cal M}_{q}=\operatorname{diag}(m_{u},m_{d},m_{s})$, and $\chi_{V,A}$ are $3\times 3$ Hermitian matrices of the form

[TABLE]

with $c_{d_{V,A}}^{ij}$ defined in (46). The QCD$+Z^{\prime}$ Lagrangian is formally invariant under $SU(3)_{R}\times SU(3)_{L}$ transformations, $q_{R,L}\to g_{R,L}(x)q_{R,L}$, if $\chi_{V,A}Z_{\mu}$, and ${\cal M}_{q}$ are promoted to spurions that transform as $v_{\mu}+a_{\mu}\to g_{R}(v_{\mu}+a_{\mu})g_{R}^{\dagger}+ig_{R}\partial_{\mu}g_{R}^{\dagger}$, $v_{\mu}-a_{\mu}\to g_{L}(v_{\mu}-a_{\mu})g_{L}^{\dagger}+ig_{L}\partial_{\mu}g_{L}^{\dagger}$, $s+ip\to g_{R}(s+ip)g_{L}^{\dagger}$, and that take the values $v^{\mu}=\chi_{V}Z_{\mu}^{\prime}$, $a^{\mu}=\chi_{A}Z_{\mu}^{\prime}$, $s={\cal M}_{q}$, $p=0$.

The LO ChPT Lagrangian, including $Z^{\prime}$ as the light degree of freedom, is therefore (we work in the unitary gauge for the $Z^{\prime}$ gauge boson)

[TABLE]

where $U=\exp(i\lambda^{a}\pi^{a}/f)$, with $\lambda^{a}$ the Gell-Mann matrices, $f$ is related to the meson decay constant,555We use the normalization $\langle 0|\bar{q}_{1}\gamma_{\mu}q_{2}|P(p)\rangle=ip_{\mu}f_{P}$. $f\simeq f_{\pi}/\sqrt{2}=93$ MeV, the covariant derivative is $\nabla_{\mu}U=\partial_{\mu}U-i(\chi_{V}+\chi_{A})Z_{\mu}^{\prime}U+iU(\chi_{V}-\chi_{A})Z_{\mu}^{\prime}$, while $B_{0}$ is the low energy constant. Under chiral rotations the $U$ field transforms as $U\to g_{R}Ug_{L}^{\dagger}$, $\nabla_{\mu}U\to g_{R}\nabla_{\mu}Ug_{L}^{\dagger}$.

Since we consider radiative corrections due to the $Z^{\prime}$ gauge bosons, the $Z^{\prime}$ 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}=\chi_{V}\pm\chi_{A}$, that transform as $Q_{R,L}\to g_{R,L}(x)Q_{R,L}\,g_{R,L}^{\dagger}(x)$. This gives the following contributions to the LO ChPT Lagrangian,

[TABLE]

where we only show the terms that will be relevant for $K-\bar{K}$ mixing, where they will serve as counter-terms to one-loop corrections. For this we also need part of the ${\mathcal{O}}(p^{4})$ Lagrangian,

[TABLE]

The ChPT Lagrangian is invariant under parity, $U\to U^{\dagger}$, $v^{\mu}\to-v^{\mu}$, $a^{\mu}\to a^{\mu}$, $s+ip\to s-ip$, $Q_{L}\leftrightarrow Q_{R}$. In addition, since the terms containing $Q_{R,L}$ can only arise from loops of $Z^{\prime}$, the ChPT terms without external $Z^{\prime}$ fields need to have even powers of $g^{\prime}$. This is insured by requiring that the ChPT Lagrangian is invariant under the $Z_{2}$ transformation $Q_{L,R}\to-Q_{L,R}$. The chiral counting of the spurions is $Q_{R,L}\sim{\mathcal{O}}(p^{0})$, $v_{\mu},a_{\mu},\sim{\mathcal{O}}(p)$, $s+ip\sim{\mathcal{O}}(p^{2})$. In principle one could thus have an arbitrary number of $Q_{R,L}$ insertions in ${\cal L}_{Q_{R,L}}^{(2),(4)}$. However, we work only to ${\mathcal{O}}(g^{\prime}{}^{2})$ and thus merely keep the terms that are $\sim Q_{R,L}^{2}$. Note that the first term in (88), the analogue of the ChPT+QED Lagrangian Ecker:1988te ; Urech:1994hd , is in our counting ${\mathcal{O}}(p^{0})$. However, its coefficient is $C_{0}^{(2)}\sim{\mathcal{O}}(m_{Z^{\prime}}^{2})$ and we thus include it as part of ${\cal L}_{Q_{R,L}}^{(2)}$.

In the usual ChPT the one-loop counterterms reside in the ${\mathcal{O}}(p^{4})$ Lagrangian, while in the light $Z^{\prime}$ case we need both the ${\cal L}_{Q_{R,L}}^{(2)}$ and ${\cal L}_{Q_{R,L}}^{(4)}$ Lagrangians. This is easy to see from the ChPT amplitude scaling, $M\sim p^{\nu}(p/m_{Z^{\prime}})^{2I_{Z^{\prime}}}$, where $\nu=2+2L+\sum_{i}V_{i}(d_{i}-2)$, with $L$ the number of loops, $d_{i}$ the number of derivatives in $V_{i}$ vertices of type $i$, and $I_{Z^{\prime}}$ the number of $Z^{\prime}$ internal lines. Each of the $p$ momenta in this scaling can be either of order ${\mathcal{O}}(m_{\pi,K,\eta})$ or ${\mathcal{O}}(m_{Z^{\prime}})$, where the latter we can in general take to be parametrically smaller. The LO vertices with $Z^{\prime}$ or $Z^{\prime 2}$ have $d_{i}=1$ and $d_{i}=0$. The one loop $Z^{\prime}$ contribution, $L=1,I_{Z^{\prime}}=1$, thus scales as $M\sim p^{4}/m_{Z^{\prime}}^{2}$, with $p\sim{\mathcal{O}}(m_{Z^{\prime}})$ or $p\sim{\mathcal{O}}(m_{K})$, showing that we need both ${\cal L}_{Q_{R,L}}^{(2)}$ and ${\cal L}_{Q_{R,L}}^{(4)}$ Lagrangians in order to capture all the counterterms.

We are now ready to show the one-loop results for $K-\bar{K}$ mixing. Expanding (87) and (88) in meson fields gives

[TABLE]

where we abbreviated $m_{uds}=m_{u}+m_{d}+m_{s}$, and $m_{ds}=m_{d}+m_{s}$. Expanding the NLO Lagrangian (87) in meson fields gives

[TABLE]

Above, we kept only the terms relevant for $K^{0}-\bar{K}^{0}$ mixing, defined $\phi_{1}\negthickspace\stackrel{{\scriptstyle\leftrightarrow}}{{\partial}}\negthickspace{}_{\mu}\phi_{2}=\phi_{1}\partial_{\mu}\phi_{2}-(\partial_{\mu}\phi_{1})\phi_{2}$, and replaced $\sqrt{2}f$ with the kaon decay constant, $f_{K}=155.6\pm 0.4$ MeV, to account for the SU(3) breaking. The first term in (90) is due to the axial vector coupling of the $Z^{\prime}$ boson and leads to the tree level contribution to the $K^{0}-\bar{K}^{0}$ mixing amplitude, Eq. (59). 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 (90) contribute to the middle diagram in Fig. 19, while the $Z^{\prime}_{\mu}Z^{\prime}{}^{\mu}$ term in the second line of (90) gives the last diagram in Fig. 19.

The one-loop contributions to $K-\bar{K}$ mixing are

[TABLE]

with $c_{\pi^{0}}=-1/\sqrt{2}$, $c_{\eta}=\sqrt{3/2}$, and $A_{0}$ and $B_{0}$ the Passarino-Veltman loop functions, see, e.g., Patel:2015tea . The counter-terms give

[TABLE]

Since there are no one-loop contributions of the form $\sim c_{d_{A}}^{12}+c_{d_{V}}^{12}$, this implies a relation between the counter-tems in the last line, m_{d}^{2}\big{(}C_{0}^{(4)}+C_{2}^{(4)}\big{)}=m_{s}^{2}C_{3}^{(4)}.

Appendix B Benchmarks

In this appendix we give the details of the two benchmarks, one for the case of the decoupled chains with $U(1)_{\rm FN}^{3}$ horizontal symmetry, and one for the coupled chains with the $U(1)_{\rm FN}$ horizontal symmetry.

B.1 The benchmark for the decoupled FN chains

The vector-like masses for our $U(1)_{\rm FN}^{3}$ benchmark are (in units of $10^{7}$ GeV)

[TABLE]

where the entries are for different nodes, for instance for the first generation FN fermions that are quark doublets, $M^{q_{1}}=\{M^{q_{1}}_{1},M^{q_{1}}_{2},M^{q_{1}}_{3}\}$, etc. For Yukawa couplings between FN fermions and the flavon, we have (in units of $10^{7}$ GeV)

[TABLE]

where again the entries correspond to values on different nodes. Note that the ratios $|(Y^{f}\langle\phi\rangle)_{ij}/(M^{f})_{ij}|$ are on average equal to $q=5$, but have a distribution that allows for ${\mathcal{O}}(1)$ deviations (in relative terms) from this value.

The Yukawa couplings on the zero node, which couple chiral fermions to the Higgs, are

[TABLE]

[TABLE]

[TABLE]

The absolute values, $r_{a}$, for each of the above entries were taken to be in the interval $r_{a}\in[0.3,0.9]$, when constructing the benchmark. Finally, the matrix of coefficients in the neutrino mass term, Eq. (25), is

[TABLE]

In the benchmark we set the kinetic mixing between different $Z_{i}^{\prime}$ to zero, as we do the kinetic mixing of $Z_{i}^{\prime}$ with hypercharge. The mass eigenstates, $Z_{i}^{\prime}$, 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. (50) are, for $Z_{1}^{\prime}$,

[TABLE]

while the nonzero imaginary entries for the couplings to quarks are $\operatorname{Im}[(c_{u_{R,1}})_{\{12,13\}}]=\{0.006,-0.001\}$, $\operatorname{Im}[(c_{d_{R,1}})_{\{12,13,23\}}]=\{0.508,0.233,0.069\}$, and for the couplings to leptons, $\operatorname{Im}[(c_{\ell_{R,1}})_{\{12,13\}}]=\{0.036,-0.002\}$, $\operatorname{Im}[(c_{\nu_{L},1})_{\{12,13,23\}}]=\{-0.624,0.521,0.196\}$. The coefficients below the diagonal, $i>j$, are given by $\operatorname{Im}(c_{a})_{ij}=-\operatorname{Im}(c_{a})_{ji}$ since $(c_{a})_{ij}=(c_{a})_{ji}^{*}$.

For $Z_{2}^{\prime}$ the real parts of the couplings in Eq. (50) are

[TABLE]

while the nonzero imaginary entries are $\operatorname{Im}[(c_{u_{L,2}})_{\{13,23\}}]=\{-0.002,0.009\}$,

$\operatorname{Im}[(c_{u_{R,2}})_{\{12,13,23\}}]=\{-0.003,-0.001,-0.036\}$, $\operatorname{Im}[(c_{d_{L,2}})_{\{13,23\}}]=\{-0.002,0.03\}$,

$\operatorname{Im}[(c_{d_{R,2}})_{\{12,13,23\}}]=\{-0.148,-0.381,0.617\}$, $\operatorname{Im}[(c_{\ell_{L,2}})_{\{12,13,23\}}]=\{0.317,-0.533,-1.014\}$, $\operatorname{Im}[(c_{\ell_{R,2}})_{\{12,23\}}]=\{-0.009,0.003\}$, $\operatorname{Im}[(c_{\nu_{L,2}})_{\{12,13,23\}}]=\{0.459,-0.842,0.793\}$, with $\operatorname{Im}(c_{a})_{ij}=-\operatorname{Im}(c_{a})_{ji}$.

For $Z_{3}^{\prime}$ the only nonzero coupling matrices are

[TABLE]

B.2 The benchmark for the coupled FN chains

The values of the vector-like mass matrices in complex plane for the $U(1)_{\rm FN}$ benchmark are shown in Fig. 20 (left) for $M^{q}_{n}$ in blue, $M_{n}^{u}$ in red, and $M_{n}^{d}$ in black, and in Fig. 20 (right) for $Y^{q}_{n}\langle\phi\rangle$ in blue, $Y_{n}^{u}\langle\phi\rangle$ in red, and $Y_{n}^{d}\langle\phi\rangle$ in black, in both cases in units of $10^{7}$ GeV. The values for the matrices on nodes $n=1(2,3)$ are denoted with a dot (star, square), with the labels denoting which element is being plotted, $(M_{n}^{f})_{ij}\to(i,j)$. The Yukawa couplings on the zero node and the matrix of coefficients in the neutrino mass term are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Fig. 20 shows that the entries are relatively uniformly distributed over the complex plane. In constructing the benchmark we restricted the values of the amplitudes to lie within the same interval, $r\in[0.3,0.9]$, that was used in the numerical scan in Section 3.2. The numerical values for the inputs are also available on request in the form of a Mathematica notebook.

After mass diagonalization, the real parts of the $Z^{\prime}$ couplings in Eq. (44) are

[TABLE]

while the imaginary parts of the couplings to quarks are $\operatorname{Im}[(c_{u_{L}})_{\{12,13,23\}}]=\{1.4,0.7,4.1\}\cdot 10^{-2}$, $\operatorname{Im}[(c_{u_{R}})_{\{12,13,23\}}]=\{-3.6,1,20.3\}\cdot 10^{-2}$, $\operatorname{Im}[(c_{d_{L}})_{\{12,13,23\}}]=\{-1.1,-1.2,-7.1\}\cdot 10^{-2}$, $\operatorname{Im}[(c_{d_{R}})_{\{12,13,23\}}]=\{1.2,3.1,-6.3\}\cdot 10^{-2}$, while the imaginary parts of the couplings to leptons are $\operatorname{Im}[(c_{\ell_{L}})_{\{12,13,23\}}]=\{-0.2,-0.4,1.1\}\cdot 10^{-2}$, $\operatorname{Im}[(c_{\ell_{R}})_{\{12,13,23\}}]=-\{0.3,-0.3,1.5\}\cdot 10^{-2}$, $\operatorname{Im}[(c_{\nu_{L}})_{\{12,13,23\}}]=\{-0.8,0.7,0.8\}\cdot 10^{-2}$. The elements below the diagonal are obtained from hermiticity of the coefficients, $(c_{a})_{ij}=(c_{a})_{ji}^{*}$, so that $\operatorname{Im}(c_{a})_{ij}=-\operatorname{Im}(c_{a})_{ji}$.