††institutetext: Department of Physics, University of Basel,
Klingelbergstrasse 82, CH-4056 Basel, Switzerland

# Froggatt-Nielsen ALP

Admir Greljo  Aleks Smolkovič  Alessandro Valenti<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The Froggatt-Nielsen (FN) mechanism, a prominent framework for explaining the
observed flavor hierarchies, generically predicts the existence of an axion-
like particle (ALP). This work examines a class of FN models based on
$\mathbb{Z}_{N}$ discrete symmetries. We chart the allowed parameter space
from a set of theoretical considerations and construct explicit renormalizable
completions with minimal field content necessary to generate consistent
textures. We then conduct comprehensive phenomenological analyses of two
particularly elegant $\mathbb{Z}_{4}$ and $\mathbb{Z}_{8}$ models,
highlighting the interplay between the effects of the ALP and the associated
UV fields. We find that the FN scale can be as low as a few TeV.

##  1 Introduction

The observed pattern of fermion masses and mixings cannot be explained within
the Standard Model (SM). Despite being generated from Yukawa interactions with
a single Higgs field, one particularly intriguing aspect is the presence of
roughly two-order-of-magnitude mass hierarchies between consecutive
generations of all three types of charged fermions: up quarks, down quarks,
and charged leptons. Additionally, the CKM mixing matrix Cabibbo:1963yz ;
Kobayashi:1973fv  is approximately a unit matrix, with hierarchies evident in
its off-diagonal elements. In stark contrast, the neutrino sector displays a
different behavior. Apart from the small overall mass scale, neutrino
oscillations reveal large mixing angles and a seemingly anarchic flavor
structure. The origin of generation hierarchies in the charged fermion sector,
in contrast to the neutrino sector, necessitates an organizing principle
beyond the SM.

The Froggatt-Nielsen (FN) mechanism Froggatt:1978nt ; Leurer:1992wg ;
Leurer:1993gy  is a prototypical example of how approximate flavor symmetries
can be used to generate Yukawa textures. This framework introduces a global
$\mathrm{U}(1)_{\text{FN}}$ symmetry, with different generations assigned
different integer charges and a symmetry-breaking spurion $\epsilon\ll 1$
carrying a unit charge. By judiciously selecting these charges, one can
construct Yukawa matrices $Y^{f}_{ij}$ that exhibit a desired suppression
factor $\epsilon^{n_{ij}}$. The simplicity of a FN model often conflicts with
its ability to accurately fit the observed flavor parameters Fedele:2020fvh ;
Cornella:2023zme . Opting for simpler charge assignments, with a relatively
small largest integer charge, results in increased variance in the
$\mathcal{O}(1)$ ultraviolet (UV) parameters. Nonetheless, such scenarios
might be preferable from a simplicity perspective, given the uncertainties in
the distribution of marginal couplings in the UV.

A straightforward way to realize $\epsilon$ is through the spontaneous
symmetry breaking of a perturbatively exact global $\mathrm{U}(1)_{\text{FN}}$
symmetry. This involves a complex scalar field $\Phi$ (a SM gauge singlet)
with charge 1, which acquires a vacuum expectation value (VEV)
$\langle\Phi\rangle=v_{\Phi}/\sqrt{2}$. Additionally, this framework assumes a
gap between $v_{\Phi}$ and the common UV scale $M$, which underlies the
generation of the FN effective field theory (EFT) operators. Consequently,
$\epsilon=v_{\Phi}/\sqrt{2}M$.

As a result of a spontaneously broken continuous global symmetry, a Goldstone
boson emerges in the spectrum, dominating infrared (IR) physics. This boson is
not exactly massless due to the chiral anomaly with QCD, where non-
perturbative instanton effects endow it with a mass, making it the QCD axion
Peccei:1977hh ; Wilczek:1977pj ; Weinberg:1977ma . The QCD axion, proposed to
address the strong CP problem, has recently been the focus of extensive
theoretical DiLuzio:2020wdo  and experimental studies Irastorza:2018dyq . The
concept of the axiflavon Calibbi:2016hwq  or flaxion Ema:2016ops  is
particularly compelling, as it simultaneously addresses two structural issues
of the Standard Model: the flavor puzzle and the strong CP problem. This leads
to intriguing phenomenology, characterized by the interplay between
astrophysical and flavor physics constraints. Notably, the flavor constraints
from $K^{+}\to\pi^{+}a$ impose the most stringent limits on the FN breaking
scale, requiring $v_{\Phi}\gtrsim\mathcal{O}(10^{12})$ GeV, as detailed in
MartinCamalich:2020dfe .

The $\mathrm{U}(1)$ symmetry can be explicitly broken in various ways, which
would undermine its solution to the strong CP problem. The most problematic
and perhaps unavoidable breaking arises from Planck-scale suppressed
operators, leading to the axion quality problem Georgi:1981pu ; Dine:1986bg ;
Barr:1992qq ; Kamionkowski:1992mf ; Holman:1992us . The primary objective of
this work is to explore a broader parameter space for the FN axion-like
particle (ALP) beyond the axiflavon limit while delegating the resolution of
the strong CP problem to another source. Explicit breaking can arise from
operators of different dimensions. For instance, consider a soft breaking by a
$m^{2}\phi^{2}$ operator, which predicts the ALP mass $m_{a}\sim m$. A
consistency condition is required to avoid disrupting the FN EFT power
counting and maintaining the desired textures: $m_{a}\lesssim v_{\Phi}$.

A systematic approach to exploring this terrain involves considering
perturbatively exact discrete subgroups. In this work, we focus on
$\mathbb{Z}_{N}$. We find that the simplest consistent model in non-
supersymmetric theories has $N=4$ and already fits the observed flavor
parameters at the level of recently proposed $\mathrm{SU}(2)$ gauged flavor
models Greljo:2023bix ; Greljo:2024zrj  or flavor deconstruction examples
Davighi:2023iks ; Greljo:2024ovt . Increasing $N$ allows for a finer fit to
the observed fermion masses and mixings. For instance, we demonstrate that
$\mathbb{Z}_{8}$ can realize the same Yukawa structure as the recently
proposed $\mathrm{U}(2)_{q+e^{c}+u^{c}}\times\mathbb{Z}_{2}$ symmetry
Antusch:2023shi .

Imposing a discrete symmetry instead of a continuous one eliminates the
presence of a Goldstone boson in the theory. Specifically, for
$\mathbb{Z}_{4}$, no light states are expected, with $m_{a}\sim v_{\Phi}$
predicted from the $\Phi^{4}$ term in the scalar potential. A light ALP is
anticipated for $N>4$, with its mass generated by a higher-dimensional
operator. In fact, the same scale $M$ responsible for the Yukawa textures can
also generate the ALP mass, as illustrated later by the wheel models (see Fig.
2). A natural outcome of this class of models is an FN ALP with mass $\gtrsim$
GeV, where the decay $K^{+}\to\pi^{+}a$ is kinematically forbidden, as
illustrated by the $\mathbb{Z}_{8}$ model in Section 4.

This raises an intriguing question: what is the lowest FN scale that remains
compatible with experimental constraints? This is a quantitative question with
an interesting interplay between ALP-induced effects and UV contributions. In
this paper, we conduct a thorough phenomenological study performing precision
calculations in explicit renormalizable models, completing the FN
$\mathbb{Z}_{4}$ and $\mathbb{Z}_{8}$ structures with minimalistic
ingredients. The philosophy here is that rather than estimating UV
contributions through FN EFT power counting, we introduce only the necessary
UV field content to realize the FN mechanism and study its irreducible
effects. As we demonstrate for the chosen models, the limit on $v_{\Phi}$ is
weaker than the one estimated using the FN EFT. In contrast to the axiflavon
case, this scenario motivates direct searches for heavy new physics at future
colliders.

The paper is organized as follows. In Section 2, we provide a theoretical
discussion of $\mathbb{Z}_{N}$ models, charting the parameter space of the FN
ALP. Section 3 focuses on constructing a concrete model based on
$\mathbb{Z}_{4}$ symmetry, fitting the SM flavor parameters, and presenting a
comprehensive survey of experimental bounds. In Section 4, we construct and
examine an explicit $\mathbb{Z}_{8}$ model. Finally, in Section 5, we conclude
our study. Technical details are discussed in the Appendices.

##  2 $\mathbb{Z}_{N}$ Froggatt-Nielsen

This section introduces a class of discrete FN models based on
$\mathbb{Z}_{N}$ symmetries and provides a theoretical discussion of the ALP
parameter space.

#### Yukawa textures

Consider an EFT comprised of the SM matter content supplemented by a singlet
complex scalar field $\Phi$. Both $\Phi$ and the SM fermions are assumed to be
charged under a global, perturbatively exact111For realistic choices of the
chiral fermion charges, the $\mathbb{Z}_{N}$ symmetry is typically anomalous.
$\mathbb{Z}_{N}$ symmetry, with $[\Phi]=1$. The mechanism generating
hierarchies in the SM Yukawas is of a Froggatt-Nielsen type:

$\displaystyle-\mathcal{L}\supset\sum_{f,F}\left[x_{ij}^{f}\left(\frac{\Phi}{M}\right)^{n_{ij}^{f}}+x_{ij}^{f\prime}\left(\frac{\Phi^{*}}{M}\right)^{n_{ij}^{f\prime}}\right]\bar{F}_{L,i}Hf_{R,j}+\\{\widetilde{H}\\}+\text{h.c.},$
(1)

where $x_{ij}^{f},x_{ij}^{f\prime}$ are a priori $\mathcal{O}(1)$ couplings,
and $n_{ij}^{f}$ and $n_{ij}^{f\prime}$ are determined by the different
$\mathbb{Z}_{N}$ charges of the SM fields. Here, we assume a common mass scale
$M$ that generates the effective operators. The term $\\{\widetilde{H}\\}$ is
a shorthand notation for the interactions involving
$\widetilde{H}\equiv\varepsilon H^{*}$, while $f$ and $F$ denote SM
$SU(2)_{\text{L}}$ singlets and doublets, respectively. $H$ is uncharged under
$\mathbb{Z}_{N}$, $[H]=0$.

It is instructive to think of $\mathbb{Z}_{N}\subset\mathrm{U}(1)$. The
$\mathrm{U}(1)$-breaking but $\mathbb{Z}_{N}$-preserving interactions
involving $\Phi^{*}$ ($\Phi$) represent an important difference compared to
the “standard” $\mathrm{U}(1)$ case. The $\mathbb{Z}_{N}$, but not the
$\mathrm{U}(1)$, allows the replacement
$\Phi^{n_{ij}^{f}}\rightarrow(\Phi^{*})^{N-n_{ij}^{f}}$ in Eq. (1), since
$[\Phi]^{*}=-[\Phi]\equiv(N-[\Phi])\,\>\text{mod}\>N$. To avoid spoiling the
hierarchies generated exclusively by $\braket{\Phi}$, it is required that
$n_{ij}^{f}\leq\left\lfloor\frac{N}{2}\right\rfloor$. Therefore, special care
must be taken with charge assignments to ensure the desired textures. A
sufficient condition is that chiral fermions of the SM are assigned charges
from the set $\\{0,\dots,\left\lfloor\frac{N}{2}\right\rfloor\\}$. For
example, the smallest group that generates hierarchies between consecutive
families by setting each family distinct charges $\\{2,1,0\\}$ is
$\mathbb{Z}_{4}$ instead of $\mathbb{Z}_{3}$.222In this work, we focus on non-
supersymmetric theories. Holomorphy, however, would make the $\mathbb{Z}_{3}$
symmetry eligible. For a study of $\mathbb{Z}_{N}$ FN models in supersymmetric
context, see Higaki:2019ojq . Therefore, we will focus only on models with
$N\geq 4$.

The numerical value for the spurion $\epsilon=\langle\Phi\rangle/\sqrt{2}M$
depends non-trivially on $N$ and the specific details of the model. However, a
simple estimate can be made by requiring the largest hierarchy to be saturated
by $\epsilon^{\left\lfloor\frac{N}{2}\right\rfloor}\sim y_{e}/y_{t}$, where
the ratio of the electron to top quark Yukawa couplings $y_{e}/y_{t}\sim
10^{-6}$. Since $\epsilon\xrightarrow[]{N\rightarrow\infty}1$, large $N$
setups are disfavored not only from the simplicity perspective but also due to
concerns regarding the consistency of the FN EFT expansion.

Figure 1: Theoretical constraints on the FN ALP, presented in the
$(m_{a},1/v_{\Phi})$ plane, arising from Eqs. (6), (10), and (11), with
$\lambda=0.5$. The plots have been drawn for $N=4$ and $8$; in particular, the
thickness of the darker gray and red lines represents the variability of the
boundary from $N=8$ (less constraining) to $N=4$ (more constraining). The
expectation for a representative $N=4$ model with
$\lambda^{\prime}_{4}=\lambda/4$ is depicted with a dotted line, while the
band enclosed by the dashed lines corresponds to the expectation around Eq.
(9) for a $N=8$ “wheel” model. For details, see Section 2.

#### Charting the ALP

Theoretical constraints that define the ALP’s parameter space of interest are
illustrated in Fig. 1. This figure presents various constraints, which will be
discussed below, in the plane of ALP mass versus the inverse FN scale. The
most general potential compatible with the $\mathbb{Z}_{N}$ symmetry reads

$\displaystyle\begin{aligned}
V(\Phi)&=-m^{2}|\Phi|^{2}+\frac{1}{2}\lambda|\Phi|^{4}-\frac{1}{4}\frac{\lambda_{N}^{\prime}}{M_{\Phi}^{N-4}}\left[\Phi^{N}+(\Phi^{*})^{N}\right]+\dots,\\\
\end{aligned}$ (2)

where without loss of generality, we take
$\lambda^{\prime}_{N}\in\mathds{R}_{>0}$. For simplicity, we are neglecting
any portals with the SM Higgs doublet and further higher-dimensional
operators, which will lead to subleading corrections in the following
discussion.

Let us employ a convenient parametrization
$\Phi=\left(\frac{v_{\Phi}+\rho}{\sqrt{2}}\right)e^{ia/v_{\Phi}}$, where
$\rho$ and $a$ represent the radial and angular (ALP) modes, respectively,
with $\langle\Phi\rangle=v_{\Phi}/\sqrt{2}$. The SM Yukawa hierarchies are
explained by a small parameter $\epsilon\equiv v_{\Phi}/(\sqrt{2}M)$ after
expanding Eq. (1). Since, in general, $M$ and $M_{\Phi}$ may differ, we define
an effective coupling

$\displaystyle\lambda_{4,\text{eff}}^{\prime}=\lambda_{N}^{\prime}\left(\frac{M}{M_{\Phi}}\right)^{N-4}\epsilon^{N-4}\,.$
(3)

The minima of the potential in Eq. (2) are set at
$v_{\Phi}^{2}=2m^{2}/(\lambda-N\lambda_{4,\text{eff}}^{\prime}/4)$ and
$\braket{a/v_{\Phi}}=2\pi k/N$, $k=1,2,\dots N$. The masses of the radial and
angular modes read

$\displaystyle
m_{\rho}^{2}=\left(\lambda-\frac{1}{8}N(N-2)\lambda_{4,\text{eff}}^{\prime}\right)v_{\Phi}^{2}\quad\text{and}\quad
m_{a}^{2}=\frac{1}{8}N^{2}\lambda_{4,\text{eff}}^{\prime}v_{\Phi}^{2}\,.$ (4)

They satisfy a sum rule:

$\displaystyle m_{\rho}^{2}+\left(\frac{N-2}{N}\right)m_{a}^{2}=\lambda
v_{\Phi}^{2}\,.$ (5)

This, combined with the stability conditions $m_{\rho}^{2}\geq 0$ and
$m_{a}^{2}\geq 0$, provides an admissible range for the $\rho$ and $a$
masses:333This also holds when taking into account possible subleading
corrections to Eq. (2) by defining effective couplings
$\lambda_{\text{eff}}=\lambda(1+\mathcal{O}(\epsilon^{2}))$ and
$\lambda_{\text{4,
eff}}^{\prime}=\lambda_{N}^{\prime}\left(M/M_{\Phi}\right)^{N-4}\epsilon^{N-4}(1+\mathcal{O}(\epsilon^{2}))$.
Higher-order harmonics cannot, however, be described in this simple picture.
They are suppressed by a factor of $\epsilon^{N}$ compared to the leading term
and can thus be safely neglected.

$\displaystyle\phantom{\qquad\textbf{[Stability condition]}}\begin{aligned}
0\leq m_{\rho}^{2}&\leq\lambda v_{\Phi}^{2},\\\ 0\leq
m_{a}^{2}&\leq\left(\frac{N}{N-2}\right)\lambda
v_{\Phi}^{2}.\end{aligned}\qquad\textbf{[Stability condition]}$ (6)

The excluded region is shown in gray in Fig. 1. With marginal coupling
$\lambda$ of $\mathcal{O}(1)$, we generally expect $m_{\rho}\sim v_{\Phi}$.
The $N=4$ case is unique because the $\mathrm{U}(1)$-breaking operator is
renormalizable, resulting in $m_{a}\sim v_{\Phi}$ for
$\lambda^{\prime}_{N}\sim\mathcal{O}(1)$ and, consequently, no ALP.444The
stability of the $\mathbb{Z}_{4}$ potential requires
$\lambda_{4}^{\prime}<\lambda$. In the phenomenology section, we will also
consider $\lambda_{N}^{\prime}\ll 1$, in which case $a$ is an ALP. However,
for $N>4$, there is an ALP with mass $m_{a}\ll v_{\Phi}$.

Let us sketch a renormalizable UV model example that generates
$\lambda_{N}^{\prime}$ for $N>4$ and predicts the ALP mass. Any completion
with only a minimal set of vector-like fermions (VLFs) required to generate
Yukawa operators in Eq. (1) with desired textures typically leads to an
accidental $\mathrm{U}(1)_{\text{FN}}$ symmetry when these fields are
integrated out, thus failing to generate $\lambda_{N}^{\prime}$. An elegant
mechanism to generate $\lambda_{N}^{\prime}$ would be to consider a _full_ set
of VLFs $\Psi_{k}$ with $\mathbb{Z}_{N}$ charges $[\Psi_{k}]=k$ and couplings
$\sum_{k}y_{k}\,\Phi\,\bar{\Psi}_{L,k+1}\Psi_{R,k}$. The “wheel” diagram in
Fig. 2 predicts

$\displaystyle\frac{\lambda_{N}^{\prime}}{M_{\Phi}^{N-4}}\sim\frac{1}{M_{\Psi}^{N-4}}\frac{\prod_{k}y_{k}}{16\pi^{2}}\,,$
(7)

under the assumption of approximately degenerate $\Psi_{i}$ masses.555A
straightforward mechanism to achieve this would be to introduce an additional
real scalar field $\chi$ such that $(\chi,\Phi,\Psi_{R,k})$ are odd under a
new $\mathbb{Z}_{2}$ symmetry. When $\chi$ acquires a non-zero vacuum
expectation value, the interactions $\chi\bar{\Psi}_{L,k}\Psi_{R,k}$ result in
comparable masses among all VLFs.

Incidentally, a subset of these VLFs, with the appropriate couplings, is
sufficient to generate the Yukawa operators in Eq. (1). Therefore, in this
class of models,

$\displaystyle M\sim M_{\Phi}\sim M_{\Psi}.$ (8)

This, finally, allows us to predict the ALP mass using Eq. (4),

$\frac{m_{a}}{v_{\Phi}}\sim\frac{1}{8\pi}N\epsilon^{\frac{N}{2}-2}\sqrt{\prod_{k}y_{k}}~{},$
(9)

where we take an estimate for
$\epsilon\sim(y_{e}/y_{t})^{\left\lfloor\frac{N}{2}\right\rfloor^{-1}}$, as
discussed earlier. The prediction for the $N=8$ case is shown in Fig. 1 as a
region within the two black dashed lines which represent the variation
$0.001<\sqrt{\prod_{k}y_{k}}<0.2$.666Choosing $y_{k}$ from a flat distribution
in the $(0,1)$ range, there is $90\%$ probability for $\sqrt{\prod_{k}y_{k}}$
to be above (below) 0.001 (0.2).

$\Phi$$\Phi$$\Phi$$\Phi$$\Phi$$\Phi$$\Phi$$\Psi_{1}$$\Psi_{2}$$\Psi_{3}$$\Psi_{N-2}$$\Psi_{N-1}$$\Psi_{N}$

Figure 2: The “wheel” diagram generating the $\mathrm{U}(1)_{\text{FN}}$
breaking term $\Phi^{N}+(\Phi^{*})^{N}$. For details see Section 2.

This mechanism might, however, not be present at all, and the
$\mathrm{U}(1)_{\text{FN}}$-breaking term in Eq. (2) could be generated by UV
physics at scales $M_{\Phi}$ significantly larger than $M$. In the most
extreme scenario, the $\mathrm{U}(1)_{\text{FN}}$ symmetry is broken solely by
quantum gravity effects $M_{\Phi}\sim M_{\text{Pl}}$, which are expected to
break any global symmetry Palti:2019pca .777In this scenario, we assume that
$\mathbb{Z}_{N}$ symmetry is protected, possibly through an appropriate
gauging mechanism Krauss:1988zc ; Banks:1989ag ; Preskill:1990bm ;
Ibanez:1991hv . If this were the case, the coupling
$\lambda_{4,\text{eff}}^{\prime}$ would be lower bounded as
$\lambda_{4,\text{eff}}^{\prime}\gtrsim
16\pi^{2}(M/M_{\text{Pl}})^{N-4}\epsilon^{N-4}$, barring any fine-tuning on
Planck-suppressed operators.888Here $M_{\text{Pl}}$ denotes the reduced Planck
mass, $M_{\text{Pl}}\simeq 2.4\times 10^{18}$ GeV. This leads, in turn, to a
lower bound on the ALP mass:

$\displaystyle\phantom{\qquad\textbf{[QG condition]}}m_{a}^{2}\gtrsim
8\pi^{2}\left(\frac{N^{2}}{2^{N/2}}\right)\left(\frac{v_{\Phi}}{M_{\text{Pl}}}\right)^{N-4}v_{\Phi}^{2}\,.\qquad\textbf{[QG
condition]}$ (10)

This lower bound is analogous in spirit to the constraint on the dimension of
PQ-breaking operators encountered in the well-known QCD axion quality problem
Georgi:1981pu ; Dine:1986bg ; Barr:1992qq ; Kamionkowski:1992mf ;
Holman:1992us . The “excluded” region for $N=8$ is shown in blue in Fig. 1.

Finally, one should not forget QCD itself and its contribution to the axion
mass, which is approximately given by $(m_{a}^{2})_{\text{QCD}}\sim
m_{\pi}^{2}f_{\pi}^{2}c_{a}^{2}/v_{\Phi}^{2}$, where $c_{a}$ is the anomaly
coefficient entering $(c_{a}g^{2}_{s}/32\pi^{2})G\widetilde{G}$, and we
approximate $c_{a}\sim N$. As the physical mass of the ALP is given by the sum
of this contribution and the one of Eq. (4), any value of $m_{a}$ lower than
that of the single contributions requires a non-trivial cancellation. Avoiding
this tuning can then be stated as

$\displaystyle\phantom{\qquad\textbf{[QCD
condition]}}m_{a}\gtrsim(m_{a})_{\text{QCD}}.\qquad\textbf{[QCD condition]}$
(11)

The red region shown in Fig. 1 does not satisfy this condition.999In the
presence of multiple ALPs, this condition applies only to the linear
combination entering the QCD potential. Hence, there might be additional
lighter fields Gavela:2023tzu .

To conclude, Eqs. (6), (10), and (11) together define an acceptable region for
$m_{a}$ and $v_{\Phi}$ that a $\mathbb{Z}_{N}$ FN model should respect. This
corresponds to the white region in Fig. 1, where we also show model
expectations for $N=\\{4,8\\}$ as discussed above.

[] Fields | $\mathrm{SU}(3)_{C}$ | $\mathrm{SU}(2)_{L}$ | $\mathrm{U}(1)_{Y}$ | $\mathbb{Z}_{4}$
---|---|---|---|---
$q_{1,2,3}$ | 3 | 2 | $\frac{1}{6}$ | $\\{2,1,0\\}$
$\ell_{i}$ | 1 | 2 | $\vbox{\hbox{\scalebox{0.6}[1.0]{$-$}}}\frac{1}{2}$ | $0$
$e_{1,2,3}$ | 1 | 1 | $\vbox{\hbox{\scalebox{0.6}[1.0]{$-$}}}1$ | $\\{\vbox{\hbox{\scalebox{0.6}[1.0]{$-$}}}2,\vbox{\hbox{\scalebox{0.6}[1.0]{$-$}}}1,0\\}$
$u_{i}$ | 3 | 1 | $\frac{2}{3}$ | $0$
$d_{i}$ | 3 | 1 | $\vbox{\hbox{\scalebox{0.6}[1.0]{$-$}}}\frac{1}{3}$ | $0$
$Q^{a}_{2}$ | 3 | 2 | $\frac{1}{6}$ | $0$
$Q_{1}$ | 3 | 2 | $\frac{1}{6}$ | $1$
$E^{a}_{2}$ | 1 | 1 | $\vbox{\hbox{\scalebox{0.6}[1.0]{$-$}}}1$ | $0$
$E_{1}$ | 1 | 1 | $\vbox{\hbox{\scalebox{0.6}[1.0]{$-$}}}1$ | $1$
$H$ | 1 | 2 | $\frac{1}{2}$ | $0$
$\Phi$ | 1 | 1 | $0$ | $1$
Table 1: Field content of the UV model as described in Section 3. The rows in
this table are grouped into four categories. The first and second represent
left and right chiral fermions, respectively. The third consists of vector-
like fermions, and the last comprises scalars. Flavor indices are denoted as
$i=1,2,3$ and $a=1,2$. See Section 3.1 for details.

##  3 The minimal model: $\mathbb{Z}_{4}$

This section constructs and thoroughly analyses an explicit $\mathbb{Z}_{4}$
FN model. We begin in Section 3.1 by defining a concrete renormalizable model
that realizes the $\mathbb{Z}_{4}$ FN texture, and we examine the EFT
structure it generates. Following this, Section 3.2 offers a comprehensive
analysis of the various observables predicted by the model. Section 3.3
collects constraints in $(m_{a},v_{\Phi}^{-1})$ plane and discusses their
broader implications.

###  3.1 Setup

The $\mathbb{Z}_{4}$ symmetry offers the simplest realization capable of
generating the SM fermion hierarchies. It features a renormalizable potential
given by the $N\rightarrow 4$ limit of Eq. (2), ensuring no massless spin-0
states exist. Therefore, this case can be considered minimal in terms of UV
matter content, as it only requires completing the Yukawa operators in Eq.
(1).

##### UV model.

The field content of a particular model incarnating this idea is listed in
Table 1. Apart from $\Phi$ and the SM fields, the model includes only the
addition of vector-like quarks $Q_{1},Q_{2}^{a=1,2}$ and vector-like leptons
$E_{1},E_{2}^{a=1,2}$.101010Alternatively, an even more minimal completion
involves replacing VLFs with three additional Higgs doublets carrying
$\mathbb{Z}_{4}$ charges of $\\{1,2,3\\}$. Their masses are $\mathcal{O}(M)$,
and they do not develop a VEV. The non-trivial FN charges under the
$\mathbb{Z}_{4}$ symmetry are $[Q_{1}]=-[E_{1}]=1$. The only charged SM fields
are $q$ and $e$, with $[q_{1}]=-[e_{1}]=2$ and $[q_{2}]=-[e_{2}]=1$.111111This
charge assignment is inspired by the recently proposed $\mathrm{U}(2)_{q+e}$
flavor symmetry Antusch:2023shi . Charging under the flavor group only left-
handed quarks and right-handed leptons produces the observed mass hierarchies
of electrically charged fermions and the mixing hierarchies of left-handed
quarks while maintaining an anarchic flavor structure for neutrinos, ensuring
large neutrino mixing. Our $\mathbb{Z}_{4}$ is the smallest flavor symmetry
resulting in similar Yukawa textures as predicted by $\mathrm{U}(2)_{q+e}$ of
Antusch:2023shi . All the other fields are taken as $\mathbb{Z}_{4}$-neutral,
except $\Phi$. The Yukawa sector of the most general renormalizable Lagrangian
compatible with this symmetry reads121212Note that Eq. (12) respects an
accidental $\mathrm{U}(1)_{\text{FN}}$. If we would have chosen a
$\mathbb{Z}_{3}$ symmetry instead of $\mathbb{Z}_{4}$, additional couplings,
such as $\overline{q}_{1}\Phi^{*}Q^{a}_{2}$, would appear, violating
$\mathrm{U}(1)_{\text{FN}}$. This would disrupt the desired Yukawa textures,
as discussed below Eq. (1).

$\begin{split}\mathcal{L}_{\text{UV}}&\supset
x_{1}^{q}\Phi\overline{q}_{1}Q_{1}+x_{12}^{qa}\Phi\overline{Q}_{1}Q_{2}^{a}+x_{2}^{qa}\Phi\overline{q}_{2}Q_{2}^{a}-y_{j}^{da}\overline{Q}_{2}^{a}Hd_{j}-y_{j}^{ua}\overline{Q}_{2}^{a}\tilde{H}u_{j}\\\
&+x_{1}^{e}\Phi\overline{E}_{1}e_{1}+x_{12}^{ea}\Phi\overline{E}_{2}^{a}E_{1}+x_{2}^{ea}\Phi\overline{E}_{2}^{a}e_{2}-y_{i}^{ea}\overline{\ell}_{i}HE_{2}^{a}\\\
&-z_{j}^{d}\bar{q}_{3}Hd_{j}-z_{j}^{u}\bar{q}_{3}\tilde{H}u_{j}-z_{i}^{e}\bar{\ell}_{i}He_{3}+\text{h.c.}\,.\end{split}$
(12)

Flavor symmetries of the kinetic terms are used to perform field redefinitions
and eliminate mass mixings between chiral and vector-like fermions of the same
representation without any loss of generality. The specific basis and
parameterization used are detailed in A.1.2.

##### Tower of EFTs.

A systematic analysis of this model involves constructing a tower of effective
field theories, as illustrated in Fig. 3. In the first step, the VLF states
are integrated out, matching to the SM+$\Phi$ EFT. The resulting set of local
higher-dimensional operators can be categorized into two broad classes based
on their involvement of the $\Phi$ field. First are operators that include the
$\Phi$ field. A subset of these include Eq. (1), which will generate the SM
Yukawa couplings at the FN scale $v_{\Phi}$ upon matching to the next EFT, the
SM+$a$ EFT. This class of operators also generates $a$ couplings. As already
discussed, $m_{a}$ derived in Eq. (4) is governed by $\lambda_{N}^{\prime}$.
Although the natural expectation $\lambda_{N}^{\prime}\sim 1$ predicts
$m_{a}\sim v_{\Phi}$, in this section, we treat $\lambda_{N}^{\prime}$ as a
free parameter that can be arbitrarily small.131313It is technically natural
to take $\lambda_{4}^{\prime}\ll 1$, though this opposes the spirit of our
narrative. Our goal is to explore a broader parameter space of ALP masses,
deriving general lessons that relate to $N>4$ models where ALP is expected to
be light. Second, another class of operators exists that does not involve the
$\Phi$ field. These operators violate the accidental flavor symmetries of the
SM and induce significant UV contributions to rare and forbidden SM
transitions, as we will see in a moment.

Figure 3: Illustration of the hierarchical scales within the model and the
corresponding layers of effective field theory descriptions.

##### Generating the SM Yukawas.

Integrating out the VLF at the tree-level leads to the effective Yukawa
Lagrangian of Eq. (1) with $x_{ij}^{f\prime}=0$ and

$n_{ij}^{d}=n_{ij}^{u}=n_{ij}^{e}=\begin{pmatrix}2&2&2\\\ 1&1&1\\\
0&0&0\end{pmatrix}.$ (13)

Representative Feynman diagrams are shown in Fig. 4. The matching of
$x_{ij}^{f}$ to the UV parameters in Eq. (12) is provided in Appendix A,
together with an explicit benchmark point in the parameter space correctly
reproducing the SM fermion masses and mixings, see Eqs. (89) and (93).
Interestingly, their values can be correctly predicted with most of the UV
couplings in Eq. (12) within a factor of 10 from each other. This stems from
the fact that the light SM Yukawas are expressed as products of many UV
couplings, which can naturally contribute to the hierarchy if these are
somewhat small, say $\mathcal{O}(0.3)$.141414A similar effect is observed in
Fig. 3 of Greljo:2024zrj . The only exceptions are $y_{b}$ and $y_{\tau}$,
which are here not dynamically generated, and hence we require them to be of
$\mathcal{O}(10^{-2})$.

Figure 4: Leading diagrams in the $\mathbb{Z}_{4}$ model generating effective
down-quark Yukawa operators in Eq. (1), and analogously for up quarks and
charged leptons. For further details, refer to Section 3.1.

##### ALP couplings.

For the upcoming flavor phenomenology discussion, we express the couplings of
the radial and angular modes of $\Phi$ in the SM fermion mass basis as

$\mathcal{L}\supset-\sum_{f=u,d,e}c_{ij}^{f}\bar{f}_{Li}\left(\rho+ia\right)f_{Rj}+\text{h.c.}+\text{h.o.}\,,$
(14)

where the coupling matrices can be generally expressed as

$c^{f}=\frac{1}{v_{\Phi}}\left[\left(U_{L}^{f\dagger}Q^{F}U_{L}^{f}\right)\hat{m}^{f}-\hat{m}^{f}\left(U_{R}^{f\dagger}Q^{f}U_{R}^{f}\right)\right]\,.$
(15)

Here $U_{L}$ and $U_{R}$ are the left- and right-handed fermion rotation
matrices that diagonalize the corresponding mass terms, $\hat{m}$ are the
diagonal fermion mass matrices, and $Q$ are diagonal matrices encoding the SM
fermion FN charges. The derivation of this formula is presented in Appendix B.
In our concrete realization of the FN model, only the left-handed quarks are
charged, and the down-type quark coupling matrix can be written as

$c^{d}=\frac{1}{v_{\phi}}\left(V_{\mathrm{CKM}}^{\dagger}Q^{d}V_{\mathrm{CKM}}\right)\hat{m}^{d}\,,$
(16)

where we assumed that the CKM matrix is fully given by
$U_{L}^{d}=V_{\mathrm{CKM}}$. This is a sensible assumption as the bigger
hierarchies in the up sector lead to much more suppressed mixing angles in
$U_{L}^{u}$, as shown in App. A. Interestingly, the down-quark flavor-changing
neutral currents due to $a,\rho$ are then fully determined by the measured
values of the down-quark masses and the CKM matrix, up to the overall
normalization of $v_{\Phi}$. Using the Wolfenstein parameterization, we get

$c^{d}=\frac{1}{v_{\phi}}\left(\begin{array}[]{ccc}2m_{d}&\lambda
m_{s}&A\lambda^{3}m_{b}(-2i\eta+2\rho-1)\\\ \lambda
m_{d}&m_{s}&A\lambda^{2}m_{b}\\\
A\lambda^{3}m_{d}(2i\eta+2\rho-1)&A\lambda^{2}m_{s}&A^{2}\lambda^{4}m_{b}\\\
\end{array}\right)\,,$ (17)

where in each entry, we keep only the leading nontrivial term in $\lambda$.
The $c^{d}_{33}$ entry is small as a result of the fact that both the left and
right components of the third generation are uncharged under $\mathbb{Z}_{4}$.
Hence, the coupling only arises due to a small mixing. This is a peculiarity
of the $\mathbb{Z}_{4}$ symmetry, which does not hold for larger $N$ where one
attempts to address the smallness of $y_{b}$ and $y_{\tau}$ dynamically. For
the purposes of phenomenology, we also report the lowest order in which we get
complex off-diagonal entries:

$\displaystyle\begin{split}v_{\Phi}c^{d}_{12}&=\lambda
m_{s}-\lambda^{3}m_{s}/2+\lambda^{5}m_{s}(A^{2}(i\eta-\rho+1)-1/8)\,,\\\
v_{\Phi}c^{d}_{21}&=\lambda
m_{d}-\lambda^{3}m_{d}/2+\lambda^{5}m_{d}(A^{2}(-i\eta-\rho+1)-1/8)\,,\\\
v_{\Phi}c^{d}_{23}&=A\lambda^{2}m_{b}+A\lambda^{4}m_{b}(-2i\eta+2\rho-1/2)\,,\\\
v_{\Phi}c^{d}_{32}&=A\lambda^{2}m_{s}+A\lambda^{4}m_{s}(2i\eta+2\rho-1/2)\,.\end{split}$
(18)

Similarly, for $c^{u}$ and $c^{e}$ we find

$\displaystyle\begin{aligned}
c^{u}&=\frac{1}{v_{\Phi}}\begin{pmatrix}2m_{u}&m_{c}\frac{m_{u}}{m_{c}}\frac{z_{u_{2}}}{z_{u_{1}}}&m_{t}\frac{m_{u}}{m_{t}}\frac{2y_{u_{2}}z_{u_{3}}-y_{u_{3}}z_{u_{2}}}{y_{u_{2}}z_{u_{1}}}\\\
m_{u}\frac{m_{u}}{m_{c}}\frac{z_{u_{2}}}{z_{u_{1}}}&m_{c}&m_{t}\frac{m_{c}}{m_{t}}\frac{y_{u_{3}}}{y_{u_{2}}}\\\
m_{u}\frac{m_{u}}{m_{t}}\frac{2y_{u_{2}}z_{u_{3}}^{*}-y_{u_{3}}z_{u_{2}}}{y_{u_{2}}z_{u_{1}}}&m_{c}\frac{m_{c}}{m_{t}}\frac{y_{u_{3}}}{y_{u_{2}}}&m_{c}\frac{m_{c}}{m_{t}}\frac{y_{u_{3}}^{2}}{y_{u_{2}}^{2}}\end{pmatrix},\\\
c^{e}&=\frac{1}{v_{\Phi}}\begin{pmatrix}2m_{e}&m_{e}\frac{m_{e}}{m_{\mu}}\frac{z_{e_{2}}}{z_{e_{1}}}&m_{e}\frac{m_{e}}{m_{\tau}}\frac{2y_{e_{2}}z_{e_{3}}^{*}-y_{e_{3}}z_{e_{2}}}{y_{e_{2}}z_{e_{1}}}\\\
m_{\mu}\frac{m_{e}}{m_{\mu}}\frac{z_{e_{2}}}{z_{e_{1}}}&m_{\mu}&m_{\mu}\frac{m_{\mu}}{m_{\tau}}\frac{y_{e_{3}}}{y_{e_{2}}}\\\
m_{\tau}\frac{m_{e}}{m_{\tau}}\frac{2y_{e_{2}}z_{e_{3}}-y_{e_{3}}z_{e_{2}}}{y_{e_{2}}z_{e_{1}}}&m_{\tau}\frac{m_{\mu}}{m_{\tau}}\frac{y_{e_{3}}}{y_{e_{2}}}&m_{\mu}\frac{m_{\mu}}{m_{\tau}}\frac{y_{e_{3}}^{2}}{y_{e_{2}}^{2}}\end{pmatrix},\end{aligned}$
(19)

where $y_{u_{i}},z_{u_{i}},y_{e_{i}},z_{e_{i}}$ are parameters that are not
determined within the EFT but can be correlated to the UV parameters in Eq.
(12) as shown in App. A. The matrices have been written factorizing the
mixing, which renders evident the fact that ALP-induced flavor violation is
suppressed by the ratio of the fermion masses. Again, the couplings to the
third generation are small due to the fact that these are not charged under
$\mathbb{Z}_{4}$. We stress that $c^{d}$ of Eq. (17) features the same pattern
upon recognizing that $\lambda\sim m_{d}/m_{s},A\lambda^{2}\sim m_{s}/m_{b}$
and $A\lambda^{3}\sim m_{d}/m_{b}$.

#####  $\mathbb{Z}_{4}$-symmetric UV operators.

Integrating out VLFs at the UV matching scale directly generates certain SMEFT
operators that violate approximate flavor symmetries and are crucial for the
phenomenological discussion. We refer to these as “VLF contributions”, and we
will carefully study their interplay with the ALP contributions. Already at
the tree level, the model generates

$(H^{\dagger}i\overset{\text{\footnotesize$\leftrightarrow$}}{D}_{\mu}H)(\bar{\ell}_{i}\gamma^{\mu}\ell_{j}),\quad(H^{\dagger}i\overset{\text{\footnotesize$\leftrightarrow$}}{D}_{\mu}H)(\bar{d}_{i}\gamma^{\mu}d_{j}),\quad(H^{\dagger}i\overset{\text{\footnotesize$\leftrightarrow$}}{D}_{\mu}H)(\bar{u}_{i}\gamma^{\mu}u_{j}),$
(20)

with arbitrary flavor structure and Wilson coefficients $\sim M^{-2}$. These
operators provide a leading contribution to $\Delta F=1$ transitions such as
$\epsilon^{\prime}/\epsilon$, $\mu\to e$ conversion, $K_{L}\to\mu^{+}\mu^{-}$,
etc. Additionally, under renormalization, they mix into $\Delta F=2$
operators. The primary $\Delta F=2$ contributions, however, arise from four-
fermion operators generated through one-loop matching at the UV scale. These
are

$(\bar{d}_{i}\gamma^{\mu}d_{j})(\bar{d}_{k}\gamma_{\mu}d_{l}),\quad(\bar{u}_{i}\gamma^{\mu}u_{j})(\bar{u}_{k}\gamma_{\mu}u_{l}),$
(21)

with Wilson coefficients $\sim(4\pi M)^{-2}$. The flavor violation in Eq. (21)
is expected to be of $\mathcal{O}(1)$ since our model does not impose any
selection rules. This will result in stringent constraints from neutral meson
mixing, most notably $\epsilon_{K}$. Finally, lepton flavor violating four-
fermion operators

$(\bar{\ell}_{i}\gamma^{\mu}\ell_{j})(\bar{d}_{k}\gamma_{\mu}d_{l}),\quad(\bar{\ell}_{i}\gamma^{\mu}\ell_{j})(\bar{u}_{k}\gamma_{\mu}u_{l}),$
(22)

are generated at the one-loop level and they represent the leading
contribution to LFV meson decays, such as $K_{L}\to e\mu$.

In the following subsection, we will review all the relevant observables and
provide precise details to support these statements.

###  3.2 Observables

The $\mathbb{Z}_{4}$ model provides a concrete framework to explore the
interplay between the VLF and the ALP contributions. As discussed above, we
treat the $\mathrm{U}(1)_{\text{FN}}$-breaking parameter
$\lambda^{\prime}_{N}$ as a free parameter and consider a broad range of ALP
masses, temporarily setting aside the expectation of $m_{a}\sim v_{\Phi}$ for
$\lambda^{\prime}_{N}\sim 1$. By allowing $\lambda^{\prime}_{N}\ll 1$, we can
explore the ALP phenomenology in a concrete model. As we will argue later, the
results should remain qualitatively consistent across the entire class of FN
ALP models. Our presentation in this subsection is organized by observables,
highlighting contributions from each source.151515For related phenomenological
ALP studies see Bjorkeroth:2018dzu ; Bauer:2019gfk ; MartinCamalich:2020dfe ;
Bauer:2021mvw ; Bonilla:2022qgm ; Cornella:2023kjq ; Biekotter:2023mpd ;
Calibbi:2024rcm ; DiLuzio:2024jip ; Knapen:2023zgi .

#### Rare meson decays

_ALP._ In the $\mathrm{U}(1)_{\text{FN}}$ limit, $a$ becomes very light,
making the “invisible” meson decay $K\rightarrow\pi a$ the most significant
constraint on $v_{\Phi}$, followed by $B\rightarrow K^{(*)}a$. The decay rate
is given by

$\displaystyle\begin{aligned} \Gamma(M_{1}\rightarrow
M_{2}\,a)=\frac{m_{M_{1}}^{3}}{64\pi}\lambda^{1/2}\left(\frac{m_{M_{2}}}{m_{M_{1}}},\frac{m_{a}}{m_{M_{1}}}\right)&\left(\frac{1-m_{M_{2}}^{2}/m_{M_{1}}^{2}}{m_{q_{i}}-m_{q_{j}}}\right)^{2}\left|f_{M_{1}\rightarrow
M_{2}}(m_{a}^{2})\right|^{2}\\\
&\times\left[|c^{f}_{ij}|^{2}+|c^{f}_{ji}|^{2}-2\,\text{Re}\,c^{f}_{ij}c^{f}_{ji}\right],\end{aligned}$
(23)

where $\lambda(x,y)=(1-(x+y)^{2})(1-(x-y)^{2})$ and $q_{i},q_{j}$ refer to
quarks from the parton-level transition. The quantity $f_{M_{1}\rightarrow
M_{2}}(m_{a}^{2})$ is the form factor parameterizing the hadronic matrix
element evaluated at $q^{2}=m_{a}^{2}$, where
$\braket{M_{2}(p_{2})}{\bar{q}_{i}q_{j}}{M_{1}(p_{1})}=f_{M_{1}\rightarrow
M_{2}}(q^{2})(m_{M_{1}}^{2}-m_{M_{2}}^{2})/(m_{q_{i}}-m_{q_{j}})$ with
$q=p_{1}-p_{2}$. For these, we rely on their determinations as reported in
Refs. FlavourLatticeAveragingGroupFLAG:2021npn ; Gubernari:2023puw ;
Parrott:2022rgu ; Bailey:2015dka ; Bouchard:2013eph ; Horgan:2013hoa ;
Horgan:2015vla . The constraints we consider include the combined E787 and
E949 analysis $\text{Br}(K^{+}\rightarrow\pi^{+}E_{\text{miss}})<9.5\times
10^{-11}$ E949:2007xyy , the BaBar bound $\text{Br}(B\rightarrow
K^{*}E_{\text{miss}})<1.0\times 10^{-4}$ BaBar:2013npw  and the recent Belle
II measurement $\text{Br}(B^{+}\rightarrow K^{+}\bar{\nu}\nu)=(2.3\pm
0.7)\times 10^{-5}$ Belle-II:2023esi . For the latter, given the current
$\approx 2.4\sigma$ tension between the SM prediction and the experimental
value, we employ the $3\sigma$ bound on $\text{Br}(B\rightarrow Ka)$ derived
in Altmannshofer:2023hkn  (see also Bolton:2024egx ). For $m_{a}<m_{K,B}$,
these observables provide the most significant bounds on the FN scale, with
$v_{\Phi}\gtrsim 10^{11}$ GeV for kaon decay and $v_{\Phi}\gtrsim 10^{7}$ GeV
for $B$ decay.

Importantly, these bounds apply only if the ALP behaves as an invisible final
state. This requires that the distance traveled by the ALP after production
and before an eventual decay exceeds the length of the detector,
$l_{a}>l_{\text{detector}}$, where

$\displaystyle l_{a}=\frac{c|p_{a}|}{m_{a}\Gamma_{a}}.$ (24)

In this expression, $|p_{a}|$ represents the axion momentum in the laboratory
frame, and $\Gamma_{a}$ is its total decay width. We approximate
$l_{\text{detector}}\approx 1$ m and assume the decaying meson is at
rest.161616These considerations also apply to $M_{1}\rightarrow M_{2}\,\rho$
decays for a very light $\rho$, with the replacements $m_{a}\rightarrow
m_{\rho}$ and $c^{f}_{ij}\rightarrow-ic^{f}_{ij}$. This applies to all
processes discussed in this section; hence, we will report explicit
expressions involving only $a$, particularly since a light $\rho$ is less
motivated.

If the ALP can decay inside the detector, the signatures differ. Prominent
constraints in our model involve leptonic and semi-leptonic decays of flavored
mesons Cornella:2019uxs . The rate for lepton-flavor-violating decays is
straightforward and reads

$\displaystyle\begin{aligned}
\Gamma\left(M\rightarrow\ell_{i}\bar{\ell}_{j}\right)=\frac{f_{M}^{2}m_{M}}{64\pi}&\frac{\lambda^{1/2}\left(\frac{m_{\ell_{i}}}{m_{M}},\frac{m_{\ell_{j}}}{m_{M}}\right)}{(m_{q_{k}}+m_{q_{l}})^{2}}\frac{\left[|c^{q}_{kl}|^{2}+|c^{q}_{lk}|+2\,\text{Re}\,c^{q}_{kl}c^{q}_{lk}\right]}{\left(1-m_{a}^{2}/m_{M}^{2}\right)^{2}+m_{a}^{2}\Gamma_{a}^{2}/m_{M}^{4}}\\\
&\times\left[\left(|c^{e}_{ij}|^{2}+|c^{e}_{ji}|^{2}\right)\left(1-\frac{m_{\ell_{i}}^{2}+m_{\ell_{j}}^{2}}{m_{M}^{2}}\right)+4\,\text{Re}\,c^{e}_{ij}c^{e}_{ji}\frac{m_{\ell_{i}}m_{\ell_{i}}}{m_{M}^{2}}\right],\end{aligned}$
(25)

where $q=u,d$. Subscripts $q_{k}$ and $q_{l}$ refer to the underlying quark
transition, $\Gamma_{a}$ is the total width of $a$ and
$\braket{0}{\bar{q}_{k}\gamma^{5}q_{l}}{M}=f_{M}m_{M}^{2}/(m_{q_{k}}+m_{q_{l}})$.
We take the $95\%$ CL bounds on the branching ratios from Table 2 of
Plakias:2023esq . As expected, the strongest constraint comes from
$\text{Br}(K_{L}\rightarrow e\mu)<6.3\times 10^{-12}$
ParticleDataGroup:2022pth , which plateaus for $m_{a}\lesssim m_{K_{L}}$,
setting $v_{\Phi}\gtrsim 10$ TeV. Bounds from other scalar meson decays are at
least two orders of magnitude less stringent, while vector meson decays are
negligible due to their large width. For completeness, our numerical analysis
includes constraints from unflavored meson decays. The expressions differ
slightly due to the anomaly contribution Cornella:2019uxs . However, the large
width of these mesons leads to much weaker bounds.

Lepton flavor conserving decays also play an important role, as the diagonal
couplings in Eq. (19) are larger than the off-diagonal ones, which are further
suppressed by the small hierarchical mixing. The expression for the total
decay rate is complicated by the interference between the SM and ALP
contributions. To handle this properly, we compute the corresponding bounds
using flavio Straub:2018kue , leveraging the fact that in these processes, the
virtual momentum of the ALP is essentially fixed by the meson mass. Our
findings show that the strongest constraint arises from
$K_{L}\rightarrow\mu\mu$, which plateaus at roughly $v_{\Phi}\gtrsim 70$ TeV.

Semileptonic decays are enhanced in the on-shell ALP production regime
$m_{M_{1}}-m_{M_{2}}>m_{a}>m_{\ell_{i}}+m_{\ell_{j}}$, in which the branching
ratio can be factorized as

$\displaystyle\text{Br}(M_{1}\rightarrow
M_{2}^{(*)}\,\ell_{i}\bar{\ell}_{j})\simeq\text{Br}(M_{1}\rightarrow
M_{2}^{(*)}\,a)\text{Br}(a\rightarrow\ell_{i}\bar{\ell}_{j}),$ (26)

where

$\displaystyle\begin{aligned} \Gamma(a\rightarrow
f_{i}\bar{f}_{j})=m_{a}\frac{N_{c}}{16\pi}&\lambda^{1/2}\left(\frac{m_{i}}{m_{a}},\frac{m_{j}}{m_{a}}\right)\\\
&\times\left[\left(|c_{ij}^{f}|^{2}+|c_{ji}^{f}|^{2}\right)\left(1-\frac{m_{i}^{2}+m_{j}^{2}}{m_{a}^{2}}\right)+4\,\text{Re}\,c_{ij}^{f}c_{ji}^{f}\frac{m_{i}m_{j}}{m_{a}^{2}}\right].\end{aligned}$
(27)

Defining
$\braket{M_{2}^{*}(p_{2})}{\bar{q}_{k}\gamma^{5}q_{l}}{M_{1}(p_{1})}=-2if^{*}(q^{2})_{M_{1}\rightarrow
M_{2}^{*}}(\epsilon\cdot q)m_{M_{2}^{*}}/(m_{q_{k}}+m_{q_{l}})$ with
$q=p_{1}-p_{2}$ and $\epsilon_{\mu}=\epsilon_{\mu}(p_{2})$ the polarization
vector of $M_{2}^{*}$, one finds

$\displaystyle\begin{aligned} \Gamma(M_{1}\rightarrow
M_{2}^{*}a)=\frac{m_{M_{1}}^{3}}{64\pi}\lambda^{3/2}\left(\frac{m_{M_{2}^{*}}}{m_{M_{1}}},\frac{m_{a}}{m_{M_{1}}}\right)&\frac{|f^{*}_{M_{1}\rightarrow
M_{2}^{*}}(m_{a}^{2})|^{2}}{(m_{q_{k}}+m_{q_{l}})^{2}}\\\
&\times\left[|c_{kl}^{f}|^{2}+|c_{lk}^{f}|^{2}+2\,\text{Re}\,c_{kl}^{f}c_{lk}^{f}\right].\end{aligned}$
(28)

In this category, we can identify two classes of experimental searches.
Searches for prompt decays reconstruct the 3-body vertex with a typical
resolution of approximately $\mathcal{O}(100)\,\mu$m. When applying
constraints from such searches, we require $l_{a}\lesssim 100\,\mu\text{m}$,
unless stated otherwise. Searches for displaced vertices, instead, provide
bounds based on the lifetime of the long-lived particle. We will distinguish
between these two types of searches and their corresponding bounds in the
following discussion.

Let us first consider prompt decays. For lepton flavor-violating ones, we
revisit the processes studied in Ref. Plakias:2023esq . The strongest
constraint comes from $\text{Br}(B\rightarrow K\mu\tau)<7.7\times 10^{-6}$
Belle:2022pcr , followed by $\text{Br}(B\rightarrow\pi\mu\tau)<5.9\times
10^{-5}$ BaBar:2012azg , and $\text{Br}(B\rightarrow K^{*}\mu\tau)<9.8\times
10^{-6}$ LHCb:2022wrs . For lepton flavor conserving decays, the situation is
complicated by current tensions between SM predictions and experimental values
of various branching ratios and angular observables involving the underlying
$b\to s\mu\mu$ transition Greljo:2022jac ; Alguero:2023jeh ; Ciuchini:2022wbq
; Hurth:2023jwr ; Guadagnoli:2023ddc . In contrast, the theoretically cleaner
$R_{K^{(*)}}$ ratios are measured to be consistent with SM expectations
LHCb:2022qnv ; LHCb:2022vje . Here, we only consider these ratios to derive
constraints on the ALP contributions, roughly setting $v_{\Phi}\gtrsim$ TeV.

For displaced vertex signatures, we consider the bounds from LHCb on
$B\rightarrow K^{(*)}\mu\mu$ LHCb:2015nkv ; LHCb:2016awg , which sets
exclusion limits in the region $10^{3}\text{ GeV}\lesssim v_{\Phi}\lesssim
10^{7}$ GeV, see also Ref. Gavela:2019wzg . Additionally, we consider the
decays $B\rightarrow K\gamma\gamma$ and $K\rightarrow\pi\gamma\gamma$,
computed analogously using the well-known expressions for
$\text{Br}(a\rightarrow\gamma\gamma)$ found in Bauer:2017ris . We employ the
bound on $\text{Br}(K\rightarrow\pi\gamma\gamma)$ for long-lived ALPs from
NA64 NA62:2023olg  and the bound on $\text{Br}(B\rightarrow K\gamma\gamma)$
from BaBar BaBar:2021ich . These constraints are similar, setting
$v_{\Phi}\gtrsim$ TeV, with the former being slightly stronger but applicable
over a much narrower interval of ALP masses.

_VLF_. Let us now turn to the UV contributions obtained by integrating out
VLFs to generate local operators that directly affect rare meson decays. As
mentioned earlier, the leading effects come from operators not involving the
$\Phi$ field. In particular, the Higgs-current operators

$\begin{split}\mathcal{O}_{Hd}^{ij}&=(H^{\dagger}i\overset{\leftrightarrow}{D}_{\mu}H)(\bar{d}_{R}^{i}\gamma^{\mu}d_{R}^{j})\,,\\\
\mathcal{O}_{Hu}^{ij}&=(H^{\dagger}i\overset{\leftrightarrow}{D}_{\mu}H)(\bar{u}_{R}^{i}\gamma^{\mu}u_{R}^{j})\,,\end{split}$
(29)

are generated already at the tree level integrating out the heavy $Q_{2}^{a}$.
The corresponding Wilson coefficients read

$C_{Hd}^{ij}=\frac{y_{j}^{da}y_{i}^{da\ast}}{2M_{Q_{2}}^{2}}\,,\qquad
C_{Hu}^{ij}=\frac{y_{j}^{ua}y_{i}^{ua\ast}}{2M_{Q_{2}}^{2}}\,.$ (30)

These operators induce flavor-changing $Z$-boson couplings to right-handed
quarks after electroweak symmetry breaking Jenkins:2017jig . Focusing on down
quark transitions, we define the weak effective Hamiltonian for $d_{j}\to
d_{i}\ell\ell$ and $d_{j}\to d_{i}\nu\nu$ as
$\mathcal{H}_{\mathrm{eff}}=-\mathcal{N}\sum_{i}C_{i}O_{i}+\text{h.c.}$, where
$\mathcal{N}=\frac{4G_{F}}{\sqrt{2}}\frac{e^{2}}{16\pi^{2}}V_{tj}V_{ti}^{\ast}$
and the sum runs over the semileptonic local operators. Then, the
$\mathcal{O}_{Hd}$ operator is matched by a $Z$ exchange to

$C_{10}^{\prime ij\ell\ell}=\frac{1}{2\mathcal{N}}C_{Hd}^{ij}\,,\qquad
C_{R}^{ij\nu\nu}=\frac{1}{\mathcal{N}}C_{Hd}^{ij}\,,$ (31)

with $O_{10}^{\prime
ij\ell\ell}=(\bar{d}_{i}\gamma_{\mu}P_{R}d_{j})(\bar{\ell}\gamma^{\mu}\gamma_{5}\ell)$171717The
related operator involving a vector leptonic current $O_{9}^{\prime}$ is
suppressed due to the small SM $Z$ couplings to these leptonic currents. and
$O_{R}^{ij\nu\nu}=(\bar{d}_{i}\gamma_{\mu}P_{R}d_{j})(\bar{\nu}\gamma_{\mu}(1-\gamma_{5})\nu)$.
We stress that these contributions are lepton flavor universal. We use flavio
Straub:2018kue  to obtain the bounds from meson decays with charged or neutral
leptons in the final states. The leading constraints are due to kaon decays,
in particular to $K\to\pi\nu\nu$ setting the leading constraint
$v_{\Phi}\gtrsim 200~{}\mathrm{GeV}$. Purely leptonic kaon decays $K_{L}\to
ee,\mu\mu$ offer a slightly worse sensitivity, $v_{\Phi}\gtrsim
100~{}\mathrm{GeV}$. $B$ and $D$ decays lead to subleading constraints.

Lepton flavor-violating contributions can arise in two ways. The first is via
the exchange of a $Z$ boson, with one insertion of Eq. (29) and one of

$\displaystyle\begin{aligned}
\mathcal{O}_{H\ell}^{(1)ij}&=(H^{\dagger}i\overset{\leftrightarrow}{D}_{\mu}H)(\bar{\ell}_{R}^{i}\gamma^{\mu}\ell_{R}^{j})\\\
\mathcal{O}_{H\ell}^{(3)ij}&=(H^{\dagger}i\overset{\leftrightarrow}{D}_{\mu}^{I}H)(\bar{\ell}_{R}^{i}\tau^{I}\gamma^{\mu}\ell_{R}^{j})\,,\end{aligned}$
(32)

generated at tree level via the heavy $E_{2}^{a}$. The corresponding Wilson
coefficients are

$C_{H\ell}^{(1)ij}=C_{H\ell}^{(3)ij}=-\frac{y_{i}^{ea}y_{j}^{ea\ast}}{4M_{E}^{2}}\,.$
(33)

The second way involves integrating out the VLFs at 1-loop, directly
generating the 4-fermion SMEFT operator

$\displaystyle\mathcal{O}^{ijkl}_{\ell
d}=(\bar{\ell}_{i}\gamma^{\mu}\ell_{j})(\bar{d}_{k}\gamma_{\mu}d_{l})\,,$ (34)

with Wilson coefficient181818We used Matchete Fuentes-Martin:2022jrf  to
obtain this and the other 1-loop Wilson coefficients presented in the paper,
cross-checking the expressions with known results in the literature.

$\displaystyle\begin{aligned} C^{ijkl}_{\ell
d}=\frac{1}{128\pi^{2}}&\left[2y_{i}^{ea}(y_{j}^{ea})^{*}(y_{k}^{da})^{*}y_{l}^{da}\frac{\log
M_{E_{2}}^{2}/M_{Q_{2}}^{2}}{M_{E_{2}}^{2}-M_{Q_{2}}^{2}}\right.\\\
&\left.-\frac{3}{M_{E_{2}}^{2}}y_{i}^{ea}(y_{j}^{ea})^{*}(z_{k}^{d})^{*}z_{l}^{d}-\frac{3}{M_{Q_{2}}^{2}}z_{i}^{ea}(z_{j}^{ea})^{*}(y_{k}^{da})^{*}y_{l}^{da}\right]\,.\end{aligned}$
(35)

Despite the 1-loop factor, this provides the dominant constraint since
$M_{Q_{2},E_{2}}\gtrsim 4\pi v$. In principle, there is a third way via SMEFT
(renormalization group) RG effects on Eq. (32). However, the RG effect is
always subleading, as flavor mixing would be CKM-suppressed and thus less
significant compared to the anarchic nature of the coefficient in Eq. (35). To
establish the bounds, we map the contribution of these Wilson coefficients to
the same set of lepton flavor violating decays analyzed in the ALP-mediated
case, following Plakias:2023esq . The resulting constraints are significantly
weaker compared to the lepton flavor-conserving ones obtained earlier, with
the strongest being associated with $K_{L}\rightarrow e\mu$, leading to
$v_{\Phi}\gtrsim 10$ GeV.

#### Neutral meson mixing

One of the leading constraints for both flavorful ALP and VLQs is due to their
contributions to the CP-even and CP-odd observables measured in
$P^{0}-\bar{P}^{0}$ oscillations, with $P=K,B,B_{s},D$. We use flavio to
obtain constraints on both the ALP and the VLQ contributions to the following
observables: $\epsilon_{K}$, measuring CP violation in kaon mixing
Buras:2010pza ; Blum:2011ng ; Brod:2011ty ; Garron:2016mva ;
ParticleDataGroup:2022pth ; FlavourLatticeAveragingGroupFLAG:2021npn
,191919Note that we do not use $\Delta m_{K}$ as a constraint due to the large
uncertainties in its SM prediction from long-distance contributions. Instead,
we use its experimental value as an input when predicting $\epsilon_{K}$.
Nevertheless, we checked that saturating the experimental value of $\Delta
m_{K}$ with the NP contributions leads at most to bounds comparable to that
from $\epsilon_{K}$. $\Delta m_{s}$ ($\Delta m_{d}$) and $S_{\psi\phi}$
($S_{\psi K}$) measuring the neutral meson mass difference and mixing induced
CP violation in $B_{s}$ ($B_{d}$) mixing King:2019lal ; Dowdall:2019bea ;
Greljo:2022jac ; HFLAV:2022esi , and $x_{12}^{\mathrm{Im},D}$ measuring CP
violation in $D$ mixing Carrasco:2015pra ; HFLAV:2022esi .

_ALP_. We write the $\Delta F=2$ Hamiltonian in the quark mass basis as

$\displaystyle\begin{aligned} \mathcal{H}^{\Delta
F=2}_{\mathrm{NP}}&=C_{1}^{ij}(\bar{q}_{L}^{i}\gamma^{\mu}q_{L}^{j})^{2}+\widetilde{C}_{1}^{ij}(\bar{q}_{R}^{i}\gamma^{\mu}q_{R}^{j})^{2}+C_{2}^{ij}(\bar{q}_{R}^{i}q_{L}^{j})^{2}+\widetilde{C}_{2}^{ij}(\bar{q}_{L}^{i}q_{R}^{j})^{2}\\\
&+C_{4}^{ij}(\bar{q}_{R}^{i}q_{L}^{j})(\bar{q}_{L}^{i}q_{R}^{j})+C_{5}^{ij}(\bar{q}_{L}^{i}\gamma^{\mu}q_{L}^{j})(\bar{q}_{R}^{i}\gamma_{\mu}q_{R}^{j})+\text{h.c.}\,,\end{aligned}$
(36)

with $q=u,d$ and $i,j$ are flavor indices, whereas the color indices are
contracted inside parentheses. Integrating out the scalar and pseudoscalar
components of $\Phi$ generates $\Delta F=2$ operators already at the tree
level as follows

$\displaystyle\begin{aligned}
C_{2}^{q,ij}&=-\frac{(c_{ji}^{q*})^{2}}{2}\left(\frac{1}{m_{\rho}^{2}}-\frac{1}{m_{a}^{2}}\right)\,,\\\
\widetilde{C}_{2}^{q,ij}&=-\frac{(c_{ij}^{q})^{2}}{2}\left(\frac{1}{m_{\rho}^{2}}-\frac{1}{m_{a}^{2}}\right)\,,\\\
C_{4}^{q,ij}&=-c_{ij}^{q}c_{ji}^{q*}\left(\frac{1}{m_{\rho}^{2}}+\frac{1}{m_{a}^{2}}\right)\,.\\\
\end{aligned}$ (37)

Here we have assumed the ALP to be heavier than the appropriate scale of each
$P-\bar{P}$ process. In the case of a light ALP, this is no longer applicable,
and ALP propagator effects should be included Bauer:2021mvw . However, in this
work, we do not consider the meson mixing constraints in the case of a light
ALP, as we find that parameter space is predominantly constrained by other
processes, such as rare meson decays. Outside of that region, we anticipate
that meson mixing observables are responsible for the strongest bound coming
from ALP physics. This roughly sets $v_{\Phi}\gtrsim 10\text{
TeV}\times(10\text{ GeV}/m_{a})$.

_VLF_. Among the VLQs with the interaction Lagrangian defined in Eq. (12),
$Q_{2}^{a}$ can generate $\Delta F=2$ operators via two effects. The first
effect is obtained by integrating out $Q_{2}^{a}$ at the tree level, as
discussed in the rare meson decay case; see Eq. (30). The generated operators
could, in principle, generate $\Delta F=2$ operators via two insertions of the
modified $Z$ couplings. Such contributions are highly suppressed. Important
effects can, however, be generated by 1-loop SMEFT RG. Performing one
insertion of (30) and closing the loop with two Higgses and a fermion, we
obtain $y_{t}$-enhanced effects by generating the left-right 4-fermion
$\mathcal{O}_{qd}$ operator (in the notation of Grzadkowski:2010es ). These
effects have been found to be phenomenologically relevant in Ref.
Bobeth:2016llm , and we take them into account in our analysis.

The second way to generate $\Delta F=2$ effects is from one-loop matching
contributions. This directly generates 4-fermion operators required to induce
such effects:

$\displaystyle\begin{aligned}
\mathcal{O}_{dd}^{ijkl}&=(\bar{d}_{i}\gamma^{\mu}d_{j})(\bar{d}_{k}\gamma_{\mu}d_{l}),\\\
\mathcal{O}_{uu}^{ijkl}&=(\bar{u}_{i}\gamma^{\mu}u_{j})(\bar{u}_{k}\gamma_{\mu}u_{l}).\end{aligned}$
(38)

The Wilson coefficients associated with these operators read

$\displaystyle\begin{aligned}
C^{ijkl}_{dd}&=\frac{1}{64\pi^{2}}\frac{1}{M_{Q_{2}}^{2}}(y_{i}^{da})^{*}y_{j}^{da}\left[-(y_{k}^{db})^{*}y_{l}^{db}+3(z_{s}^{d})^{*}z_{t}^{d}\right],\\\
C^{ijkl}_{uu}&=\frac{1}{64\pi^{2}}\frac{1}{M_{Q_{2}}^{2}}(y_{i}^{ua})^{*}y_{j}^{ua}\left[-(y_{k}^{ub})^{*}y_{l}^{ub}+3(z_{s}^{u})^{*}z_{t}^{u}\right].\end{aligned}$
(39)

The $\Delta F=2$ operators are obtained by selecting $i=k$ and $j=l$. In the
minimal basis defined in App. A, the terms involving $z$ couplings do not
contribute.

Employing the benchmark discussed in App. A, we find that these one-loop
matching contributions offer the dominant constraint, which comes from
$\epsilon_{K}$, setting $v_{\Phi}\gtrsim 2.8~{}\mathrm{TeV}$. The rest of the
$\Delta F=2$ constraints, namely the CP-even and CP-odd observables in
$B_{d}$, $B_{s}$ and $D$-mixing, show various degrees of interplay between
tree-level and one-loop contributions of VLQs, but offer only subleading
constrains of $v_{\Phi}\gtrsim 100~{}\mathrm{GeV}$.

It is worth mentioning that also the operators
$\mathcal{O}_{qd},\mathcal{O}_{qu}$ can be generated at the one-loop level and
could, in principle, lead to important effects. Indeed, the components
diagonal in $q$ do not need an insertion of $\Phi$ and hence can be generated
already at dimension 6. After electroweak symmetry breaking, the rotation to
the mass basis induces additional flavor violation in the down-strange sector,
suppressed only by the CKM. This suppression could, in principle, be overcome
by an enhancement in the observables due to the different chirality structures
of the operators in such a way that the bounds become potentially important.
Accidentally, in $\mathbb{Z}_{4}$ models the operator
$(\bar{q}_{2}\gamma^{\mu}q_{2})(\bar{d}_{i}\gamma_{\mu}d_{j})$ is generated at
1-loop with coefficient $\sim(y^{da}_{i}x_{2}^{qa})^{*}(y^{da}_{j}x_{2}^{qa})$
(and similarly for up-quarks), which vanishes whenever the first generator is
involved (see App. A). Hence, this effect could only appear in $B_{s}$-mixing,
which is not very constraining.

#### $\bm{\epsilon^{\prime}/\epsilon}$

The $\epsilon^{\prime}/\epsilon$ observable acts as a probe of direct CP
violation in $K_{L}\to\pi\pi$. The current best estimate of the SM predictions
is $(\epsilon^{\prime}/\epsilon)_{\mathrm{SM}}=(13.9\pm 5.2)\times 10^{-4}$
RBC:2020kdj ; Buras:2020pjp ; Aebischer:2020jto ; Aebischer:2021hws ;
Cirigliano:2019cpi , while the experimental world average is
$(\epsilon^{\prime}/\epsilon)_{\mathrm{exp}}=(16.6\pm 2.3)\times 10^{-4}$
ParticleDataGroup:2022pth . We use the master formula presented in Refs.
Aebischer:2021hws ; Aebischer:2018quc  to compute the BSM contributions to
$\epsilon^{\prime}/\epsilon$ as

$\left(\frac{\epsilon^{\prime}}{\epsilon}\right)_{\mathrm{BSM}}=\sum_{b}P_{b}(\mu_{\mathrm{ew}})\mathrm{Im}\left[C_{b}(\mu_{\mathrm{ew}})-C_{b}^{\prime}(\mu_{\mathrm{ew}})\right]\,,$
(40)

where $P_{b}$ parameterizes the matrix elements of the relevant four quark
operators, which we consider in the JMS basis Jenkins:2017jig . We use the
numerical values reported in Table 3 of Ref. Aebischer:2021hws .

_ALP_. The ALP contributes to $\epsilon^{\prime}/\epsilon$ at the tree-level
through scalar four-fermion $\Delta F=1$ operators. In the notation of Ref.
Aebischer:2021hws  these are $[O_{dd}^{S1,RR}]_{2111}$,
$[O_{ud}^{S1,RR}]_{1121}$, and their $L\leftrightarrow R$ counterparts.
However, its contributions are suppressed by first-generation quark masses. At
most, this renders only subleading constraint up to $v_{\Phi}\gtrsim 100\text{
GeV}$.

_VLF_. VLQs induce flavor-violating $Z$-boson couplings to right-handed quarks
via Eq. (30). As can be seen from Table $3$ of Ref. Aebischer:2021hws , the
largest matrix elements are due to left-right vector operators. In our setup,
$[O_{dd}^{V1,LR}]_{2111}$, $[O_{du}^{V1,LR}]_{2111}$ are generated by
integrating out the $Z$ boson, with VLQ-induced flavor violating right-handed
down quark vector currents, contracted with either left-handed up or down
quark currents due to the SM $Z$ couplings. We obtain the bound
$v_{\Phi}\gtrsim 0.7~{}\text{TeV}$.

#### $\bm{\ell_{i}\to\ell_{j}\gamma}$

Both the ALP and the VLFs generate lepton flavor-violating couplings that are
strongly constrained by purely leptonic observables. Dipole transitions can be
encoded in the WET Lagrangian via

$\mathcal{L}_{D}\supset
C_{e\gamma}^{ij}\left(\bar{e}_{i}\sigma^{\mu\nu}P_{L}e_{j}\right)F^{\mu\nu}+\mathrm{h.c.}\,.$
(41)

Assuming $m_{\ell_{j}}>m_{\ell_{i}}$, the branching ratio of
$\ell_{j}\to\ell_{i}\gamma$ then reads Lavoura:2003xp

$\mathrm{Br}(\ell_{j}\to\ell_{i}\gamma)=\frac{m_{\ell_{i}}^{3}}{4\pi\Gamma_{\ell_{j}}}\left(|C_{e\gamma}^{ij}|^{2}+|C_{e\gamma}^{ji}|^{2}\right)\,.$
(42)

_ALP_. The ALP $a$ generates these observables at the one loop. In the regime
in which the ALP is much more massive than $\ell_{j}$, the Wilson coefficient
associated with the operator (41) reads

$\begin{split}C_{e\gamma}^{ij}&=\frac{e}{32\pi^{2}m_{a}^{2}}\sum_{k=1,2,3}\left[\frac{1}{6}\left(m_{\ell_{j}}c_{ki}^{e\ast}c_{kj}^{e}+m_{\ell_{i}}c_{ik}^{e}c_{jk}^{e\ast}\right)-m_{\ell_{k}}c_{ki}^{e\ast}c_{jk}^{e\ast}\left(\frac{3}{2}+\log{\frac{m_{\ell_{j}}^{2}}{m_{a}^{2}}}\right)\right]\,\end{split}$
(43)

whereas the contribution of $\rho$ can be obtained by the replacement
$m_{a}^{2}\to m_{\rho}^{2}$, $c\to-ic$. The contribution from Barr-Zee
diagrams is relevant only for ALP masses below the charm mass, given the
smallness of our ALP couplings with the top and bottom quarks. The expression
for these contributions can be found in Bauer:2021mvw , and we take them into
account in our final plot. The relevant current experimental bounds are
$\text{Br}(\mu\rightarrow e\gamma)<4.2\times 10^{-13}$ MEG:2013oxv ;
MEG:2016leq  and $\text{Br}(\tau\rightarrow\mu\gamma)<4.4\times 10^{-8}$
BaBar:2009hkt ; Belle:2007qih . Both lead to comparable bounds due to the
hierarchical structure of $c^{e}_{ij}$ up to $v_{\Phi}\gtrsim$ TeV.

_VLF_. The operator (41) can also be generated by integrating out the heavy
VLFs at 1-loop. The corresponding $\mathbb{Z}_{4}$-invariant SM$+\Phi$ terms
are

$\displaystyle
C_{eB\Phi}^{ij}\,\Phi^{n^{e}_{ij}}\bar{\ell}_{i}\sigma^{\mu\nu}He_{j}B_{\mu\nu}+C_{eW\Phi}^{ij}\,\Phi^{n^{e}_{ij}}\bar{\ell}_{i}\sigma^{\mu\nu}\tau^{a}He_{j}W_{\mu\nu}^{a}\,,$
(44)

where by simple dimensional analysis $C_{eB\Phi,eW\Phi}^{ij}\sim
g^{(\prime)}y_{*}^{3+n^{e}_{ij}}/16\pi^{2}M^{2+n^{e}_{ij}}$, with $y_{*}$ some
representative UV coupling. After $\Phi$ condensation, this leads to the SMEFT
operators $\mathcal{O}_{eB,eW}$, which combine to Eq. (41). Using Matchete
Fuentes-Martin:2022jrf , we were able to exactly extract the Wilson
coefficients relevant for $\mu\rightarrow e\gamma$ and
$\tau\rightarrow\mu\gamma$:

$\displaystyle\begin{aligned}
C_{eB\Phi}^{i2}&=\frac{g^{\prime}}{32\pi^{2}}\frac{1}{M_{E_{2}}^{2}}\left[-\frac{1}{M_{E_{1}}}(x^{ea}_{12})^{*}x_{2}^{ea}x_{12}^{eb}y^{eb}_{i}\right.\\\
&\qquad\qquad\qquad\left.+\frac{1}{12}\frac{1}{M_{E_{2}}}\left(x_{2}^{ea}(x_{2}^{ea})^{*}x_{2}^{eb}y_{i}^{eb}+2x_{2}^{ea}y_{k}^{ea}(y_{k}^{eb})^{*}y_{i}^{eb}\right)\right],\\\
C_{eW\Phi}^{i2}&=\frac{g}{384\pi^{2}}\frac{1}{M_{E_{2}}^{3}}x_{2}^{ea}y_{k}^{ea}(y_{k}^{eb})^{*}y_{i}^{eb}\,.\end{aligned}$
(45)

The coefficients $C_{eB\Phi,eW\Phi}^{i3}$ are null because $e_{3}$ does not
couple to any VLF, while $C_{eB\Phi,eW\Phi}^{i1}$ arise at dimension 8 since
they involve a further power of $\Phi$. Hence, $C_{eB\Phi,eW\Phi}^{i2}$
provide the leading contributions to the aforementioned transitions. Note that
the couplings in Eq. (45) are not perfectly aligned with the effective $Y_{e}$
(see App. A). This leads to a further mixing between the components of the
dipoles after the diagonalization of the Yukawa with a $e_{i}$ rotation. We
match these coefficient to $C_{e\gamma}$ after $\Phi$ condensation via
$C_{e\gamma}^{ij}=v(c_{w}C_{eB}-s_{w}C_{eW})/\sqrt{2}$. The coefficient
qualitatively scales as $C_{e\gamma}^{i2}\sim
evy_{*}^{4}\epsilon^{3}/16\pi^{2}v_{\Phi}^{2}$, that can lead to bounds on
$v_{\Phi}$ of $\mathcal{O}(\text{TeV})$ for $y_{*}\sim\mathcal{O}(1)$ due to
the stringent constraint on $\mu\rightarrow e\gamma$. Precisely, employing the
benchmark in App. A, we find $v_{\Phi}\gtrsim 150$ GeV.

####  $\bm{\mu\to e}$ conversion in nuclei

$\mu\rightarrow e$ conversion in nuclei is one of the most constraining
observables on new physics models predicting lepton flavor violation.
Following Ref. Crivellin:2017rmk , the operators relevant to our discussion
comprise the following set

$\begin{split}O_{qq}^{VLL}&=(\bar{e}\gamma_{\mu}P_{L}\mu)(\bar{q}\gamma^{\mu}P_{L}q)\,,\\\
O_{qq}^{VLR}&=(\bar{e}\gamma_{\mu}P_{L}\mu)(\bar{q}\gamma^{\mu}P_{R}q)\,,\\\
O_{qq}^{SLL}&=(\bar{e}P_{L}\mu)(\bar{q}P_{L}q)\,,\\\
O_{qq}^{SLR}&=(\bar{e}P_{L}\mu)(\bar{q}P_{R}q)\,,\\\
O_{gg}^{L}&=\alpha_{s}m_{\mu}G_{F}(\bar{e}P_{L}\mu)G_{\mu\nu}^{a}G^{a\mu\nu}\,,\end{split}$
(46)

together with their $L\leftrightarrow R$ equivalents and the dipole operators
of Eq. (41).202020As for $\ell_{i}\rightarrow\ell_{j}\gamma$ processes, the
contributions from Barr-Zee diagrams are relevant only for $m_{a}\lesssim
m_{c}$. We take this into account in our final plot, taking the expression for
these contributions from Bauer:2021mvw . The conversion rate on a nucleus $N$
is Kitano:2002mt ; Cirigliano:2009bz ; Crivellin:2017rmk ; Shifman:1978zn

$\Gamma_{\mu\to
e}^{N}=\frac{m_{\mu}^{5}}{4}\left|C_{L}^{D}D_{N}+2\left(G_{F}m_{\mu}m_{p}\tilde{C}_{(p)}^{SL}S_{N}^{(p)}+\tilde{C}_{(p)}^{VR}V_{N}^{(p)}+p\to
n\right)\right|^{2}+L\leftrightarrow R\,,$ (47)

where $C_{L}^{D}\equiv C_{e\gamma}^{e\mu}/m_{\mu}$ and

$\begin{split}\tilde{C}_{p/n}^{VR}&=\sum_{q=u,d}\left(C_{qq}^{VRL}+C_{qq}^{VRR}\right)f_{Vp/n}^{(q)}\,,\\\
\tilde{C}_{(p/n)}^{SL}&=\sum_{q=u,d,s}\frac{\left(C_{qq}^{SLL}+C_{qq}^{SLR}\right)}{m_{\mu}m_{q}G_{F}}f_{Sp/n}^{(q)}+\tilde{C}_{gg}^{L}f_{Gp/n}\,,\\\
\tilde{C}_{gg}^{L}&=C_{gg}^{L}-\frac{1}{12\pi}\sum_{q=c,b}\frac{\left(C_{qq}^{SLL}+C_{qq}^{SLR}\right)}{m_{\mu}m_{q}G_{F}}\,,\end{split}$
(48)

and similarly for $L\leftrightarrow R$. The vector operator nucleon form
factors are $f_{Vp}^{(u)}=2$, $f_{Vn}^{(u)}=1$, $f_{Vp}^{(d)}=1$,
$f_{Vn}^{(d)}=2$, while for the scalar and the gluonic form factors, we take
the values as reported in Ref. Crivellin:2017rmk , see also Refs.
Crivellin:2014cta ; Hoferichter:2015dsa ; Crivellin:2013ipa ; Junnarkar:2013ac
. The most sensitive current measurement is on $\mathrm{Au}$ for which we have
the overlap integrals $V_{\mathrm{Au}}^{(p)}=0.0974$,
$V_{\mathrm{Au}}^{(n)}=0.146$, $S_{\mathrm{Au}}^{(p)}=0.0614$,
$S_{\mathrm{Au}}^{(n)}=0.0918$, $D_{\mathrm{Au}}=0.189$ Kitano:2002mt .
Finally, the prediction should be compared to the 90% CL limit on $\mu\to e$
conversion in gold by SINDRUM II SINDRUMII:2006dvw

$\mathrm{Br}(\mu\to e)=\frac{\Gamma(\mu^{-}\mathrm{Au}\to
e^{-}\mathrm{Au})}{\Gamma_{\mathrm{capt}}(\mu^{-}\mathrm{Au})}<7\times
10^{-13}\,,$ (49)

where the estimate for the capture rate can be found in Suzuki:1987jf ;
Kitano:2002mt .

_ALP_. The ALP contributes to $\mu\to e$ conversion in multiple ways: directly
generating the scalar four fermion operators $O_{qq}^{S}$ at the tree level,
contributing through the dipole operator as in Eq. (43), or generating
$O_{gg}$ after integrating out the heavy SM quarks. The tree-level matching
results in

$\begin{split}C_{qq}^{SLL}&=c^{q\ast}_{qq}c^{e\ast}_{21}\left(\frac{1}{m_{\rho}^{2}}-\frac{1}{m_{a}^{2}}\right)\,,\\\
C_{qq}^{SLR}&=c^{q}_{qq}c^{e\ast}_{21}\left(\frac{1}{m_{\rho}^{2}}+\frac{1}{m_{a}^{2}}\right)\,,\\\
C_{qq}^{SRL}&=c^{q\ast}_{qq}c^{e}_{12}\left(\frac{1}{m_{\rho}^{2}}+\frac{1}{m_{a}^{2}}\right)\,,\\\
C_{qq}^{SRR}&=c^{q}_{qq}c^{e}_{12}\left(\frac{1}{m_{\rho}^{2}}-\frac{1}{m_{a}^{2}}\right)\,,\\\
\end{split}$ (50)

Note that for coherent conversion in Eq. (47) only the scalar quark current
operator enters, thus the pseudoscalar $a$ contributions above cancel as
$c_{qq}^{q}$ are real. We obtain bounds of up to $v_{\Phi}\gtrsim$ TeV.

_VLF_. Muon conversion into an electron in the field of a nucleus can also
proceed through SM $Z$ exchange, with flavor-violating couplings to leptons
induced by the operator in Eq. (32) and exploiting the SM coupling to quarks.
The vector operators in Eq. (46) are generated with Crivellin:2020ebi

$\displaystyle\begin{aligned} C_{qq}^{VLL}&=\Gamma_{e\mu}^{\ell
L}\frac{1}{M_{Z}^{2}}\Gamma_{qq}^{L}\,,\\\ C_{qq}^{VLR}&=\Gamma_{e\mu}^{\ell
L}\frac{1}{M_{Z}^{2}}\Gamma_{qq}^{R}\,.\end{aligned}$ (51)

The $Z$ boson couples to quarks via $\Gamma_{qq}^{L,R}$ in the standard way,
whereas the $Z\mu e$ coupling reads

$\Gamma_{e\mu}^{\ell
L}=\frac{g_{2}}{2c_{W}}v^{2}\left(C_{H\ell}^{(1)12}+C_{H\ell}^{(3)12}\right)\,.$
(52)

Using our benchmark, we obtain the bound $v_{\Phi}\gtrsim 1.1~{}\mathrm{TeV}$.

#### Rare lepton decays

Rare lepton decays comprise a wide set of observables, including flavor-
violating $\ell_{i}\rightarrow\ell_{j}+\text{invisible}$,
$\ell_{i}\rightarrow\ell_{j}\ell_{k}\ell_{k}$ and $\ell_{i}\rightarrow
M^{(*)}\ell_{j}$. As usual, we will focus on ALP and VLF contributions
separately.

_ALP_. As with mesons, we begin by discussing flavor-violating decays with an
invisible ALP as a final state. Bounds on the decay of a muon into an electron
and an invisible scalar have been set at TRIUMF Jodidio:1986mz  and from TWIST
TWIST:2014ymv , while analogous bounds for $\tau$ decays have been reported by
Belle II Belle-II:2022heu . The relevant formula for such decays is

$\displaystyle\Gamma(f_{i}\rightarrow
f_{j}a)=m_{i}\frac{N_{c}}{32\pi}\lambda^{1/2}\left(\frac{m_{a}}{m_{i}},\frac{m_{j}}{m_{i}}\right)\left[\left(|c_{ij}^{f}|^{2}+|c_{ji}^{f}|^{2}\right)\left(1-\frac{m_{a}^{2}-m_{j}^{2}}{m_{i}^{2}}\right)-4\,\text{Re}\,c_{ij}^{f}c_{ji}^{f}\frac{m_{j}}{m_{i}}\right],$
(53)

and we again impose $\ell_{a}\gtrsim 1$ m for the ALP to be considered
invisible. The bounds on $v_{\Phi}$ that we obtain are almost always
superseded by the one coming from $K\rightarrow\pi a$.

Next we consider the $\ell_{i}\to M^{(*)}\ell_{j}$ processes, where $M$ is a
scalar (vector) meson. In the ALP-mediated case, the decay rate for decays in
flavored scalar mesons closely resembles one of the leptonic flavored meson
decays (25):

$\displaystyle\begin{aligned} \Gamma\left(\ell_{i}\rightarrow
M\ell_{j}\right)=\frac{f_{M}^{2}m_{\ell_{i}}}{128\pi}&\frac{\lambda^{1/2}\left(\frac{m_{M}}{m_{\ell_{i}}},\frac{m_{\ell_{j}}}{m_{\ell_{i}}}\right)}{(m_{q_{k}}+m_{q_{l}})^{2}}\frac{\left[|c^{q}_{kl}|^{2}+|c^{q}_{lk}|+2\,\text{Re}\,c^{q}_{kl}c^{q}_{lk}\right]}{\left(1-m_{a}^{2}/m_{M}^{2}\right)^{2}+m_{a}^{2}\Gamma_{a}^{2}/m_{M}^{4}}\\\
&\times\left[\left(|c^{e}_{ij}|^{2}+|c^{e}_{ji}|^{2}\right)\left(1-\frac{m_{M}^{2}-m_{\ell_{j}}^{2}}{m_{\ell_{i}}^{2}}\right)-4\,\text{Re}\,c^{e}_{ij}c^{e}_{ji}\frac{m_{\ell_{j}}}{m_{\ell_{i}}}\right]\,.\end{aligned}$
(54)

The expression for decays to vector mesons vanishes at tree level because of
the scalar nature of the interaction mediated by $a$, while for unflavored
mesons, the situation is once again complicated by the overlap with
$G\widetilde{G}$ Cornella:2019uxs . Since the experimental constraints on the
branching ratio of these decays are all quite similar, the most constraining
ones are those involving decays of taus into muons ParticleDataGroup:2022pth
and all lead to comparable bounds $v_{\Phi}\gtrsim$ TeV.

For $\ell_{i}\rightarrow\ell_{j}\ell_{k}\ell_{k}$, the branching ratio in the
presence of an ALP is enhanced in the on-shell production regime
$2m_{k}<m_{a}<m_{i}-m_{j}$. Here, the ratio can be factorized as

$\displaystyle\text{Br}(\ell_{i}\rightarrow\ell_{i}\bar{\ell}_{k}\ell_{k})\simeq\text{Br}(\ell_{i}\rightarrow\ell_{j}a)\text{Br}(a\rightarrow\bar{\ell}_{k}\ell_{k}).$
(55)

with relevant decay rates given in Eqs. (27) and (53). The situation for these
decays is similar to that of semileptonic meson decays in that we can
distinguish prompt and displaced vertex searches. Among the former, we find
that the strongest constraint is due to $\text{Br}(\tau\rightarrow
3\mu)<2.5\times 10^{-8}$ at $95\%$ CL Belle-II:2024sce  constraining up to
$v_{\Phi}\gtrsim 10^{5}$ GeV. Decays into electrons are suppressed due to the
hierarchical nature of $c^{e}$; even the strong bounds from $\mu\to 3e$
SINDRUM:1987nra  do not beat $\tau\rightarrow 3\mu$. Searches for displaced
vertices are instead lacking in this context, although new searches have been
proposed over the last years Heeck:2017xmg ; Cheung:2021mol . To our
knowledge, the only available bounds are on $\mu\rightarrow 3e$
SINDRUM:1986klz  and $\mu\rightarrow e\gamma\gamma$ Bolton:1988af ;
Natori:2012gga ; MEG:2020zxk . These, however, lead to completely negligible
bounds in our model of the order $v_{\Phi}\gtrsim$ GeV.

_VLF_. The VLLs can induce the $\ell_{i}\to M^{(*)}\ell_{j}$ decays in
multiple ways: either by a double insertion of the tree-level operators in
Eqs. (29), (32), a single insertion of the 1-loop generated operator in Eq.
(34), or, for unflavored $M$, by a single LFV modification of the $Z$ boson
couplings from Eq. (32). We find that the latter effect is phenomenologically
the most important. Using flavio Straub:2018kue  we obtain constraints from
$\tau\to\rho\ell$ and $\tau\to\phi\ell$ Aebischer:2018iyb , which are however
of order $v_{\Phi}\gtrsim 20~{}\mathrm{GeV}$ and hence subleading.

Also the $\ell\to 3\ell$ decays can proceed through a $Z$ exchange212121It is
worth mentioning that in our model the operator
$\mathcal{O}_{\ell\ell}^{ijkl}=(\bar{\ell}_{i}\gamma^{\mu}\ell_{j})(\bar{\ell}_{k}\gamma^{\mu}\ell_{l})$
is also generated at 1-loop integrating out the VLLs. This effect is, however,
subleading compared to the tree-level one. with the tree-level LFV coupling
from Eq. (32) and an SM $Z$ coupling. As the couplings in the VLL benchmark
are anarchic, $\mu\to 3e$ is now the only phenomenologically important
channel. The branching ratio prediction reads Crivellin:2020ebi

$\mathrm{Br}(\mu\to
3e)=\frac{m_{\mu}^{5}}{1536\pi^{3}m_{Z}^{4}\Gamma_{\mu}}\left(2|\Gamma_{e\mu}^{\ell
L}\Gamma_{ee}^{\ell L}|^{2}+|\Gamma_{e\mu}^{\ell L}\Gamma_{ee}^{\ell
R}|^{2}\right)\,,$ (56)

with $\Gamma_{e\mu}^{\ell L}$ defined in Eq. (52) and with the SM $Z$
couplings to electrons. We obtain the bound $v_{\Phi}\gtrsim 400\text{ GeV}$.

#### EDMs

Electric dipole moments are among the most sensitive probes of CP violating
new physics. The strongest constraints are the ones on the electron EDM,
$|d_{e}|/e<5.4\times 10^{-30}$ cm at $95\%$ CL Roussy:2022cmp , and on the
neutron EDM, $|d_{n}|/e<2.6\times 10^{-26}$ cm at $95\%$ CL Abel:2020pzs . In
particular, the electron EDM at leading order is directly associated with the
dipole operator $\mathcal{O}_{e\gamma}$ introduced earlier Kley:2021yhn

$\displaystyle d_{e}\simeq-2\,\text{Im}\,C_{e\gamma}^{11}\,,$ (57)

while the expression for the neutron EDM, in the notation of Kley:2021yhn ,
reads222222We omit the contribution from $C_{Hud}$, since in any case below
the EW scale this operator can be mapped to 4-fermion operators obtained
integrating out $W^{\pm}$ bosons.

$\displaystyle\begin{aligned}
d_{n}\simeq&-0.2\,d_{u}+0.8\,d_{d}-0.003\,d_{s}+0.05\,\hat{d}_{u}+0.1\,\hat{d}_{d}-51e\cdot\text{MeV}\,C_{G\widetilde{G}}\\\
&-0.6e\times\text{GeV}\,\text{Im}\,\left[C_{ud}^{S1,RR}+C_{ud}^{S8,RR}-C_{duud}^{S1,RR}-C_{duud}^{S8,RR}\right]^{1111}.\end{aligned}$
(58)

In principle, our model also predicts contributions to the magnetic dipole
moments, which are much less constraining than those from the EDMs.

_ALP_. In the heavy ALP case ($m_{a}\gtrsim$ GeV), the leading contribution to
the eEDM can be directly read off from Eq. (43). The imaginary part of
$C_{e\gamma}^{11}$ vanishes, meaning that no contribution to the eEDM is
generated at this order. To check that NLO corrections do not lead to
signification bounds, we can estimate the hypothetical constraint if an
$\mathcal{O}(1)$ imaginary part in $C_{e\gamma}$ was there, obtaining the weak
bound $v_{\Phi}>\mathcal{O}(10)\text{ GeV}\times(\text{GeV}/m_{a})$,
indicating that this observable is under control. For the neutron EDM, the
situation is similar: the dipoles
($d_{u},d_{d},d_{s},\hat{d}_{u},\hat{d}_{d}$) feature an expression close to
that of $C_{e\gamma}$ and are hence real. The 4-fermion operator coefficient
$C_{ud}^{S1,RR}$ entering Eq. (58) is generated at the tree level, but
features no imaginary part. An analysis analogous to the one for the eEDM
reveals that the bound for a hypothetical $\mathcal{O}(1)$ phase reads
$v_{\Phi}>\mathcal{O}(10^{2})\text{ GeV}\times(\text{GeV}/m_{a})$, once again
negligible compared to the flavor bounds analyzed earlier. In the very light
ALP case ($m_{a}\lesssim$ GeV), contributions to the nEDM approximate to a
constant close to the previous estimate when $m_{a}\sim m_{u}$ DiLuzio:2020oah
. While this could potentially lead to a bound comparable to the ones from
other observables analyzed earlier, in reality, we know that this contribution
is still not there due to the absence of an imaginary part in the relevant
Wilson coefficients. Therefore, whatever effect may arise at NLO must be
phenomenologically irrelevant. The operator $C_{G\widetilde{G}}$ can be
generated only at very high loop orders.

_VLF_. The VLFs generate the operators Eq. (44) directly contributing to the
eEDM. We can estimate the size of the bound employing Eq. (45), that after
$e_{i}$ rotation generates a non-zero entry $C_{e\gamma}^{11}\sim
evy_{*}^{5}\epsilon^{4}/16\pi^{2}v_{\Phi}^{2}$ with an $\mathcal{O}(1)$
imaginary part. This is a single contribution to the eEDM, which, in
principle, should be added to that of the dimension-8
$C_{eB\Phi,eW\Phi}^{i,1}$, whose induced $C_{e\gamma}^{11}$ features the same
scaling. However, barring accidental cancellations, this should be enough to
estimate the constraint coming from the eEDM. We find the bound
$v_{\Phi}\gtrsim 230$ GeV for our benchmark. This is interesting but again
subleading compared to $\epsilon_{K}$ and other bounds derived earlier. For
the nEDM, the bound from the dipole operators can be obtained in an analogous
way to that of the eEDM. As the experimental bound on the former is 4 orders
of magnitude weaker than the one on the latter, and given the same functional
dependence on the UV scales and couplings (but in the quark sector), this does
not lead to a competitive bound on $v_{\Phi}$. The 4-fermion operators
entering Eq. (58) are associated to $\mathcal{O}_{quqd}$ in the SMEFT, but in
the SM+$\Phi$ EFT are uplifted to
$\mathcal{O}_{quqd\Phi}^{ijkl}\sim\Phi^{n_{ik}^{q}}(\bar{q}_{i}u_{j})\epsilon(\bar{q}_{k}d_{l})$.
Due to the large $\epsilon$ suppression, this does not lead to relevant
constraints. Finally, the full expression for Eq. (58) should include also the
contribution from $C_{Hud}^{ij}=-(y_{i}^{ua})^{*}y_{j}^{da}/M_{Q_{2}}^{2}$
Kley:2021yhn . We find $\,\text{Im}\,C_{Hud}^{11}=0$; hence this cannot
contribute to leading order. A hypothetical $\mathcal{O}(1)$ imaginary part
would, in any case, result in the subleading
$v_{\Phi}>|y_{*}|\times\mathcal{O}(10^{2})$ GeV, hence NLO effects on this
operator can be safely neglected. The Weinberg operator is generated only at a
very high loop order.

In conclusion, the current constraints on EDMs generate only subleading bounds
compared to the other observables we analyzed.

#### Cosmology, astrophysics, and haloscopes

We mention the impact of haloscope searches, astrophysics, and cosmology
bounds on our model for completeness. Haloscope searches constrain $v_{\Phi}$
through the axion-photon coupling. We take the bounds from AxionLimits ,
rescaling the plot to adapt it to our model correctly. The usual astrophysics
bounds (from AxionLimits ) are completely superseded by the bound due to
$K\rightarrow\pi a$ and hence are irrelevant. For cosmology, the left part of
the blue region in Fig. 5 is taken again from AxionLimits  and is dominated by
the $X$-rays bound, assuming that the ALP constitutes a significant portion of
the dark matter observed in our universe. The rightmost part is obtained via
the requirement that the total ALP decay rate satisfies
$\Gamma_{\text{a,tot}}\geq 3H(T_{\text{BBN}})$, where $T_{\text{BBN}}=4$ MeV,
in order to comply with the constraints from Big Bang Nucleosynthesis
deSalas:2015glj . The rate $\Gamma_{\text{a,tot}}$ has been computed taking
into account all the decay channels of the ALP in our model. For ALP masses
around the QCD confinement scale, this has been approximated with a linear
interpolation between $m_{a}\sim 3m_{\pi}$ and $m_{a}\simeq$ 2 GeV, given the
complications due to flavor and CP violating ALP interactions with hadrons
present in our model. The impact of the heavy $\rho$ and the VLFs is
irrelevant as they quickly decay to SM particles via Eqs. (12) and (14).

#### Direct searches

Low-scale FN models present an exciting opportunity for direct searches at
high-energy colliders. Collider searches for heavy neutral scalars and ALPs
have recently received significant attention Mimasu:2014nea ; Brivio:2017ije ;
Carmona:2022jid ; Esser:2023fdo ; Blasi:2023hvb ; Rygaard:2023dlx ;
Phan:2023dqw ; Li:2024thq . Light ALPs ($m_{a}\lesssim 10$ GeV) either escape
the detector or lead to displaced vertices. Mono$-\gamma$ and mono-jet
searches for light ALPs Mimasu:2014nea  are, however, not competitive with
meson decays. Moving on to EW scale masses, phenomonology in this model is
driven by couplings to second-generation fermions, see Eqs. (17) and (19).
This is because, in the $\mathbb{Z}_{4}$ case, the bottom and tau Yukawas size
are not dynamically explained (unlike in the $\mathbb{Z}_{8}$ model of Section
4). Resonant axion production at hadron colliders has been discussed in Ref.
Bauer:2016rxs . We expect the dominant channels in our case to be
$gc\rightarrow ta$ and $\bar{c}c\rightarrow a$. The cross sections, computed
in Fig. 11 of Ref. Bauer:2016rxs  for $v_{\Phi}=0.5$ TeV, show that one should
expect a large number of events. However, as discussed in Bauer:2016rxs , due
to large backgrounds a dedicated detailed experimental analysis is required.

The flavor-violating ALP-top-charm coupling leads to the $t\rightarrow ca$
decay. A study of this anomalous decay has been carried out in Bauer:2016rxs ,
considering the current and the potential future bounds on
$\text{Br}(t\rightarrow Hc)$ with $H\rightarrow\bar{b}b$ and assuming
$\text{Br}(a\rightarrow\bar{b}b)>80\%$. A major difference between ours and
their scenario is the branching ratio of the decay $a\rightarrow\bar{b}b$.
This would reduce the impact of this search since we expect $a\to c\bar{c}$
decay channel to dominate. A proper assessment would require repeating the
numerical study from Bauer:2016rxs  within our model, which is beyond the
scope of our paper. Therefore, we will not include these bounds in our
figures. Even considering the nominal value they find, the bound from $\Delta
F=2$ processes remains dominant in that region. However, with future
measurements, the situation might change. For light ALPs, the displaced
signatures from top decays could lead to important bounds in the future
Carmona:2022jid .

For completeness, searches for vectorlike quarks and leptons at LHC set lower
bounds on their masses of $\mathcal{O}$(TeV) ATLAS:2024gyc ; CMS:2024bni ,
which cannot compete against bounds from flavor physics in our scenario,
forcing $M\approx v_{\Phi}/\epsilon\gtrsim$ PeV.

To sum up, while the current direct searches provide subleading constraints,
the situation will improve with proposed future high-energy colliders.

Figure 5: A global overview of the phenomenology of a FN ALP in the
$(m_{a},1/v_{\Phi})$ plane. The gray region is excluded by the stability of
the potential, the dotted line shows the $N=4$ expectation, and the dashed
line shows the QCD expectation, as discussed around Figure 1. The rest of the
contours are due to the constraints discussed in detail in Section 3.2. The
red contour is due to the constraints from the UV VLFs, the orange contour is
due to the combined constraints on the FN ALP from flavor observables, while
the green, yellow, and blue contours are due to astrophysics, haloscopes, and
cosmology, respectively. See the surrounding text for more details.

###  3.3 Summary of the results

This subsection summarizes and illustrates the main phenomenological results
derived in the previous subsection. In Fig. 5, we show the exclusion contours
in the $(m_{a},1/v_{\Phi})$ plane for a broad range of ALP masses,
$m_{a}\in[10^{-18},10^{9}]~{}\mathrm{GeV}$. The stability region in gray, the
dotted line for the $N=4$ expectation, and the dashed line for the QCD
expectation are the same as already discussed in Fig. 1. Importantly, we can
immediately notice that following the $N=4$ expectation line, the main
constraint on the $\mathbb{Z}_{4}$ model is due to the vector-like fermions,
illustrating the significance of a UV completion. The contour in red depicts
the combined constraint due to the presence of these states in the UV, which
is dominated by the CP violating kaon mixing observable $\epsilon_{K}$ as
discussed in the previous subsection, setting $v_{\Phi}\gtrsim
2.8~{}\mathrm{TeV}$. The answer to the question posed in the introduction,
“How low can the FN symmetry-breaking scale be?”, aligns with the current
energy frontier.

The ALP flavor phenomenology is less important in this particular scenario,
even though it plays a major role in a large part of the parameter space as
depicted by the orange contour, which combines all of the flavorful ALP
constraints from the previous subsection. Notably, it sets the leading
constraint on the flaxion/axiflavon scenario Calibbi:2016hwq ; Ema:2016ops
through the kinematically allowed $K\to\pi a$ decays. As already mentioned,
the astrophysics constraints depicted in green are less important in scenarios
with flavorful ALPs. In blue, we show the constraints from cosmology mentioned
earlier, with the left part taken from Ref. AxionLimits  and the right coming
from the requirement of not spoiling the success of BBN. Although these
constraints are important in a large part of the parameter space, we now zoom
in to a region that is arguably more interesting for low-scale FN models.

Figure 6: A zoomed-in version of Figure 5, focusing on the parameter space
interesting for low-scale FN models. Again, the gray region is excluded by
potential stability, while the dotted line shows the expectation of $N=4$. We
show the leading constraint from VLFs in red and a compendium of constraints
on the flavorful ALP in the rest of the contours. See the surrounding text for
a detailed discussion.

In Fig. 6, we show a zoomed-in version of Fig. 5, now focusing on the ALP
masses in the range $m_{a}\in[0.1,10^{4}]~{}\mathrm{GeV}$. This region of the
parameter space is of main interest for the FN models studied in this work and
is constrained mainly from flavor physics observables. We emphasize again that
the constraint due to VLFs dominates the natural $\mathbb{Z}_{4}$ model.
However, by allowing for $\lambda_{N}^{\prime}\ll 1$ we can in principle
populate the whole parameter space. Moreover, as we will see in the next
section, the qualitative picture varies little, even for more elaborate
$\mathbb{Z}_{N}$ models. Since the ALP flavor constraints dominate a large
portion of this parameter space, we now split them into multiple contours to
differentiate between observables. In the subsequent discussion, we move from
high to low ALP masses, commenting on the leading constraints for the
different regimes.

In purple, we show the combined constraints from the ALP contributions to
neutral meson mixing observables, here dominated by $\epsilon_{K}$. These
provide the leading constraint for $m_{B}\lesssim m_{a}\lesssim
60~{}\mathrm{GeV}$, setting $v_{\Phi}\gtrsim 10-100~{}\mathrm{TeV}$. Rare
meson decays become kinematically allowed for ALP masses below $m_{B}$, as
depicted in the yellow contour. At first, $B$-meson decays to visible final
states provide the leading constraint, dominated by displaced vertex analyses
of $B\to K^{(\ast)}\mu\mu$ by LHCb, setting $v_{\Phi}\gtrsim
100-1000~{}\mathrm{TeV}$. The region below that is due to $B\to Ka$, setting
an upper and a lower bound on $v_{\Phi}$ as we demand the ALP to behave as an
invisible final state, as discussed in the previous subsection. Finally, for
ALP masses below $m_{K}$, $K\to\pi a$ decays completely dominate the
phenomenology.

For completeness we comment on the subleading LFV phenomenology: in brown
contours, we depict the combined constraints from rare lepton decays,
dominated first by $\tau\to 3\mu$ and then by $\tau\to\mu+\mathrm{invisible}$,
with the flat region being associated to the plateau of $\tau\to
M^{(\ast)}\mu$ where $M^{(\ast)}$ are unflavored mesons. In water green, we
show the constraints from $\ell_{i}\to\ell_{j}\gamma$, with the sensitivity of
up to $v_{\Phi}\gtrsim 10~{}\mathrm{TeV}$. Finally, in pink, we show the
contour from $\mu\to e$ conversion, currently setting bounds up to
$v_{\Phi}\gtrsim$ TeV.

Looking ahead, several experimental advancements will significantly probe the
allowed parameter space, see e.g. Calibbi:2020jvd . First, future bounds on
charged lepton flavor violation from $\mu$ to $e$ conversion Bernstein:2019fyh
; Moritsu:2022lem  will impose stricter constraints on the UV sector, pushing
the FN scale by another order of magnitude. Second, there is a pressing need
for dedicated experimental programs to search for many $M\to M^{\prime}a$ and
$\ell\to\ell^{\prime}a$ processes with displaced ALP vertices
$a\to\ell_{i}\ell_{j},\gamma\gamma$, etc. Finally, both Belle-II and LHCb will
improve the existing limits on rare decays, further pushing into the
unexplored territory.

##  4 A finer model: $\mathbb{Z}_{8}$

The simplicity of the $\mathbb{Z}_{4}$ model comes at the expense of not
perfectly explaining the quark and lepton flavor hierarchies. While the
existence of hierarchies between generations is consistently predicted, this
model fails to account for the smallness of the bottom and tau Yukawas, which
are observed to be $\mathcal{O}(0.01)$ instead of the expected
$\mathcal{O}(1)$. Furthermore, the mass hierarchy between different
generations of down quarks and charged leptons is more compressed compared to
that of the up quarks, despite the $\mathbb{Z}_{4}$ model predicting a
universal $\epsilon$ factor between generations for all three gauge
representations. These additional finer features of the SM flavor puzzle call
for a larger symmetry.

#####  $\mathbb{Z}_{8}$ EFT.

The up-quark sector exhibits a double hierarchy compared to the down sector.
This feature was recently successfully achieved by the
$\mathrm{U}(2)_{q+e^{c}+u^{c}}\times\mathbb{Z}_{2}$ symmetry in Ref.
Antusch:2023shi , which retains the favorable characteristics of
$\mathrm{U}(2)_{q+e^{c}}$, including the neutrino anarchy, while producing a
single (double) hierarchy in the down (up) quarks. The $\mathbb{Z}_{2}$ factor
serves to suppress overall the down sector, which explains the smallness of
$y_{b,\tau}$. This section demonstrates that the smallest discreet symmetry
capable of reproducing the same texture as in Ref. Antusch:2023shi  is
$\mathbb{Z}_{8}$.

To illustrate this, if the hierarchy parameter in the down sector is
$\epsilon$, then in the up sector, it should be $\epsilon^{2}$, leading to
$y_{u}/y_{t}\sim(y_{d}/y_{b})^{2}\sim\epsilon^{4}$. Following the discussion
below Eq. (1), given the largest $n^{f}_{ij}=4$, the minimal symmetry group
required is therefore $\mathbb{Z}_{8}$. The non-trivially charged fields in
this case are $[q_{1}]=-[e_{3}]=-[u_{1}]=2$, $[q_{2}]=-[e_{2}]=-[u_{2}]=1$,
$[d_{i}]=-2$ and $[\ell_{i}]=4$, where $i=1,2,3$. The leading patterns are

$\displaystyle n_{ij}^{d}=\begin{pmatrix}4&4&4\\\ 3&3&3\\\
2&2&2\end{pmatrix},\quad n_{ij}^{e\prime}=\begin{pmatrix}4&3&2\\\ 4&3&2\\\
4&3&2\end{pmatrix},\quad n_{ij}^{u}=\begin{pmatrix}4&3&2\\\ 3&2&1\\\
2&1&0\end{pmatrix},$ (59)

while $n_{ij}^{d\prime}$, $n_{ij}^{e}$ and $n_{ij}^{u\prime}$ give comparable
or subleading contributions to $Y^{f}_{ij}$.

A few comments are in order. First, we can see explicitly why the
$\mathbb{Z}_{8}$ symmetry is necessary. A replacement
$\Phi^{M}\rightarrow(\Phi^{*})^{N-M}$ would disrupt the hierarchies if $N<8$.
In addition, the universal charges for $d_{i}$ and $\ell_{i}$ are chosen to be
different, even though they lead to the same singular values
$\hat{y}^{d,e}_{ii}\sim\\{\epsilon^{4},\epsilon^{3},\epsilon^{2}\\}$. This
choice allows the CKM matrix to be approximately a unit matrix while implying
no selection rules on the neutrino sector. In other words, the Weinberg
operator

$-\mathcal{L}\supset\frac{y^{\nu}_{ij}}{M_{\nu}}\ell_{i}\ell_{j}HH\,,$ (60)

is neutral under the $\mathbb{Z}_{8}$ symmetry. Therefore, it features
anarchic coefficients, which is consistent with the large observed mixings in
the PMNS matrix.232323In the type-I seesaw scenario, completing this operator
with heavy right-handed neutrinos $\nu_{i}$ means $[\nu_{i}]=4$, ensuring that
both the Yukawa interaction and Majorana mass terms are allowed by the
$\mathbb{Z}_{8}$ symmetry.

The details of fitting $x^{f}_{ij}$ and $x^{f\prime}_{ij}$ from Eq. (1) to the
observed SM flavor parameters are presented in Appendix A.2. As demonstrated
by the numerical benchmark point, this FN EFT offers an excellent description
of the SM flavor parameters for $\epsilon\simeq 6.6\times 10^{-2}$ with all
fit parameters being $\mathcal{O}(1)$, as anticipated.

##### UV completion.

The set of VLF with masses at the scale $M$, completing Eq. (1) for the
$\mathbb{Z}_{8}$ charges listed above, is somewhat more complex than in the
previous section. Let us start with the quark sector. The charges of the new
VLF fields are

$\begin{split}[Q^{a}_{2}]&=0,\quad[Q_{1}]=1,\\\
[U^{a}_{2}]&=0,\quad[U_{1}]=-1,\\\
[D^{i}_{2}]&=0,\quad[D^{i}_{1}]=-1.\end{split}$ (61)

The gauge representations correspond to the SM $q$, $u$, and $d$ fields,
respectively, with $a=1,2$ and $i=1,2,3$. Adding these states in the UV allows
for the tree-level completion of the effective quark Yukawa operators. The
corresponding diagrams are analogous to those in Fig. 4, featuring
$n^{f}_{ij}$ VLF propagators. The effective Yukawa matrix entries
parametrically scale as

$Y^{f}_{ij}\sim(y_{*})^{n^{f}_{ij}+1}\epsilon^{n^{f}_{ij}}.$ (62)

Here, $y_{*}$ generically represents a dimension-4 coupling in the UV
Lagrangian. As previously discussed, the product of several $\mathcal{O}(0.3)$
couplings can contribute to hierarchies. Consequently, we expect the variance
in the UV couplings to be even smaller than in the EFT fit presented above,
similar to what was observed for the $\mathbb{Z}_{4}$ model.

Let us now turn our discussion to the charged lepton sector. The required set
of VLF, all in the gauge representation of the SM field $\ell$, is

$[L_{N}]=0,\quad[L^{a}_{N-1}]=-1,\quad[L^{i}_{N-2}]=-2,\quad[L^{i}_{N-3}]=-3,$
(63)

where $N=8$. Integrating them out at the tree level produces the desired
charged lepton effective Yukawa operators. An example diagram involves
connecting $e_{1}$ (which has charge zero) with $\ell_{i}$ (charge 4) through
the aforementioned chain of VLF, attaching $H$, and subsequently adding four
$\Phi^{*}$.

The final model-building task is to endow the ALP with a non-null mass. The
simplest option is to consider the wheel diagram already introduced in Fig. 2,
which requires a full set of VLF representation charges. Since we already have
half of the $L_{k}$ charges, we only need, in addition,242424Due to the large
number of $\mathrm{SU}(2)_{\rm L}$ doublets in this model, the gauge coupling
$g_{2}$ develops a Landau pole at approximately $10^{14}$ GeV for $M=1$ PeV
(1-loop estimate). In contrast, for the $\mathbb{Z}_{4}$ model, the couplings
remain perturbative up to the Planck scale.

$[L_{1}]=1,\quad[L_{2}]=2,\quad[L_{3}]=3,\quad[L_{4}]=4,$ (64)

which together with Eq. (63) gives the diagram in Fig. 2. The prediction for
the ALP mass aligns precisely with the one discussed in Section 2 for the
$N=8$ case, as shown in Fig. 1. We expect
$m_{a}/v_{\Phi}\sim\\{10^{-6},10^{-4}\\}$, see Eq. (9) and the discussion
below it. Thus, unlike in the $\mathbb{Z}_{4}$ model, here we expect a light
ALP, potentially leading to interesting effects described by the SM$+a$ EFT.

##### Phenomenology.

First, let us discuss the UV contributions to operators without $\Phi$
insertions before moving on to the ALP.

_VLF_. The leading effect arises from $\mathcal{L}\supset
y^{ij}_{D_{1}}\bar{D}^{i}_{1}\Phi d^{j}$, which, after integrating out
$D^{i}_{1}$ at the one-loop level, matches to the
$(\bar{d}_{1}\gamma_{\mu}d_{2})^{2}$ operator. This is analogous to the
$\mathbb{Z}_{4}$ model, where the $H$ and $Q^{1}_{2}$ fields appeared in the
box diagram, whereas in this case, we have the $\Phi$ and $D_{1}^{i}$ fields
(see Eq. (38) and the subsequent discussion). Therefore, the leading
constraint comes from $\epsilon_{K}$ and can be easily derived by rescaling
the $\mathbb{Z}_{4}$ case by the $\epsilon$ ratios in the two models. Since
$\epsilon$ is larger in the $\mathbb{Z}_{8}$ model, the bound becomes
$v_{\Phi}\gtrsim 40$ TeV, up to $\mathcal{O}(1)$ variation in
$y^{ij}_{D_{1}}$.

Figure 7: Constraints in the $(m_{a},1/v_{\Phi})$ plane, similar to Figure 6,
but now for the $\mathbb{Z}_{8}$ model. The gray region is excluded by vacuum
stability, while the dashed black lines show the expected range in a $N=8$
wheel model. The red and purple (full) contours show the major differences
with respect to the $\mathbb{Z}_{4}$ case (dashed). The rest of the
constraints are qualitatively the same as in the $\mathbb{Z}_{4}$ case. See
the surrounding text for more details.

_ALP_. Assuming $v_{\Phi}$ saturates this bound, we expect $m_{a}\in(0.04,4)$
GeV. This mass range coincides with numerous significant flavor constraints in
the $\mathbb{Z}_{4}$ model, providing leading bounds on the FN scale. We first
need to examine the ALP coupling matrices to compare the predictions between
the two models. In this case, the hierarchies in the left-handed rotation
matrices $U_{L}^{d}$ and $U_{L}^{u}$ are similar. Hence, the CKM receives, in
principle, a sizeable contribution from both. Furthermore, this case has many
more free parameters for $Y_{u}$. For this reason, we find it more instructive
to report the general structure of the couplings in the $\mathbb{Z}_{8}$ model
up to $\mathcal{O}(1)$ factors:

$\displaystyle\begin{gathered}c^{d}\sim\begin{pmatrix}m_{d}&m_{s}\frac{m_{d}}{m_{s}}&m_{b}\frac{m_{d}}{m_{b}}\\\
m_{d}\frac{m_{d}}{m_{s}}&m_{s}&m_{b}\frac{m_{s}}{m_{b}}\\\
m_{d}\frac{m_{d}}{m_{b}}&m_{s}\frac{m_{s}}{m_{b}}&m_{b}\\\
\end{pmatrix},\qquad
c^{u}\sim\begin{pmatrix}m_{u}&m_{c}\sqrt{\frac{m_{u}}{m_{c}}}&m_{t}\sqrt{\frac{m_{u}}{m_{t}}}\\\
m_{c}\sqrt{\frac{m_{u}}{m_{c}}}&m_{c}&m_{t}\sqrt{\frac{m_{c}}{m_{t}}}\\\
m_{t}\sqrt{\frac{m_{u}}{m_{t}}}&m_{t}\sqrt{\frac{m_{c}}{m_{t}}}&m_{t}\frac{m_{c}}{m_{t}}\\\
\end{pmatrix},\\\
c^{e}\sim\begin{pmatrix}m_{e}&m_{e}\frac{m_{e}}{m_{\mu}}&m_{e}\frac{m_{e}}{m_{\tau}}\\\
m_{\mu}\frac{m_{e}}{m_{\mu}}&m_{\mu}&m_{\mu}\frac{m_{\mu}}{m_{\tau}}\\\
m_{\tau}\frac{m_{e}}{m_{\tau}}&m_{\tau}\frac{m_{\mu}}{m_{\tau}}&m_{\tau}\\\
\end{pmatrix}.\end{gathered}$ (67)

The most prominent differences compared to the $\mathbb{Z}_{4}$ case are
$c^{e}_{33},c^{d}_{33}$, which now are unsuppressed due to the fact that the
third generation is charged, and of course, the whole $c^{u}$. The fact that
$c^{u}$ mixings scale as the square root of the ratios of masses is due to the
fact that neither $q_{i}$ nor $u_{i}$ are “universally” charged, unlike
$d_{i}$ and $e_{i}$. In comparison with the latter two, the up sector has
doubled suppression. Hence, the mixings are “split” between the left and right
rotation matrices, which give comparable contributions to $c^{u}$, see Eq.
(108). Because of this, ALP contributions to $D-\bar{D}$ mixing become now
more important Carmona:2021seb . In particular, if the CKM matrix is still
fully given by $U_{L}^{d}$, as is for simplicity assumed in the example of
App. A.2, the constraint from $\mathbb{Z}_{4}$ due to $\epsilon_{K}$ remains
essentially the same. Hence, in this case, CP violation in D meson mixing is
actually the most important $\Delta F=2$ probe. In the general case in which
the CKM receives sizeable contributions from both the diagonalizing matrices,
we find that the $\epsilon_{K}$ bound becomes slightly stronger than in the
$\mathbb{Z}_{4}$ case, leading to bounds comparable to that from CP violation
in $D-\bar{D}$ mixing.

We illustrate the phenomenology of the $\mathbb{Z}_{8}$ model in Figure 7. The
dashed black lines now show the expectation in the $(m_{a},1/v_{\Phi}$) plane
for a $N=8$ wheel model, as already discussed around Figure 1. With dashed red
and purple lines we show the $\Delta F=2$ and VLF constraints in the case of
$\mathbb{Z}_{4}$, as in Figure 6. The full contours show the respective bounds
in the $\mathbb{Z}_{8}$ case, both stronger in line with the discussion above.
It is interesting to note how in this scenario, contrary to the case of
$\mathbb{Z}_{4}$, the leading phenomenology is due to the light flavorful ALP
and no longer due to the VLFs. In particular, the $\Delta F=2$ constraints
together with displaced vertex searches in $B\to K^{(*)}\mu\mu$
(quantitatively the same as in $\mathbb{Z}_{4}$), set the most stringent
constraints of order $v_{\Phi}\gtrsim 10^{5}-10^{7}~{}\mathrm{GeV}$ for the
corresponding expected mass range of $m_{a}\in[1,20]~{}\mathrm{GeV}$.

##  5 Conclusion

The presence of an ALP in quark and lepton flavor physics opens up a diverse
array of new phenomena, making it an exciting area of research for
phenomenologists. The appropriate theoretical framework to study these effects
is the SM+$a$ EFT, which constructs a tower of local operators that respect
both the SM symmetry and the approximate ALP shift symmetry. However,
depending on the coefficients of the operators and the ALP mass, which are
input parameters a priori, the expected patterns of deviations vary widely,
thereby reducing the predictive power of this approach.

In the bulk of this paper, we surveyed the flavor physics phenomenology in the
presence of an ALP, adding, however, a deeper UV context to the study.
Specifically, we imposed a discrete FN symmetry motivated by the SM flavor
puzzle. This symmetry provides power counting and selection rules that
correlate various observables, offering a structured and focused framework for
studying ALP effects in flavor physics. The main results are summarized in
Figures 6 and 7, which depict the parameter space of the ALP mass versus the
FN scale. The two plots interpret a plethora of flavor physics measurements
discussed in Section 3.2, organized into suitable categories, within two
concrete models based on $\mathbb{Z}_{4}$ and $\mathbb{Z}_{8}$ FN symmetries,
respectively. The two choices allow us to quantify the uncertainty due to
variations in the realization of the FN mechanism.

Prior to the comprehensive phenomenological work, we introduced a whole class
of FN models based on $\mathbb{Z}_{N}$ symmetries. While a commonly considered
perturbatively-exact $\mathrm{U}(1)_{\text{FN}}$ predicts a very light
(pseudo) Goldstone boson, the advantage of the discrete symmetries is that
$m_{a}$ can be considerably heavier, opening up new phenomenological
possibilities as discussed above. We began in Section 2 by charting the ALP
parameter space based on theoretical considerations, as summarized in Fig. 1,
and by presenting concrete predictive models, such as the wheel model shown in
Fig. 2. This setup is more predictive than an explicitly broken
$\mathrm{U}(1)_{\text{FN}}$, which can produce correct textures provided
$m_{a}\lesssim v_{\Phi}$ but does not correlate the two parameters.
Nevertheless, our results qualitatively apply to the latter case as well.
Assuming marginal couplings of $\mathcal{O}(1)$, our models predict a region
in $(m_{a},v_{\Phi})$ as illustrated in the plots for the two chosen examples.
Notably, the $\mathbb{Z}_{4}$ model is unique in predicting $m_{a}\sim
v_{\Phi}$ due to the presence of a $\mathrm{U}(1)_{\text{FN}}$-breaking
dimension-4 term in the potential.

Another advantage of the $\mathbb{Z}_{N}$ class is the ability to balance
between simplicity versus “baroqueness”, controlled by the value of $N$.
Unsurprisingly, a large enough $N$ allows for a perfect fit of the SM flavor
parameters. However, we prioritized simplicity, accepting an increased
variance in the UV parameters. Following this path, we uncovered a remarkably
simple and elegant model based on the $\mathbb{Z}_{4}$ group (Section 3.1). By
charging only left-handed quarks and right-handed charged leptons, this model
generates hierarchies in the charged fermion masses and left-handed quark
mixing matrices while maintaining an anarchic neutrino flavor structure. This
minimal symmetry framework produces similar Yukawa textures as the full
$\mathrm{U}(2)_{q+e}$ symmetry recently proposed in Antusch:2023shi . To
address subtleties of the SM flavor puzzle, such as the smallness of
$y_{b,\tau}$ and the difference in up versus down quark mass hierarchies, we
introduced a refined $\mathbb{Z}_{8}$ model in Section 4. This model is
carefully crafted, drawing inspiration from the novel
$\mathrm{U}(2)_{q+e^{c}+u^{c}}\times\mathbb{Z}_{2}$ symmetry proposed recently
in Antusch:2023shi .

A crucial aspect of our work is that in both scenarios, we constructed
explicit renormalizable completions featuring new vector-like fermions with
masses above the FN scale, given by $v_{\Phi}/\epsilon$ where $\epsilon$ is
the FN spurion controlling the Yukawa textures. These new fields generate
effective Yukawa operators in Eq. (1), but can also lead to deviations in rare
flavor-changing neutral currents. An intriguing question arises about how the
constraints on the minimal UV sector, which is essential for completing the FN
mechanism, compare to those from a light ALP. The explicit models provide a
basis for conducting precision calculations and studying this interplay
between the ALP and the irreducible UV-induced effects. The so-called VLF
bounds, shown in red in Figures 6 and 7, are dominant for a heavy ALP
($m_{a}\gtrsim 100$ GeV). This parameter space region permits the lowest FN
scale $v_{\Phi}$. For instance, in the case of $\mathbb{Z}_{4}$ predicting
$m_{a}\sim v_{\phi}$, the FN symmetry could be restored even at the TeV scale.

In conclusion, with its rich phenomenology and solid theoretical foundations,
the FN ALP represents an exciting physics frontier for ongoing and future
experiments in flavor and collider physics.

###### Acknowledgements.

This work received funding from the Swiss National Science Foundation (SNF)
through the Eccellenza Professorial Fellowship, project number 186866, “Flavor
Physics at the High Energy Frontier.”

##  Appendix A UV and EFT parameters matching

###  A.1 $\mathbb{Z}_{4}$ model

####  A.1.1 EFT

In the $\mathbb{Z}_{4}$ model, other than $\Phi$, only $q_{i}$ and $e_{i}$ are
charged under the flavor symmetry, with charges $[q_{1}]=-[e_{1}]=2$,
$[q_{2}]=-[e_{2}]=1$. Hence, one can always perform a set of unitary rotations
of $u_{i},d_{i},\ell_{i}$ such that the SM Yukawas read, after $\Phi$
condensation,

$\displaystyle
Y_{d}=\begin{pmatrix}z_{d_{1}}\epsilon^{2}&z_{d_{2}}\epsilon^{2}&z_{d_{3}}\epsilon^{2}\\\
0&y_{d_{2}}\epsilon&y_{d_{3}}\epsilon\\\ 0&0&x_{d_{3}}\end{pmatrix},\quad
Y_{u}=\begin{pmatrix}z_{u_{1}}\epsilon^{2}&z_{u_{2}}\epsilon^{2}&z_{u_{3}}\epsilon^{2}\\\
0&y_{u_{2}}\epsilon&y_{u_{3}}\epsilon\\\ 0&0&x_{u_{3}}\end{pmatrix},\quad
Y_{e}=\begin{pmatrix}z_{e_{1}}\epsilon^{2}&0&0\\\
z_{e_{2}}\epsilon^{2}&y_{e_{2}}\epsilon&0\\\
z_{e_{3}}\epsilon^{2}&y_{e_{3}}\epsilon&x_{e_{3}}\end{pmatrix},$ (68)

with $\epsilon\equiv v_{\Phi}/\sqrt{2}M$. The parameters in Eq. (68)
correspond to $x_{ij}^{f}$ in Eq. (1) and exploiting the residual rephasing
symmetry of the SM fields can all be made real, except five of them which we
take to be $z_{d_{3}},z_{u_{2}},z_{u_{3}},y_{u_{3}}$ and $z_{e_{3}}$. A
perturbative singular value decomposition of Eq. (68) reveals that the SM
Yukawa matrices in the mass basis are simply given by the diagonal elements,

$\displaystyle\begin{aligned}
(y_{d},y_{s},y_{b})&\simeq(z_{d_{1}}\epsilon^{2},y_{d_{2}}\epsilon,x_{d_{3}}),\\\
(y_{u},y_{c},y_{t})&\simeq(z_{u_{1}}\epsilon^{2},y_{u_{2}}\epsilon,x_{u_{3}}),\\\
(y_{e},y_{\mu},y_{\tau})&\simeq(z_{e_{1}}\epsilon^{2},y_{e_{2}}\epsilon,x_{e_{3}}),\end{aligned}$
(69)

with the leading rotation matrices

$\displaystyle\begin{gathered}U^{d}_{L}\simeq\begin{pmatrix}1&\frac{z_{d_{2}}}{y_{d_{2}}}\epsilon&\frac{z_{d_{3}}}{x_{d_{3}}}\epsilon^{2}\\\
-\frac{z_{d_{2}}}{y_{d_{2}}}\epsilon&1&\frac{y_{d_{3}}}{x_{d_{3}}}\epsilon\\\
\frac{y_{d_{3}}z_{d_{2}}-y_{d_{2}}z_{d_{3}}^{*}}{x_{d_{3}}y_{d_{2}}}\epsilon^{2}&-\frac{y_{d_{3}}}{x_{d_{3}}}\epsilon&1\end{pmatrix}\simeq\begin{pmatrix}1&\frac{m_{d}}{m_{s}}\frac{z_{d_{2}}}{z_{d_{1}}}&\frac{m_{d}}{m_{b}}\frac{z_{d_{3}}}{z_{d_{1}}}\\\
-\frac{m_{d}}{m_{s}}\frac{z_{d_{2}}}{z_{d_{1}}}&1&\frac{m_{s}}{m_{b}}\frac{y_{d_{3}}}{y_{d_{2}}}\\\
\frac{m_{d}}{m_{b}}\frac{y_{d_{3}}z_{d_{2}}-y_{d_{2}}z_{d_{3}}^{*}}{y_{d_{2}}z_{d_{1}}}&-\frac{m_{s}}{m_{b}}\frac{y_{d_{3}}}{y_{d_{2}}}&1\end{pmatrix},\end{gathered}$
(71)

and similarly for $U^{u}_{L}$ and $(U^{e}_{R})^{*}$. The Eq. (71) shows that
for $\mathcal{O}(1)$ parameters, the CKM mixings can be reproduced quite well,
and clarifies why in Sec. 3 we assumed $V_{\text{CKM}}\simeq U^{d}_{L}$: the
bigger hierarchies in the up sector suppress the off-diagonal entries of
$U^{u}_{L}$ to the point that a significant fine-tuning would be required for
it to actually contribute to $V_{\text{CKM}}$.

Eqs. (69) and (71) can be used to fix some of the parameters in Eq. (68),
including $\epsilon$. The optimal benchmark, where the parameters are close to
$\mathcal{O}(1)$ and which is used to derive the bounds in the main text, is
achieved for $\epsilon\approx 4.4\times 10^{-3}$ with

$\displaystyle\begin{gathered}Y_{d}\simeq\begin{pmatrix}0.55\epsilon^{2}&2.5\epsilon^{2}&(0.73-1.8i)\epsilon^{2}\\\
0&0.049\epsilon&0.10\epsilon\\\ 0&0&0.011\end{pmatrix},\qquad
Y_{u}\simeq\begin{pmatrix}0.25\epsilon^{2}&z_{u_{2}}\epsilon^{2}&z_{u_{3}}\epsilon^{2}\\\
0&0.57\epsilon&y_{u_{3}}\epsilon\\\ 0&0&0.71\end{pmatrix},\\\
Y_{e}\simeq\begin{pmatrix}0.15\epsilon^{2}&0&0\\\
z_{e_{2}}\epsilon&0.14\epsilon&0\\\
z_{e_{3}}&y_{e_{3}}&0.01\end{pmatrix}.\end{gathered}$ (74)

where for the fit, we used the value of the SM Yukawas at 100 TeV from
Huang:2020hdv . The presence of some $\mathcal{O}(10^{-1})$ hierarchies in the
lighter generations is not a problem, as it can naturally be explained by a
product of several UV couplings of order $\mathcal{O}(0.3)$. Conversely, the
small values of the bottom and tau Yukawas, approximately
$\mathcal{O}(10^{-2})$, as indicated by $x_{d_{3}}$ and $x_{e_{3}}$, cannot be
accounted for within this model and must be accepted as given.

####  A.1.2 UV

Integrating out the VLF in Eq. (12) at tree-level leads to a FN EFT with
Yukawas

$Y^{(u,d)}=\begin{pmatrix}\epsilon_{Q}^{2}r_{Q}x_{1}^{q}x_{12}^{qa}\bm{y}^{(u,d)a}\\\
\epsilon_{Q}x_{2}^{qa}\bm{y}^{(u,d)a}\\\ \bm{z}^{(u,d)}\end{pmatrix}\,,\qquad
Y^{e}=\begin{pmatrix}\epsilon_{E}^{2}r_{E}x_{1}^{e}x_{12}^{ea}\bm{y}^{ea}~{}~{}\epsilon_{E}x_{2}^{ea}\bm{y}^{ea}~{}~{}\bm{z}^{e\,T}\end{pmatrix},$
(75)

where $\epsilon_{Q}\equiv v_{\Phi}/\sqrt{2}M_{Q_{2}},\epsilon_{E}\equiv
v_{\Phi}/\sqrt{2}M_{E_{2}}$. For simplicity, we have assumed degenerate
$Q_{2}^{a}$ and $E_{2}^{a}$ masses and $r_{Q}=M_{Q_{2}}/M_{Q_{1}}$,
$r_{E}=M_{E_{2}}/M_{E_{1}}$. Boldface quantities denote vectors in flavor
space. To match Eq. (75) to (68) and ultimately the SM masses and mixings, it
is convenient to pick a basis in the UV theory, which simplifies calculation
as much as possible. Rotation of the fields in Eq. (12) allows us to put the
couplings in the following schematic form

$\displaystyle\begin{gathered}M_{Q_{1}}\sim\mathds{R},\quad
M_{Q_{2}}\sim\text{diag}(\mathds{R},\mathds{R}),\quad
M_{E_{1}}\sim\mathds{R},\quad
M_{E_{2}}\sim\text{diag}(\mathds{R},\mathds{R}),\\\
{\bm{z}}_{d}\sim(0,0,\mathds{R}),\quad{\bm{z}}_{u}\sim(0,0,\mathds{R}),\quad{\bm{z}}_{e}\sim(0,0,\mathds{R}),\\\
{\bm{y}}^{da}\sim\begin{pmatrix}\mathds{R}&\mathds{C}&\mathds{C}\\\
0&\mathds{R}&\mathds{R}\end{pmatrix},\quad{\bm{y}}^{ua}\sim\begin{pmatrix}\mathds{R}&\mathds{C}&\mathds{C}\\\
0&\mathds{R}&\mathds{C}\end{pmatrix},\quad{\bm{y}}^{ea}\sim\begin{pmatrix}0&\mathds{R}\\\
\mathds{R}&\mathds{C}\\\ \mathds{R}&\mathds{C}\end{pmatrix},\\\
x_{2}^{qa}\sim(0,\mathds{R}),\quad
x_{12}^{qa}\sim(\mathds{R},\mathds{C}),\quad x_{1}^{q}\sim\mathds{R},\\\
x_{2}^{ea}\sim(\mathds{R},0),\quad
x_{12}^{ea}\sim(\mathds{R},\mathds{R}),\quad
x_{1}^{e}\sim\mathds{R},\end{gathered}$ (81)

which leads to Yukawas with the triangular structure of (68). Unfortunately,
the matching is still not straightforward as not all the phases in (81) are
physical in the EFT. A spurious phase in
$\delta^{d}=\,\text{arg}\,(Y_{d})_{12}$ must be removed with the rotation
$d_{2,3}\rightarrow e^{-i\delta^{d}}d_{2,3},\bar{q}_{2,3}\rightarrow
e^{i\delta^{d}}\bar{q}_{2,3},u_{2,3}\rightarrow e^{-i\delta^{d}}u_{2,3}$ and
must be taken into account when matching the UV parameters to (69), (71). Note
that this rotation also affects the off-diagonal entries of $Y_{u}$ (unfixed
in the EFT) and the Wilson coefficient of higher-dimensional operators in a
non-trivial manner. Similarly, a spurious phase $\delta^{e}$ in $(Y_{e})_{21}$
also appears, which can be removed with
$\bar{\ell}_{2,3}\rightarrow\bar{\ell}_{2,3}e^{-i\delta^{e}},e_{2,3}\rightarrow
e_{2,3}e^{i\delta^{e}}$.

The matching conditions for the down and quarks and the leptons read:

$\displaystyle\begin{cases}r_{Q}x_{1}^{q}x_{12}^{q1}y^{d1}_{1}=z_{d_{1}}=\hat{y}_{d}/\epsilon^{2}\\\
x_{2}^{q2}y_{2}^{d2}=y_{d_{2}}=\hat{y}_{s}/\epsilon\\\
x_{2}^{q2}y_{3}^{d2}=y_{d_{3}}=A\lambda_{c}^{2}\hat{y}_{b}/\epsilon\\\
r_{Q}x_{1}^{q}[x_{12}^{q1}y^{d1}_{2}+x_{12}^{q2}y_{2}^{d2}]=\rho_{d}e^{i\delta_{d}},\qquad\rho_{d}=z_{d_{2}}=\lambda_{c}y_{d_{2}}/\epsilon=\lambda_{c}\hat{y}_{s}/\epsilon^{2}\\\
r_{Q}x_{1}^{q}[x_{12}^{q1}y^{d1}_{3}+x_{12}^{q2}y_{3}^{d2}]e^{-i\delta_{d}}=z_{d_{3}}=A\lambda_{c}^{3}(\rho-i\eta)\hat{y}_{b}/\epsilon^{2}\end{cases}\qquad$
$\displaystyle(\text{down quarks}),$ (82)
$\displaystyle\begin{cases}r_{Q}x_{1}^{q}x_{12}^{q1}y^{u1}_{1}=z_{u_{1}}=\hat{y}_{u}/\epsilon^{2}\\\
x_{2}^{q2}y_{2}^{u2}=y_{u_{2}}=\hat{y}_{c}/\epsilon\\\
x_{2}^{q2}y_{3}^{u2}=y_{u_{3}}\\\
r_{Q}x_{1}^{q}[x_{12}^{q1}y^{u1}_{2}+x_{12}^{q2}y_{2}^{u2}]e^{-i\delta^{d}}=z_{u_{2}}\\\
r_{Q}x_{1}^{q}[x_{12}^{q1}y^{u1}_{3}+x_{12}^{q2}y_{3}^{u2}]e^{-i\delta^{d}}=z_{u_{3}}\end{cases}\qquad$
$\displaystyle(\text{up quarks}),$ (83)
$\displaystyle\begin{cases}r_{E}x_{1}^{e}x_{12}^{e2}y^{e2}_{1}=z_{e_{1}}=\hat{y}_{e}/\epsilon^{2}\\\
x_{2}^{e1}y_{2}^{e1}=y_{e_{2}}=\hat{y}_{\mu}/\epsilon\\\
x_{2}^{e1}y_{3}^{e1}=y_{e_{3}}\\\
r_{E}x_{1}^{e}[x_{12}^{e2}y^{e2}_{2}+x_{12}^{e1}y_{2}^{e1}]=\rho_{e}e^{i\delta^{e}},\qquad\rho_{e}=z_{e_{2}}\\\
r_{E}x_{1}^{e}[x_{12}^{e2}y^{e2}_{3}+x_{12}^{e1}y_{3}^{e1}]e^{-i\delta^{e}}=z_{e_{3}}\end{cases}\qquad$
$\displaystyle(\text{leptons}).$ (84)

The solution to this system of equations is clearly not unique. Unfortunately,
we were not able to completely solve this system analytically due to the last
two equations in (82), which involve non-polynomial functions of the UV
parameters and can be expressed as

$\displaystyle\begin{cases}|x_{2}+z|=b\\\ |x_{3}+az|=d\\\
\,\text{arg}\,[x_{3}+az]-\,\text{arg}\,[x_{2}+z]=c\end{cases}$ (85)

where
$a=A\lambda_{c}^{2}\hat{y}_{b}/\hat{y}_{s},b=\lambda_{c}\hat{y}_{s}/\hat{y}_{d},c=\,\text{arg}\,A\lambda_{c}^{3}(\rho-i\eta)$
and $d=|A\lambda_{c}^{3}(\rho-i\eta)|\hat{y}_{b}/\hat{y}_{d}$, and we
introduced the combinations
$x_{2}=y^{d1}_{2}/y^{d1}_{1},x_{3}=y^{d1}_{3}/y^{d1}_{1},z=x_{12}^{q2}y^{d2}_{2}/(x_{12}^{q1}y^{d1}_{1})$.
To derive the bounds in Section 3, we pick a representative benchmark obtained
by numerically solving (85) (again, solutions are not unique) and choosing the
remaining parameters in such a way to reproduce (74). This reads for the
quarks

$\displaystyle\begin{gathered}{\bm{z}}_{d}=(0,0,y_{b}),\quad{\bm{z}}_{d}=(0,0,y_{t}),\\\
{\bm{y}}^{da}=\begin{pmatrix}0.26&-0.28-0.96i&-0.95+0.24i\\\
0&0.10&0.21\end{pmatrix},\quad{\bm{y}}^{ua}=\begin{pmatrix}0.12&-0.57-0.20i&0.24+0.16i\\\
0&1.2&-0.86+0.14i\end{pmatrix},\\\ r_{Q}=4,\quad x_{1}^{q}=0.96,\quad
x_{2}^{qa}=\left(0,0.47\right),\quad
x_{12}^{qa}=\left(0.55,0.07-0.99i\right),\end{gathered}$ (89)

and for the leptons

$\displaystyle\begin{gathered}{\bm{z}}_{e}=(0,0,y_{\tau}),\\\
{\bm{y}}^{ea}=\begin{pmatrix}0&0.69\\\ 0.38&0.47+0.23i\\\
0.49&0.19+0.52i\end{pmatrix},\\\ r_{E}=0.5,\quad x_{2}^{ea}=(0.39,0),\quad
x_{12}^{ea}=(0.45,0.54),\quad x_{1}^{e}=0.8.\end{gathered}$ (93)

with $\delta^{d}=-1.2$ and $\delta^{e}=0.28$. This benchmark fixes the
remaining entries of (75), after the $\delta^{d},\delta^{u}$ rotation, to
$z_{e_{2}}=0.18,z_{e_{3}}=0.15+0.07i,y_{e_{3}}=0.19$ and
$z_{u_{2}}=2.3+0.06i,z_{u_{3}}=-2.1+0.56i,y_{u_{3}}=-0.42+0.07i$.

###  A.2 $\mathbb{Z}_{8}$ model

In this case, with flavor rotations, we can put the Yukawa matrices in the
form

$\displaystyle
Y_{d}=\begin{pmatrix}z_{d_{1}}\epsilon^{4}&z_{d_{2}}\epsilon^{4}&z_{d_{3}}\epsilon^{4}\\\
0&y_{d_{2}}\epsilon^{3}&y_{d_{3}}\epsilon^{3}\\\
0&0&x_{d_{3}}\epsilon^{2}\end{pmatrix},\quad
Y_{u}=\begin{pmatrix}z_{u_{1}}\epsilon^{4}&z_{u_{2}}\epsilon^{3}&z_{u_{3}}\epsilon^{2}\\\
y_{u_{1}}\epsilon^{3}&y_{u_{2}}\epsilon^{2}&y_{u_{3}}\epsilon\\\
x_{u_{1}}\epsilon^{2}&x_{u_{2}}\epsilon&x_{u_{3}}\end{pmatrix},\quad
Y_{e}=\begin{pmatrix}z_{e_{1}}\epsilon^{4}&0&0\\\
z_{e_{2}}\epsilon^{4}&y_{e_{2}}\epsilon^{3}&0\\\
z_{e_{3}}\epsilon^{4}&y_{e_{3}}\epsilon^{3}&x_{e_{3}}\epsilon^{2}\end{pmatrix},$
(94)

where in $Y_{d}$ and $Y_{e}$ the only complex entries are $z_{d_{3}}$ and
$z_{e_{3}}$, while in $Y_{u}$ all entries are complex except the diagonal
ones. Assuming for simplicity that the off-diagonal entries of $Y_{u}$ are
small, the SM Yukawas in the mass basis are given as

$\displaystyle\begin{aligned}
(y_{d},y_{s},y_{b})&\simeq(z_{d_{1}}\epsilon^{4},y_{d_{2}}\epsilon^{3},x_{d_{3}}\epsilon^{2}),\\\
(y_{u},y_{c},y_{t})&\simeq(z_{u_{1}}\epsilon^{4},y_{u_{2}}\epsilon^{2},x_{u_{3}}),\\\
(y_{e},y_{\mu},y_{\tau})&\simeq(z_{e_{1}}\epsilon^{4},y_{e_{2}}\epsilon^{3},x_{e_{3}}\epsilon^{2}),\end{aligned}$
(95)

with the CKM matrix still given by Eq. (71). We find that the best benchmark
in this case is given by $\epsilon\simeq 6.6\times 10^{-2}$, leading to

$\displaystyle\begin{gathered}Y_{d}\simeq\begin{pmatrix}0.55\epsilon^{4}&2.5\epsilon^{4}&(0.73-1.84i)\epsilon^{4}\\\
0&0.74\epsilon^{3}&1.51\epsilon^{3}\\\ 0&0&2.4\epsilon^{2}\end{pmatrix},\quad
Y_{u}\simeq\begin{pmatrix}0.25\epsilon^{4}&z_{u_{2}}\epsilon^{3}&z_{u_{3}}\epsilon^{2}\\\
y_{u_{1}}\epsilon^{3}&0.57\epsilon^{2}&y_{u_{3}}\epsilon\\\
x_{u_{1}}\epsilon^{2}&x_{u_{2}}\epsilon&0.71\end{pmatrix},\\\
Y_{e}\simeq\begin{pmatrix}0.15\epsilon^{4}&0&0\\\
z_{e_{2}}\epsilon^{3}&2.1\epsilon^{3}&0\\\
z_{e_{3}}\epsilon^{2}&y_{e_{3}}\epsilon^{2}&2.3\epsilon^{2}\end{pmatrix}.\end{gathered}$
(98)

Compared to Eq. (68) the fit in this case is better in terms of the parameters
being close to $\mathcal{O}(1)$, with a maximum spread between the smallest
and the biggest ones of $\mathcal{O}(10)$.

##  Appendix B Derivation of the ALP couplings

Assuming a generic Yukawa Lagrangian

$\mathcal{L}\supset-
y_{ij}^{f}\left(\frac{\Phi}{M}\right)^{n_{ij}^{f}}\bar{F}_{i}Hf_{j}+\textrm{h.c.}\,,$
(99)

where $F$ is a left-handed $\mathrm{SU}(2)_{L}$ doublet and $f$ is a right-
handed singlet, we now compute the $\rho$ and $a$ couplings to the fermions.
After SSB, we write

$\Phi=\left(\frac{v_{\Phi}+\rho}{\sqrt{2}}\right)e^{ia/v_{\Phi}},\qquad
H=\begin{pmatrix}0\\\ \frac{v_{\mathrm{EW}}+h}{\sqrt{2}}\end{pmatrix}\,.$
(100)

Plugging this into the Lagrangian above and expanding the obtained expression
to leading order in $\rho/v_{\Phi}$ and $a/v_{\Phi}$, we obtain

$\mathcal{L}\supset-
y_{ij}^{f}\left(\frac{v_{\Phi}}{\sqrt{2}M}\right)^{n_{ij}^{f}}\left(1+{n_{ij}^{f}}\frac{\rho+ia}{v_{\Phi}}\right)\bar{f}_{Li}\left(\frac{v_{\mathrm{EW}}+h}{\sqrt{2}}\right)f_{Rj}+\textrm{h.c.}+\textrm{h.o.}\,.$
(101)

where $f_{L}$ is the appropriate component from the doublet $F$. For up-type
quarks, $\tilde{H}$ would be used instead. We define the FN expansion
parameter and the fermion mass matrices

$\epsilon=\frac{v_{\Phi}}{\sqrt{2}M}\,,\qquad
m_{ij}^{f}=\frac{v_{\mathrm{EW}}}{\sqrt{2}}y_{ij}^{f}\epsilon^{n_{ij}^{f}}\,,$
(102)

seeing the FN mechanism at work once again. With this, we can write the
Lagrangian as

$\mathcal{L}\supset-\sum_{f=u,d,e}\left[m_{ij}^{f}\bar{f}_{Li}\left(1+\frac{h}{v_{\mathrm{EW}}}\right)f_{Rj}+\frac{m_{ij}^{f}n_{ij}^{f}}{v_{\Phi}}\bar{f}_{Li}\left(\rho+ia\right)f_{Rj}\right]+\textrm{h.c.}\,.$
(103)

Next, the fermion masses should be diagonalized as usual

$f_{R}\rightarrow U_{R}^{f}f_{R}\,,\quad f_{L}\rightarrow
U_{L}^{f}f_{L}\,,\quad U_{L}^{f\dagger}m^{f}U_{R}^{f}=\hat{m}^{f}\,,$ (104)

where $\hat{m}^{f}$ is a diagonal matrix, $\hat{m}_{ij}^{f}\equiv
m_{i}^{f}\delta_{ij}$. With this, the Higgs-Yukawa couplings are also
diagonalized, whereas the $a$ and $\rho$ couplings are not. This induces, in
general, flavor-changing neutral currents. We then have in the fermion mass
basis

$\mathcal{L}\supset-\sum_{f=u,d,e}\left[m_{i}^{f}\bar{f}_{Li}\left(1+\frac{h}{v_{\mathrm{EW}}}\right)f_{Ri}+c_{ij}^{f}\bar{f}_{Li}\left(\rho+ia\right)f_{Rj}\right]+\textrm{h.c.}\,.$
(105)

with

$c_{ij}^{f}=\frac{1}{v_{\Phi}}(U_{L}^{f\dagger})_{ia}m_{ab}^{f}n_{ab}^{f}U_{Rbj}^{f}\,.$
(106)

We can further rewrite this by noting that $n_{ab}^{f}=[F_{a}]-[f_{b}]$ is a
matrix built from (integer) FN charges. By defining diagonal charge matrices
$Q^{f}_{ij}\equiv[f_{i}]\delta_{ij}$, the combination of matrices appearing in
the middle can be expressed as

$m_{ab}^{f}n_{ab}^{f}=(Q^{F}m^{f}-m^{f}Q^{f})_{ab}\,.$ (107)

With this, the $c^{f}$ coupling matrices can be expressed as

$\begin{split}c^{f}&=\frac{1}{v_{\Phi}}\left[U_{L}^{f\dagger}\left(Q^{F}m^{f}-m^{f}Q^{f}\right)U_{R}^{f}\right]\\\
&=\frac{1}{v_{\Phi}}\left[U_{L}^{f\dagger}Q^{F}\left(U_{L}^{f}U_{L}^{f\dagger}\right)m^{f}U_{R}^{f}-U_{L}^{f\dagger}m^{f}\left(U_{R}^{f}U_{R}^{f\dagger}\right)Q^{f}U_{R}^{f}\right]\\\
&=\frac{1}{v_{\Phi}}\left[\left(U_{L}^{f\dagger}Q^{F}U_{L}^{f}\right)\hat{m}^{f}-\hat{m}^{f}\left(U_{R}^{f\dagger}Q^{f}U_{R}^{f}\right)\right]~{}.\end{split}$
(108)

It is worth reemphasizing that this is now expressed in terms of diagonal
fermion mass matrices $\hat{m}^{f}$ and diagonal FN charge matrices $Q^{F}$
and $Q^{f}$ for the $\mathrm{SU}(2)$ doublet and singlet, respectively.

## References

  * (1) N. Cabibbo, _Unitary Symmetry and Leptonic Decays_ , _Phys. Rev. Lett._ 10 (1963) 531.
  * (2) M. Kobayashi and T. Maskawa, _CP Violation in the Renormalizable Theory of Weak Interaction_ , _Prog. Theor. Phys._ 49 (1973) 652.
  * (3) C.D. Froggatt and H.B. Nielsen, _Hierarchy of Quark Masses, Cabibbo Angles and CP Violation_ , _Nucl. Phys. B_ 147 (1979) 277.
  * (4) M. Leurer, Y. Nir and N. Seiberg, _Mass matrix models_ , _Nucl. Phys. B_ 398 (1993) 319 [hep-ph/9212278].
  * (5) M. Leurer, Y. Nir and N. Seiberg, _Mass matrix models: The Sequel_ , _Nucl. Phys. B_ 420 (1994) 468 [hep-ph/9310320].
  * (6) M. Fedele, A. Mastroddi and M. Valli, _Minimal Froggatt-Nielsen textures_ , _JHEP_ 03 (2021) 135 [2009.05587].
  * (7) C. Cornella, D. Curtin, E.T. Neil and J.O. Thompson, _Mapping and Probing Froggatt-Nielsen Solutions to the Quark Flavor Puzzle_ , 2306.08026.
  * (8) R.D. Peccei and H.R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * (9) F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279.
  * (10) S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * (11) L. Di Luzio, M. Giannotti, E. Nardi and L. Visinelli, _The landscape of QCD axion models_ , _Phys. Rept._ 870 (2020) 1 [2003.01100].
  * (12) I.G. Irastorza and J. Redondo, _New experimental approaches in the search for axion-like particles_ , _Prog. Part. Nucl. Phys._ 102 (2018) 89 [1801.08127].
  * (13) L. Calibbi, F. Goertz, D. Redigolo, R. Ziegler and J. Zupan, _Minimal axion model from flavor_ , _Phys. Rev. D_ 95 (2017) 095009 [1612.08040].
  * (14) Y. Ema, K. Hamaguchi, T. Moroi and K. Nakayama, _Flaxion: a minimal extension to solve puzzles in the standard model_ , _JHEP_ 01 (2017) 096 [1612.05492].
  * (15) J. Martin Camalich, M. Pospelov, P.N.H. Vuong, R. Ziegler and J. Zupan, _Quark Flavor Phenomenology of the QCD Axion_ , _Phys. Rev. D_ 102 (2020) 015023 [2002.04623].
  * (16) H.M. Georgi, L.J. Hall and M.B. Wise, _Grand Unified Models With an Automatic Peccei-Quinn Symmetry_ , _Nucl. Phys. B_ 192 (1981) 409.
  * (17) M. Dine and N. Seiberg, _String Theory and the Strong CP Problem_ , _Nucl. Phys. B_ 273 (1986) 109.
  * (18) S.M. Barr and D. Seckel, _Planck scale corrections to axion models_ , _Phys. Rev. D_ 46 (1992) 539.
  * (19) M. Kamionkowski and J. March-Russell, _Planck scale physics and the Peccei-Quinn mechanism_ , _Phys. Lett. B_ 282 (1992) 137 [hep-th/9202003].
  * (20) R. Holman, S.D.H. Hsu, T.W. Kephart, E.W. Kolb, R. Watkins and L.M. Widrow, _Solutions to the strong CP problem in a world with gravity_ , _Phys. Lett. B_ 282 (1992) 132 [hep-ph/9203206].
  * (21) A. Greljo and A.E. Thomsen, _Rising through the ranks: flavor hierarchies from a gauged $SU(2)$ symmetry_, _Eur. Phys. J. C_ 84 (2024) 213 [2309.11547].
  * (22) A. Greljo, A.E. Thomsen and H. Tiblom, _Flavor Hierarchies From SU(2) Flavor and Quark-Lepton Unification_ , 2406.02687.
  * (23) J. Davighi and G. Isidori, _Non-universal gauge interactions addressing the inescapable link between Higgs and flavour_ , _JHEP_ 07 (2023) 147 [2303.01520].
  * (24) A. Greljo and G. Isidori, _Neutrino Anarchy from Flavor Deconstruction_ , 2406.01696.
  * (25) S. Antusch, A. Greljo, B.A. Stefanek and A.E. Thomsen, _U(2) Is Right for Leptons and Left for Quarks_ , _Phys. Rev. Lett._ 132 (2024) 151802 [2311.09288].
  * (26) T. Higaki and J. Kawamura, _A low-scale flavon model with a $Z_{N}$ symmetry_, _JHEP_ 03 (2020) 129 [1911.09127].
  * (27) E. Palti, _The Swampland: Introduction and Review_ , _Fortsch. Phys._ 67 (2019) 1900037 [1903.06239].
  * (28) L.M. Krauss and F. Wilczek, _Discrete Gauge Symmetry in Continuum Theories_ , _Phys. Rev. Lett._ 62 (1989) 1221.
  * (29) T. Banks, _Effective Lagrangian Description of Discrete Gauge Symmetries_ , _Nucl. Phys. B_ 323 (1989) 90.
  * (30) J. Preskill and L.M. Krauss, _Local Discrete Symmetry and Quantum Mechanical Hair_ , _Nucl. Phys. B_ 341 (1990) 50.
  * (31) L.E. Ibanez and G.G. Ross, _Discrete gauge symmetry anomalies_ , _Phys. Lett. B_ 260 (1991) 291.
  * (32) B. Gavela, P. Quílez and M. Ramos, _The QCD axion sum rule_ , _JHEP_ 04 (2024) 056 [2305.15465].
  * (33) F. Björkeroth, E.J. Chun and S.F. King, _Flavourful Axion Phenomenology_ , _JHEP_ 08 (2018) 117 [1806.00660].
  * (34) M. Bauer, M. Neubert, S. Renner, M. Schnubel and A. Thamm, _Axionlike Particles, Lepton-Flavor Violation, and a New Explanation of $a_{\mu}$ and $a_{e}$_, _Phys. Rev. Lett._ 124 (2020) 211803 [1908.00008].
  * (35) M. Bauer, M. Neubert, S. Renner, M. Schnubel and A. Thamm, _Flavor probes of axion-like particles_ , _JHEP_ 09 (2022) 056 [2110.10698].
  * (36) J. Bonilla, A. de Giorgi, B. Gavela, L. Merlo and M. Ramos, _The cost of an ALP solution to the neutral B-anomalies_ , _JHEP_ 02 (2023) 138 [2209.11247].
  * (37) C. Cornella, A.M. Galda, M. Neubert and D. Wyler, _K ± → $\pi$±a at next-to-leading order in chiral perturbation theory and updated bounds on ALP couplings_, _JHEP_ 06 (2024) 029 [2308.16903].
  * (38) A. Biekötter, J. Fuentes-Martín, A.M. Galda and M. Neubert, _A global analysis of axion-like particle interactions using SMEFT fits_ , _JHEP_ 09 (2023) 120 [2307.10372].
  * (39) L. Calibbi, T. Li, L. Mukherjee and Y. Yang, _Probing ALP Lepton Flavour Violation at $\mu$TRISTAN_, 2406.13234.
  * (40) L. Di Luzio, A.W.M. Guerrera, X. Ponce Díaz and S. Rigolin, _Axion-Like Particles in Radiative Quarkonia Decays_ , 2402.12454.
  * (41) S. Knapen, K. Langhoff, T. Opferkuch and D. Redigolo, _A Robust Search for Lepton Flavour Violating Axions at Mu3e_ , 2311.17915.
  * (42) Flavour Lattice Averaging Group (FLAG) collaboration, _FLAG Review 2021_ , _Eur. Phys. J. C_ 82 (2022) 869 [2111.09849].
  * (43) N. Gubernari, M. Reboud, D. van Dyk and J. Virto, _Dispersive analysis of $B\to K^{(*)}$ and $B_{s}\to\phi$ form factors_, _JHEP_ 12 (2023) 153 [2305.06301].
  * (44) HPQCD collaboration, _$B\rightarrow K$ and $D\rightarrow K$ form factors from fully relativistic lattice QCD_, _Phys. Rev. D_ 107 (2023) 014510 [2207.12468].
  * (45) J.A. Bailey et al., _$B\to Kl^{+}l^{-}$ Decay Form Factors from Three-Flavor Lattice QCD_, _Phys. Rev. D_ 93 (2016) 025026 [1509.06235].
  * (46) HPQCD collaboration, _Rare decay $B\to K\ell^{+}\ell^{-}$ form factors from lattice QCD_, _Phys. Rev. D_ 88 (2013) 054509 [1306.2384].