On the IR/UV flavour connection
in non-universal axion models

Luca Di Luzio, Alfredo Walter Mario Guerrera,
Xavier Ponce Díaz, Stefano Rigolin

Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Padova,
Via F. Marzolo 8, 35131 Padova, Italy
Dipartimento di Fisica e Astronomia ‘G. Galilei’, Università di Padova,
Via F. Marzolo 8, 35131 Padova, Italy

> Non-universal axion models, with the Peccei-Quinn (PQ) symmetry acting on
> Standard Model (SM) fermions in a generation-dependent way, are typically
> accompanied by two different sources of flavour violation, dubbed here as
> infrared (IR) and ultraviolet (UV). The former is due to the flavour
> violating axion couplings to SM fermions, while the latter arises from the
> heavy degrees of freedom that UV complete the axion effective field theory.
> We point out that these two sources of flavour violation are directly
> related and exemplify this connection in a general class of non-universal
> axion model, based on a renormalizable DFSZ-like setup with two Higgs
> doublets (PQ-2HDM). We next discuss the interplay of axion flavour
> phenomenology with the signatures stemming from the heavy radial modes of
> the PQ-2HDM, including meson oscillation observables and charged lepton
> flavour violating decays. We emphasize the strong complementarity between
> flavour observables, LHC direct searches and standard axion physics.

###### Contents

  1. 1 Introduction
  2. 2 Non-universal DFSZ axion models
    1. 2.1 Flavour violation from the light axion field (IR source)
    2. 2.2 Flavour violation from the heavy radial modes (UV source)
    3. 2.3 Connection between IR/UV sources of flavour violation
    4. 2.4 Non-universal PQ-2HDM Yukawa textures
  3. 3 Flavour phenomenology
    1. 3.1 LHC bounds on the 2HDM
    2. 3.2 Flavour constraints on the non-universal 2HDM
    3. 3.3 IR/UV quark flavour connection
    4. 3.4 IR/UV lepton flavour connection
  4. 4 Conclusions
  5. A PQ-2HDM scalar potential
  6. B 4-fermion operators
  7. C Meson mixing bounds
  8. D Model-dependent bounds on the PQ-2HDM
  9. E Lepton flavour violating observables

##  1 Introduction

The Peccei-Quinn (PQ) solution to the strong CP problem [1, 2] has emerged as
a well-motivated paradigm beyond the Standard Model (SM), by delivering in the
low-energy effective field theory (EFT) the axion field [3, 4] which also
provides a natural dark matter candidate [5, 6, 7]. Invisible axion models can
be classified in two broad categories, known conventionally as DFSZ [8, 9] and
KSVZ [10, 11] models, depending on whether the SM fermions are charged or not
under the ${\rm U}(1)_{\rm PQ}$ symmetry. In the former case, only the scalar
sector of the SM is extended, by a second Higgs doublet and a SM singlet that
allows to establish the hierarchy between the PQ and electroweak breaking
scales. In the original formulation of the DFSZ model [8, 9] the SM fermions
PQ-charge assignment is flavour blind, namely SM fermions with the same
hypercharge share also the same PQ charge. On the other hand, non-universal
variants of the DFSZ model (and in particular of the original WW axion model
[3, 4]), in which the PQ symmetry acts on SM fermions in a generation-
dependent way, appeared in the late 80’s [12, 13, 14], in the effort to
explain in terms of a heavy axion the anomalous observation of a narrow
$e^{+}e^{-}$ pair peak at 1.8 MeV in heavy ion collision at the GSI [15, 16].

From a modern perspective, non-universal axion models can be motivated in
several respects: $i)$ their potential connection to the SM flavour puzzle,
namely explaining the observed pattern of SM fermion masses and mixings in
terms of a non-universal ${\rm U}(1)_{\rm PQ}$ acting as a flavour symmetry
(see e.g. [17, 18, 19, 20]); $ii)$ the possibility of suppressing the axion
couplings to nucleons and electrons, leading to the so-called “astrophobic”
axion scenario [21] (see also [22, 23, 24, 25, 26]), which allows to relax
astrophysical bounds on the axion decay constant; $iii)$ the experimental
opportunity of discovering the axion via flavoured axion searches (see e.g.
[27, 28, 29, 30, 31, 32, 33, 34]), which offer a complementary probe with
respect to standard axion detection strategies [35, 36].

One of the most striking consequences of non-universal axion models is that,
after going to the mass basis, fermion rotation matrices beyond the CKM/PMNS
combinations become physical and thus can be probed via flavour violating
axion couplings to SM fermions. We will indicate this type of flavour
violation as of low-energy or infrared (IR) type, in the sense that it stems
from the light axion field. However, this is not the only flavour violation to
be expected, since the ultraviolet (UV) completion of the axion EFT might
leave behind extra sources of flavour violation, indicated in the following as
of high-energy or UV type. While we expect this to be a generic feature of
axion EFTs, the scope of the present paper is to exemplify this connection in
a general class of renormalizable (i.e. fully calculable) non-universal DFSZ-
like models, whose scalar sector is based on a SM singlet and two Higgs
doublets which couple non-universally to different generations of up- and
down-type quarks (and charged leptons). Since the SM fermion mass matrices and
the Yukawa interactions of the two Higgs doublet model (2HDM) are non-
simultaneously diagonalizable, new flavour violating UV sources arise from the
Yukawa couplings of the 2HDM degrees of freedom. Note that the mass scale of
the non-SM-like Higgs scalars is a free parameter, ranging from a few hundreds
of GeV (experimentally bounded by LHC direct searches) up to the PQ breaking
scale, of the order of the axion decay constant, $f_{a}$. In this paper, we
will entertain the bottom-up perspective that the 2HDM is not completely
decoupled from the TeV scale,111A low-scale 2HDM within a PQ setup is
theoretically motivated by the possibility of establishing the hierarchy
between the electroweak and PQ scales in a “technically natural way”, i.e.
stable under radiative corrections thanks to the presence of “ultra-weak”
couplings in the scalar potential (see Appendix A for further details). and
thus its degrees of freedom can be directly tested at the LHC and indirectly
via flavour violating processes.

The aim of this paper, is to underline the general correlation between IR and
UV sources of flavour violation in non-universal DFSZ-like axion models,
focusing on flavour violating observables both in the quark and lepton
sector.222For a related work focusing on flavour violation in the charged-
lepton sector, see Ref. [37]. In particular, we will show that the pattern of
flavour violation involving the light axion field is directly correlated to
that of the heavy radial modes of the 2HDM. Assuming for instance a QCD axion
discovery in a golden channel like e.g. $K\to\pi a$ (which would definitely
point to a tree-level source of flavour violation in $s$-$d$ axion couplings),
allows us to predict a definite pattern of flavour violation from the 2HDM
heavy scalars in $s$-$d$ transitions, e.g. in $K^{0}$-$\overline{K}^{0}$
mixing. The bound arising from the latter observable can be then confronted
with LHC direct searches, as well as with astrophysical bounds on the axion
decay constant. Similar considerations apply also to other phenomenologically
relevant quark flavour transitions involving $B$ and $D$ mesons or to lepton
flavour violating (LFV) decays.

The paper is structured as follows. Section 2 focuses on theoretical aspects
of the IR/UV flavour connection. We first introduce a general class of non-
universal axion models (the only assumption being here a minimal field content
and renormalizability) and recall the pattern of flavour violating axion
couplings in terms of PQ charges. Next, we discuss the structure of a general
flavour non-conserving 2HDM in the presence of a ${\rm U}(1)_{\rm PQ}$
symmetry (hereby denoted as PQ-2HDM), and work out the flavour violating
couplings of the 2HDM scalar fields and the associated 4-fermion operators
that are generated after integrating out the heavy radial modes. A simplified
class of non-universal axion models is also introduced, in order to have an
explicit realization as a benchmark scenario. Phenomenological aspects,
including LHC bounds and flavour constraints on the non-universal 2HDM are
discussed in Section 3, together with the aforementioned interplay between
flavoured axion searches and other more standard flavour violating processes,
both in the hadronic and leptonic sectors. We summarize our conclusions in
Section 4, while more technical details about the scalar potential and flavour
violating observables are deferred to Appendices A–E.

##  2 Non-universal DFSZ axion models

The scalar sector of the DFSZ model includes two Higgs doublets
$H_{1,2}\sim(1,2,-1/2)$ endowed with PQ charges $\mathcal{X}_{1,2}$ and an
extra complex SM-singlet $\phi\sim(1,1,0)$, with PQ charge
$\mathcal{X}_{\phi}$, needed for decoupling the PQ and electroweak breaking
scales. Singlet and doublet scalar fields are coupled via a single non-
hermitian operator, $H_{2}^{\dagger}H_{1}\phi$, thus implying
$\mathcal{X}_{2}-\mathcal{X}_{1}=\mathcal{X}_{\phi}=1$, by normalizing to 1
the PQ charge of the singlet.333Other choices, e.g.
$H_{2}^{\dagger}H_{1}\phi^{2}$ as in the original DFSZ model [8, 9], are
possible as well at the renormalizable level. The general conclusions of this
paper do not depend on this particular choice, which affects instead the
periodicity of the axion potential via the so-called domain-wall number (see
e.g. [38]). In the following, we will also require the orthogonality of the PQ
and hypercharge currents (see e.g. [38]), thus ensuring the absence of a mixed
$Z^{\mu}\partial_{\mu}a$ kinetic term. This fixes the PQ charges of the two
Higgs doublets as $\mathcal{X}_{1}=-s^{2}_{\beta}\equiv-\sin^{2}\beta$ and
$\mathcal{X}_{2}=c^{2}_{\beta}\equiv\cos^{2}\beta$, in terms of the vacuum
angle, $t_{\beta}\equiv\tan\beta=v_{2}/v_{1}$, where $v_{1,2}$ are the
$H_{1,2}$ vacuum expectation values (VEVs) with
$v^{2}=v_{1}^{2}+v_{2}^{2}\simeq(246\ \text{GeV})^{2}$. Other details
concerning the renormalizable DFSZ scalar potential, the implementation of the
hierarchy between the PQ and the electroweak scales and the resulting scalar
spectrum are deferred to Appendix A.

The Yukawa sector of the 2HDM in the presence of a ${\rm U}(1)_{\rm PQ}$
symmetry can be written in a compact form as

$\displaystyle\mathcal{L}_{Y}^{\rm PQ\text{-}2HDM}=$
$\displaystyle-\overline{q}_{L}(Y_{1}^{u}H_{1}+Y_{2}^{u}H_{2})u_{R}+\overline{q}_{L}(Y_{1}^{d}\tilde{H}_{1}+Y_{2}^{d}\tilde{H}_{2})d_{R}$
$\displaystyle+\overline{\ell}_{L}(Y_{1}^{e}\tilde{H}_{1}+Y_{2}^{e}\tilde{H}_{2})e_{R}+\text{h.c.}\,,$
(2.1)

where $\tilde{H}_{1,2}=(i\sigma^{2})H_{1,2}^{*}$ and $Y^{f}_{1,2}$ are generic
three-dimensional matrices. We do not specify the origin of neutrino masses,
the effects of which will be neglected in the present discussion. Requiring
the ${\rm U}(1)_{\rm PQ}$ invariance on the Yukawa Lagrangian implies the
following sum rules in flavour space:

$\displaystyle-\mathcal{X}_{q}Y_{1,2}^{u}+Y_{1,2}^{u}\mathcal{X}_{u}+\mathcal{X}_{1,2}Y_{1,2}^{u}$
$\displaystyle=0\,,$ (2.2)
$\displaystyle-\mathcal{X}_{q}Y_{1,2}^{d}+Y_{1,2}^{d}\mathcal{X}_{d}-\mathcal{X}_{1,2}Y_{1,2}^{d}$
$\displaystyle=0\,,$ (2.3)
$\displaystyle-\mathcal{X}_{\ell}Y_{1,2}^{e}+Y_{1,2}^{e}\mathcal{X}_{e}-\mathcal{X}_{1,2}Y_{1,2}^{e}$
$\displaystyle=0\,,$ (2.4)

where the left ($\mathcal{X}_{q,\ell}$) and right ($\mathcal{X}_{u,d,e}$) SM
fermion PQ charges are diagonal, but in general not universal, three-
dimensional matrices. Note that the PQ symmetry requires some of the Yukawa
entries to vanish. For instance, in the universal DFSZ framework [8, 9] one
imposes $Y_{1}^{u}=Y_{2}^{d,e}=0$, falling into a type-II 2HDM scenario. In
Section 2.4 we will provide two explicit examples of non-universal PQ-2HDM
Yukawa textures.

###  2.1 Flavour violation from the light axion field (IR source)

The low-energy axion field, which spans over the neutral angular modes of the
scalar multiplets

$\phi\supset\frac{v_{\phi}}{\sqrt{2}}e^{ia_{\phi}/v_{\phi}}\,,\quad
H_{1}\supset\frac{v_{1}}{\sqrt{2}}e^{ia_{1}/v_{1}}\begin{pmatrix}1\\\
0\end{pmatrix}\,,\quad
H_{2}\supset\frac{v_{2}}{\sqrt{2}}e^{ia_{2}/v_{2}}\begin{pmatrix}1\\\
0\end{pmatrix}\,,$ (2.5)

is defined as [38]

$a=\frac{1}{v_{\rm
PQ}}\left(v_{\phi}a_{\phi}+v(-s^{2}_{\beta}c_{\beta}a_{1}+c^{2}_{\beta}s_{\beta}a_{2})\right)\simeq
a_{\phi}\,,\qquad v_{\rm
PQ}^{2}=v_{\phi}^{2}+v^{2}c^{2}_{\beta}s^{2}_{\beta}\simeq v_{\phi}^{2}\,,$
(2.6)

where the $\simeq$ relation holds once $v_{\rm\phi}\gg v$ is imposed. The
axion effective Lagrangian below the electroweak scale reads [38]

$\mathcal{L}_{a}\supset\frac{g_{s}^{2}}{32\pi^{2}}\frac{a}{f_{a}}G\tilde{G}\,+\,\frac{E}{N}\frac{e^{2}}{32\pi^{2}}\frac{a}{f_{a}}F\tilde{F}\,+\\!\\!\sum_{f=u,d,e}\frac{\partial_{\mu}a}{2f_{a}}\overline{f}_{i}\gamma^{\mu}\left((C^{V}_{f})_{ij}+(C^{A}_{f})_{ij}\gamma_{5}\right)f_{j}\,,$
(2.7)

where $E/N$ is the ratio between the electromagnetic and QCD anomaly of the PQ
current (for typical values in concrete axion models, see e.g. [39, 40]), and
the axion decay constant is defined by $f_{a}=v_{\rm PQ}/(2N)$. The axion
vector and axial couplings to fermions can be written as

$C^{V,A}_{f}=C^{R}_{f}\pm C^{L}_{f}\,,\qquad
C^{R,L}_{f}=\frac{1}{2N}V_{f_{R,L}}^{\dagger}\mathcal{X}_{f_{R,L}}V_{f_{R,L}}\,\qquad(f=u,d,e)\,,$
(2.8)

in terms of the PQ charges $\mathcal{X}_{f_{L,R}}$ and of the unitary
matrices, $V_{f_{L,R}}$, which diagonalize the SM fermion mass terms:

$\hat{M}_{f}\equiv V^{\dagger}_{f_{L}}M_{f}V_{f_{R}}\,.$ (2.9)

Therefore, the $C_{f}^{V,A}$ parameters encode all the IR sources of flavour
violation, stemming from the light axion field. In the present work we will
neglect renormalization group corrections on the $C_{f}^{V,A}$ couplings (see
[41, 42, 43, 44, 24]). This is a reasonable approximation under the assumption
of a (at most) TeV-scale 2HDM, being running effects proportional to $\log
m_{H}/m_{t}$, with $m_{H}$ labelling the typical mass scale of the 2HDM radial
modes in the decoupling limit (see Appendix A for a detailed discussion of the
scalar spectrum).

###  2.2 Flavour violation from the heavy radial modes (UV source)

The SM fermion mass matrices in the PQ-2HDM follow from Eq. (2) after setting
the Higgs fields on their VEVs, and read

$\displaystyle
M_{u}=\frac{v}{\sqrt{2}}c_{\beta}Y^{u}_{1}+\frac{v}{\sqrt{2}}s_{\beta}Y^{u}_{2}\,,\
\
M_{d}=\frac{v}{\sqrt{2}}c_{\beta}Y^{d}_{1}+\frac{v}{\sqrt{2}}s_{\beta}Y^{d}_{2}\,,\
\
M_{e}=\frac{v}{\sqrt{2}}c_{\beta}Y^{e}_{1}+\frac{v}{\sqrt{2}}s_{\beta}Y^{e}_{2}\,.$
(2.10)

Since in non-universal PQ models the mass matrices and the Yukawa interactions
are not simultaneously diagonalizable, flavour violating Higgs couplings to
fermions are generated at tree level both in the neutral and charged sector.
To derive those interactions, we first decompose the two Higgs doublets in
terms of mass and charge eigenstates, comprising $h$ (the SM-like Higgs) and
$H,A,H^{+}$ (the heavy radial modes). In the unitary gauge one has (see
Appendix A)

$H_{1}=\begin{pmatrix}\frac{1}{\sqrt{2}}(c_{\beta}v-s_{\alpha}h+c_{\alpha}H-is_{\beta}A)\\\
-s_{\beta}H^{-}\end{pmatrix}\,,\quad
H_{2}=\begin{pmatrix}\frac{1}{\sqrt{2}}(s_{\beta}v+c_{\alpha}h+s_{\alpha}H+ic_{\beta}A)\\\
c_{\beta}H^{-}\end{pmatrix}\,.$ (2.11)

After unfolding the ${\rm SU}(2)_{L}$ structure in Eq. (2) and projecting onto
the fermion mass basis, $f_{L,R}\to V_{f_{L,R}}f_{L,R}$ (with $f=u,d,e$), we
obtain the following Yukawa Lagrangian

$\displaystyle\mathcal{L}_{Y}^{\rm PQ\text{-}2HDM}\supset$
$\displaystyle-\sum_{f=u,d,e}\overline{f}_{Li}C^{h_{f}}_{ij}f_{Rj}h-\sum_{f=u,d,e}\overline{f}_{Li}C^{H_{f}}_{ij}f_{Rj}H-\sum_{f=u,d,e}\overline{f}_{Li}C^{A_{f}}_{ij}f_{Rj}A$
(2.12)
$\displaystyle-\sqrt{2}\left(\overline{d}_{Li}V^{\dagger}_{ik}C^{H^{-}_{u}}_{kj}u_{Rj}H^{-}+\overline{u}_{Li}V_{ik}C^{H^{+}_{d}}_{kj}d_{Rj}H^{+}+\overline{\nu}_{Li}U^{\dagger}_{ik}C^{H^{+}_{e}}_{kj}e_{Rj}H^{+}\right)+\text{h.c.}\,,$

where the CKM and PMNS mixing matrix entering in the charged-current
interactions are defined respectively as $V\equiv V_{\rm
CKM}=V^{\dagger}_{u_{L}}V_{d_{L}}$ and $U\equiv U_{\rm
PMNS}=V^{\dagger}_{e_{L}}V_{\nu_{L}}$, and the Higgs flavour violating
coefficients $C_{ij}$ read (see also [37])

$\displaystyle C^{h_{u}}_{ij}$
$\displaystyle=\frac{m_{u_{i}}}{vs_{\beta}}\delta_{ij}c_{\alpha}+\frac{1}{\sqrt{2}}\frac{-c_{\alpha-\beta}}{s_{\beta}}\epsilon^{u}_{ij}\,,\qquad\,\,C^{h_{d}(h_{e})}_{ij}=\frac{m_{d_{i}(e_{i})}}{vc_{\beta}}\delta_{ij}(-s_{\alpha})+\frac{1}{\sqrt{2}}\frac{c_{\alpha-\beta}}{c_{\beta}}\epsilon^{d(e)}_{ij}\,,$
(2.13) $\displaystyle C^{H_{u}}_{ij}$
$\displaystyle=\frac{m_{u_{i}}}{vs_{\beta}}\delta_{ij}s_{\alpha}+\frac{1}{\sqrt{2}}\frac{-s_{\alpha-\beta}}{s_{\beta}}\epsilon^{u}_{ij}\,,\quad\
\
\,C^{H_{d}(H_{e})}_{ij}=\frac{m_{d_{i}(e_{i})}}{vc_{\beta}}\delta_{ij}c_{\alpha}+\frac{1}{\sqrt{2}}\frac{s_{\alpha-\beta}}{c_{\beta}}\epsilon^{d(e)}_{ij}\,,$
(2.14) $\displaystyle C^{A_{u}}_{ij}$
$\displaystyle=\frac{m_{u_{i}}}{vs_{\beta}}\delta_{ij}(ic_{\beta})+\frac{1}{\sqrt{2}}\frac{-i}{s_{\beta}}\epsilon^{u}_{ij}\,,\qquad\
\
C^{A_{d}(A_{e})}_{ij}=\frac{m_{d_{i}(e_{i})}}{vc_{\beta}}\delta_{ij}(is_{\beta})+\frac{1}{\sqrt{2}}\frac{-i}{c_{\beta}}\epsilon^{d(e)}_{ij}\,,$
(2.15) $\displaystyle C^{H^{-}_{u}}_{ij}$
$\displaystyle=\frac{m_{u_{i}}}{vs_{\beta}}\delta_{ij}c_{\beta}-\frac{1}{\sqrt{2}}\frac{\epsilon^{u}_{ij}}{s_{\beta}}\,,\qquad\quad\
\ \
C^{H^{+}_{d}(H^{+}_{e})}_{ij}=\frac{m_{d_{i}(e_{i})}}{vc_{\beta}}\delta_{ij}s_{\beta}-\frac{1}{\sqrt{2}}\frac{\epsilon^{d(e)}_{ij}}{c_{\beta}}\,,$
(2.16)

with

$\epsilon^{u}_{ij}=(V^{\dagger}_{u_{L}}Y_{1}^{u}V_{u_{R}})_{ij}\,,\qquad\epsilon^{d}_{ij}=(V^{\dagger}_{d_{L}}Y_{2}^{d}V_{d_{R}})_{ij}\,,\qquad\epsilon^{e}_{ij}=(V^{\dagger}_{e_{L}}Y_{2}^{e}V_{e_{R}})_{ij}\,.$
(2.17)

Therefore, all the neutral and charged flavour violating sources in the Higgs
sector are encoded into the $\epsilon^{f}$ parameters. Note that, besides the
heavy Higgs radial modes, also the SM-like Higgs features flavour violating
couplings (cf. Eq. (2.13)). However, all the tree-level flavour changing
neutral current contributions arising from the light Higgs sector decouple
faster than the radial mode ones, since $c_{\alpha-\beta}\propto 1/m_{H}^{2}$
in the decoupling limit (cf. Eq. (A.28)).

To connect with typical low-energy flavour observables, like for example meson
oscillations, one has to integrate out both the heavy radial modes,
$H,A,H^{+}$ and the SM–like Higgs, $h$, obtaining the $d=6$ low-energy
effective Lagrangian describing the (light) quarks and leptons 4-fermion
interactions

$\mathcal{L}^{\rm 4\text{-}fermion}_{\rm EFT}=\mathcal{L}^{\rm
4\text{-}quark}_{\rm EFT}+\mathcal{L}^{\rm semi\text{-}lept}_{\rm
EFT}+\mathcal{L}^{\rm 4\text{-}lept}_{\rm EFT}\,.$ (2.18)

All the operators appearing in the 4-quark, semi-leptonic and 4-lepton
sectors, together with the expressions of the corresponding Wilson
coefficients as a function of the PQ-2HDM parameters are listed in Appendix B.

###  2.3 Connection between IR/UV sources of flavour violation

The IR sources of flavour violation, $C^{V,A}_{u,d,e}$, involving the light
axion field and the UV sources, $\epsilon^{u,d,e}$, involving the heavy radial
modes of the PQ-2HDM, turn out to be related thanks to Eqs. (2.2)–(2.4), which
encode the action of the ${\rm U}(1)_{\rm PQ}$ symmetry on the Yukawa sector
of the PQ-2HDM. By means of the SM fermion masses in Eq. (2.10) it is possible
to rewrite the Yukawa matrices in terms of the (non-universal) fermionic PQ
charges as

$\displaystyle Y_{1}^{u}$
$\displaystyle=\frac{\sqrt{2}}{vc_{\beta}}\left(-\mathcal{X}_{q}M_{u}+M_{u}\mathcal{X}_{u}+\mathcal{X}_{2}M_{u}\right)\,,\quad
Y_{2}^{u}=\frac{\sqrt{2}}{vs_{\beta}}\left(\mathcal{X}_{q}M_{u}-M_{u}\mathcal{X}_{u}-\mathcal{X}_{1}M_{u}\right)\,,$
(2.19) $\displaystyle Y_{1}^{d}$
$\displaystyle=\frac{\sqrt{2}}{vc_{\beta}}\left(\mathcal{X}_{q}M_{d}-M_{d}\mathcal{X}_{d}+\mathcal{X}_{2}M_{d}\right)\,,\quad\
\ \
Y_{2}^{d}=\frac{\sqrt{2}}{vs_{\beta}}\left(-\mathcal{X}_{q}M_{d}+M_{d}\mathcal{X}_{d}-\mathcal{X}_{1}M_{d}\right)\,,$
(2.20) $\displaystyle Y_{1}^{e}$
$\displaystyle=\frac{\sqrt{2}}{vc_{\beta}}\left(\mathcal{X}_{\ell}M_{e}-M_{e}\mathcal{X}_{e}+\mathcal{X}_{2}M_{e}\right)\,,\quad\
\ \
Y_{2}^{e}=\frac{\sqrt{2}}{vs_{\beta}}\left(-\mathcal{X}_{\ell}M_{e}+M_{e}\mathcal{X}_{e}-\mathcal{X}_{1}M_{e}\right)\,,$
(2.21)

where the PQ charge relation $\mathcal{X}_{2}-\mathcal{X}_{1}=1$ arising from
the scalar potential has also been used. Hence, by exploiting the definitions
of $C^{L,R}_{f}$ and $\epsilon^{f}$ in Eq. (2.8) and Eq. (2.17), we can write
the following expression relating the flavour violating couplings of the axion
field and the heavy radial modes:

$\displaystyle\epsilon^{u}$
$\displaystyle=\frac{2N\sqrt{2}}{vc_{\beta}}\left(-C_{u}^{L}\hat{M}_{u}+\hat{M}_{u}C_{u}^{R}+\frac{\mathcal{X}_{2}}{2N}\hat{M}_{d}\right)\,,$
(2.22) $\displaystyle\epsilon^{d}$
$\displaystyle=\frac{2N\sqrt{2}}{vs_{\beta}}\left(-C_{d}^{L}\hat{M}_{d}+\hat{M}_{d}C_{d}^{R}-\frac{\mathcal{X}_{1}}{2N}\hat{M}_{d}\right)\,,$
(2.23) $\displaystyle\epsilon^{e}$
$\displaystyle=\frac{2N\sqrt{2}}{vs_{\beta}}\left(-C_{e}^{L}\hat{M}_{e}+\hat{M}_{e}C_{e}^{R}-\frac{\mathcal{X}_{1}}{2N}\hat{M}_{e}\right)\,.$
(2.24)

Eqs. (2.22)–(2.24) provide the aforementioned connection between the IR and UV
sources of flavour violation, and they apply to the most general class of
renormalizable non-universal axion models based on two Higgs doublets. Notice
that these relations do not depend explicitly on the unknown fermion mass
rotation matrices, but only on the physical fermion masses and fermion
couplings to the axion. Moreover, notice that the terms proportional to the
Higgs doublet PQ charges, $\mathcal{X}_{1,2}$, are diagonal and hence they do
not intervene in flavour violating interactions. In addition, from the
hermiticity of the $C_{f}^{L,R}$ matrices, it follows that transposed off-
diagonal entries can be simply related through

$\displaystyle(\epsilon^{f}_{ji})^{*}=\epsilon^{f}_{ij}(m_{f_{i}}\leftrightarrow
m_{f_{j}})\,.$ (2.25)

For a more transparent connection with flavour observables, it is useful to
rewrite Eqs. (2.22)–(2.24) in terms of vector and axial flavour violating
couplings. For instance, the off-diagonal $d$-quark entries read:

$\displaystyle(\epsilon_{A}^{d})_{ij}\equiv\epsilon^{d}_{ij}+\epsilon^{d\,*}_{ji}=\sqrt{2}\,\frac{2N}{v\,s_{\beta}}(m_{d_{i}}+m_{d_{j}})(C^{A}_{d})_{ij}\,,\qquad(i\neq
j)$ (2.26)
$\displaystyle(\epsilon_{V}^{d})_{ij}\equiv\epsilon^{d}_{ij}-\epsilon^{d\,*}_{ji}=\sqrt{2}\,\frac{2N}{v\,s_{\beta}}(m_{d_{i}}-m_{d_{j}})(C^{V}_{d})_{ij}\,.\qquad(i\neq
j)$ (2.27)

The charged-lepton coefficients are simply obtained by replacing $d\to e$ in
Eqs. (2.26)–(2.27), while the up-quark ones can be obtained via the
substitutions $d\to u$ and $s_{\beta}\to c_{\beta}$. In addition, for
describing $\Delta F=2$ observables (cf. Eqs. (C.12)–(C.17)) the following
combinations are going to be useful:

$\displaystyle(\epsilon^{d}_{ij})^{2}+(\epsilon^{d}_{ij})^{*\,2}$
$\displaystyle=\frac{(\epsilon_{A}^{d})^{2}_{ij}+(\epsilon_{V}^{d})^{2}_{ij}}{2}$
$\displaystyle=\left(\frac{2N}{vs_{\beta}}\right)^{2}\left[(m_{d_{i}}+m_{d_{j}})^{2}(C^{A}_{d})^{2}_{ij}+(m_{d_{i}}-m_{d_{j}})^{2}(C^{V}_{d})^{2}_{ij}\right]\,,$
(2.28) $\displaystyle 2\epsilon^{d}_{ij}(\epsilon^{d}_{ij})^{*}$
$\displaystyle=\frac{(\epsilon_{A}^{d})^{2}_{ij}-(\epsilon_{V}^{d})^{2}_{ij}}{2}$
$\displaystyle=\left(\frac{2N}{vs_{\beta}}\right)^{2}\left[(m_{d_{i}}+m_{d_{j}})^{2}(C^{A}_{d})^{2}_{ij}-(m_{d_{i}}-m_{d_{j}})^{2}(C^{V}_{d})^{2}_{ij}\right]\,.$
(2.29)

Again, the analogous expressions for charged leptons are obtained by replacing
$d\to e$, while for the up-quark sector by sending $d\to u$ and $s_{\beta}\to
c_{\beta}$ in Eqs. (2.3)–(2.3).

###  2.4 Non-universal PQ-2HDM Yukawa textures

The relations derived in the previous section hold for a general non-universal
PQ-2HDM. On the other hand, it is useful to provide for the following
phenomenological analysis a couple of benchmark Yukawa textures that can
nicely exemplify the IR/UV flavour correlation. Here, we will consider two
explicit examples of 2+1-type models, i.e. with at least two generations
featuring the same PQ charges, originally presented in Ref. [21] (see also
[45]) in order to suppress the axion couplings to nucleons. Despite the
reasons behind their original motivation, these two textures, labelled
respectively as M1 and M4 in Ref. [21], are prototypes of models with quark
flavour violation in axion couplings, respectively in the left-handed (LH) and
in the right-handed (RH) down-quark sector. The lepton sector is more model-
dependent and we report here for definiteness just one of the possible
choices. Note that in these models it is also possible to achieve
electrophobia via a tuning involving the rotation matrices of the charged
leptons [21], thus realizing the so-called “astrobhobic” axion scenario in
which both coupling to nucleons and electrons are suppressed, and
astrophysical axion bounds are in turn suppressed. The M1 and M4 textures are
defined by the following Yukawa matrices:

  * •

M1 model

$\displaystyle Y^{u}_{1}$
$\displaystyle=\begin{pmatrix}y^{u}_{11}&y^{u}_{12}&y^{u}_{13}\\\
y^{u}_{21}&y^{u}_{22}&y^{u}_{23}\\\ 0&0&0\end{pmatrix}\,,\quad
Y^{u}_{2}=\begin{pmatrix}0&0&0\\\ 0&0&0\\\
y^{u}_{31}&y^{u}_{32}&y^{u}_{33}\end{pmatrix}\,,$ (2.30) $\displaystyle
Y^{d}_{1}$ $\displaystyle=\begin{pmatrix}0&0&0\\\ 0&0&0\\\
y^{d}_{31}&y^{d}_{32}&y^{d}_{33}\end{pmatrix}\,,\quad
Y^{d}_{2}=\begin{pmatrix}y^{d}_{11}&y^{d}_{12}&y^{d}_{13}\\\
y^{d}_{21}&y^{d}_{22}&y^{d}_{23}\\\ 0&0&0\end{pmatrix}\,,$ (2.31)
$\displaystyle Y^{e}_{1}$
$\displaystyle=\begin{pmatrix}y^{e}_{11}&y^{e}_{12}&y^{e}_{13}\\\
y^{e}_{21}&y^{e}_{22}&y^{e}_{23}\\\ 0&0&0\end{pmatrix}\,,\quad
Y^{e}_{2}=\begin{pmatrix}0&0&0\\\ 0&0&0\\\
y^{e}_{31}&y^{e}_{32}&y^{e}_{33}\end{pmatrix}\,.$ (2.32)

Upon a baryon and lepton number transformation one can set
$\mathcal{X}_{q_{\alpha}}=\mathcal{X}_{\ell_{\alpha}}=0$ ($\alpha=1,2$), so
that the PQ charge matrices satisfying the constraints in Eqs. (2.2)–(2.4)
read

$\displaystyle\mathcal{X}_{q}$
$\displaystyle=\mbox{diag}\,(0,0,1)\,,\quad\mathcal{X}_{u}=\mbox{diag}\,(s^{2}_{\beta},s^{2}_{\beta},s^{2}_{\beta})\,,\quad\mathcal{X}_{d}=\mbox{diag}\,(c^{2}_{\beta},c^{2}_{\beta},c^{2}_{\beta})\,,$
$\displaystyle\mathcal{X}_{\ell}$
$\displaystyle=-\mathcal{X}_{q}\,,\quad\quad\quad\quad\mathcal{X}_{e}=-\mathcal{X}_{u}\,.$
(2.33)

which shows that in the M1 model only the PQ charges in the quark and lepton
LH sectors are non-universal.

  * •

M4 model

$\displaystyle Y^{u}_{1}$
$\displaystyle=\begin{pmatrix}y^{u}_{11}&y^{u}_{12}&y^{u}_{13}\\\
y^{u}_{21}&y^{u}_{22}&y^{u}_{23}\\\
y^{u}_{31}&y^{u}_{32}&y^{u}_{33}\end{pmatrix}\,,\quad Y^{u}_{2}=0\,,$ (2.34)
$\displaystyle Y^{d}_{1}$
$\displaystyle=\begin{pmatrix}0&y^{d}_{12}&y^{d}_{13}\\\
0&y^{d}_{22}&y^{d}_{23}\\\ 0&y^{d}_{32}&y^{d}_{33}\end{pmatrix}\,,\quad
Y^{d}_{2}=\begin{pmatrix}y^{d}_{11}&0&0\\\ y^{d}_{21}&0&0\\\
y^{d}_{31}&0&0\end{pmatrix}\,,$ (2.35) $\displaystyle Y^{e}_{1}$
$\displaystyle=\begin{pmatrix}y^{e}_{11}&y^{e}_{12}&y^{e}_{13}\\\ 0&0&0\\\
0&0&0\end{pmatrix}\,,\quad Y^{e}_{2}=\begin{pmatrix}0&0&0\\\
y^{e}_{21}&y^{e}_{22}&y^{e}_{23}\\\
y^{e}_{31}&y^{e}_{32}&y^{e}_{33}\end{pmatrix}\,.$ (2.36)

Setting $\mathcal{X}_{q_{1}}=\mathcal{X}_{\ell_{1}}=0$ via a baryon and lepton
number transformation, the PQ charges in the M4 model are fixed to be

$\displaystyle\mathcal{X}_{q}$
$\displaystyle=\mbox{diag}\,(0,0,0)\,,\quad\quad\mathcal{X}_{u}=\mbox{diag}\,(s^{2}_{\beta},\,s^{2}_{\beta},\,s^{2}_{\beta})\,,\quad\mathcal{X}_{d}=\mbox{diag}\,(c^{2}_{\beta},\,-s^{2}_{\beta},\,-s^{2}_{\beta})\,,$
$\displaystyle\mathcal{X}_{\ell}$
$\displaystyle=-\mbox{diag}\,(0,1,1)\,,\quad\mathcal{X}_{e}=-\mathcal{X}_{u}\,,$
(2.37)

which imply non-universality in the RH down sector and in the LH lepton
sector.

With the specific choice of the trilinear $H_{2}^{\dagger}H_{1}\phi$ term in
the scalar potential, one has $2N=1$ both in the M1 and M4 models. Flavour
violating axion couplings to fermions (as defined in Eq. (2.8)) arise
exclusively in the LH sector for the M1 model, i.e.

$\displaystyle(C^{V,A}_{u})_{i\neq j}$
$\displaystyle=\pm(V_{u_{L}}^{\dagger}\mathcal{X}_{q}V_{u_{L}})_{i\neq j}\,,$
(2.38) $\displaystyle(C^{V,A}_{d})_{i\neq j}$
$\displaystyle=\pm(V_{d_{L}}^{\dagger}\mathcal{X}_{q}V_{d_{L}})_{i\neq j}\,,$
(2.39) $\displaystyle(C^{V,A}_{e})_{i\neq j}$
$\displaystyle=\pm(V_{e_{L}}^{\dagger}\mathcal{X}_{\ell}V_{e_{L}})_{i\neq
j}\,,$ (2.40)

while for the M4 model flavour violation arises in the RH down-quark sector
and LH charged-lepton sector respectively:

$\displaystyle(C^{V,A}_{u})_{i\neq j}$ $\displaystyle=0\,,$ (2.41)
$\displaystyle(C^{V,A}_{d})_{i\neq j}$
$\displaystyle=(V_{d_{R}}^{\dagger}\mathcal{X}_{d}V_{d_{R}})_{i\neq j}\,,$
(2.42) $\displaystyle(C^{V,A}_{e})_{i\neq j}$
$\displaystyle=\pm(V_{e_{L}}^{\dagger}\mathcal{X}_{\ell}V_{e_{L}})_{i\neq
j}\,.$ (2.43)

Moreover, the combination of flavour violating couplings appearing in $\Delta
F=2$ observables, defined in Eqs. (2.3)–(2.3), take the simpler form (both in
the M1 and M4 models)

$\displaystyle(\epsilon^{d}_{ij})^{2}+(\epsilon^{d}_{ij})^{*\,2}$
$\displaystyle=2\left(\frac{1}{vs_{\beta}}\right)^{2}(m_{d_{i}}^{2}+m_{d_{j}}^{2})(C^{V,A}_{d})^{2}_{ij}\,,$
(2.44) $\displaystyle 2\,\epsilon^{d}_{ij}(\epsilon^{d}_{ij})^{*}$
$\displaystyle=4\left(\frac{1}{vs_{\beta}}\right)^{2}m_{d_{i}}\,m_{d_{j}}(C^{V,A}_{d})^{2}_{ij}\,,$
(2.45)

with the analogous expressions for the up-quark sector obtained by replacing
$d\to u$ and $s_{\beta}\to c_{\beta}$.

##  3 Flavour phenomenology

With the aid of the relations obtained in the previous section we can now
investigate the interplay between axion flavour phenomenology and the
signatures arising from the heavy Higgs radial modes of the PQ-2HDM. We start
by discussing LHC constraints on low-scale, i.e. $\mathcal{O}(1\,\text{TeV})$,
2HDM and then consider all relevant flavour bounds on the non-universal
PQ-2HDM, both in the quark and lepton sectors. This phenomenological analysis
actually applies to a general 2HDM with tree-level flavour changing couplings,
regardless of the presence of the PQ symmetry. All the details regarding the
2HDM potential and its mass spectrum are deferred to Appendix A, while the
analysis of meson mixing and a collection of formulae for LFV processes can be
found in Appendix C and E. Finally, based on the IR/UV flavour connection (cf.
Eqs. (2.22)–(2.24)), we discuss the consequences of an axion discovery in
flavour violating decays on the PQ-2HDM parameter space.

###  3.1 LHC bounds on the 2HDM

The 2HDM scalar sector can be tested at the LHC in two different ways: either
indirectly via the modification of Higgs couplings or by directly hunting for
new scalar resonances. The first type of searches are only sensitive to the
mixing angles in the scalar sector, while direct searches are also sensitive
to the typical mass scale of the enlarged scalar sector.

LHC Higgs data implies that the observed properties of the $h(125)$ resonance
are very close to those predicted by the SM. In particular, observables such
as Higgs couplings to fermions and SM gauge bosons yield already a robust
upper bound on the rotation matrix that diagonalize the neutral scalar sector,
$c_{\alpha-\beta}\lesssim 0.3(0.1)$ [46] for a type-I(II) 2HDM [47].444Recall
that (universal) 2HDM’s with natural flavour conservation are conventionally
classified as $i)$ type-I model if only one of the two scalar doublets couples
(universally) to all fermions or $ii)$ type-II model if up- and down-quarks
sectors couple to different doublets. In particular, the (universal) DFSZ
axion model is compatible only with the type-II structure.

In order to automatically satisfy LHC constraints on Higgs couplings, the so-
called SM alignment limit [48] has been introduced, in which one of the two
Higgs doublets is nearly aligned with the electroweak VEV, thus reproducing
the properties of a SM-like Higgs. There are basically two ways of achieving
alignment. The first one is alignment with decoupling, in which the condition
$c_{\alpha-\beta}\to 0$ is obtained by decoupling all the heavy radial modes,
$m_{H}\simeq m_{A}\simeq m_{H^{+}}\gg m_{h}$ (cf. Eq. (A.31)). In this case we
identify $h$ and $H$ respectively as the SM-like and non-SM-like Higgs
scalars. Alignment without decoupling is obtained, instead, by taking the
scalar potential parameter $\Lambda_{6}\to 0$ which allows for a second light
Higgs in the spectrum with a mass either larger [49] (if $c_{\alpha-\beta}\to
0$) or smaller [50] (if $s_{\alpha-\beta}\to 0$) than the observed scalar
resonance at 125 GeV. The field $H$ then labels the non-SM-Higgs (SM-Higgs) in
the first (latter) eventuality. Notice that in the alignment without
decoupling scenario the mass spectrum of the scalar particles is not
necessarily degenerate and it is determined, in particular, by the parameter
$\Lambda_{5}$ of the scalar potential. Hereafter, we will only consider the
case $m_{H}\gtrsim m_{h}$ (i.e. $c_{\alpha-\beta}\simeq 0$) and therefore $h$
will always denote the SM-like Higgs.

LHC has also greatly contributed to direct searches for an extended scalar
sector, looking for both neutral [51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61,
62, 63, 64, 65, 66, 67, 68, 69] and charged scalar particles [70, 71, 72, 73].
Different mechanisms can be responsible for the production of new scalar
resonances: $i)$ Vector Boson Fusion (VBF) is proportional to
$c_{\alpha-\beta}$ and hence substantially suppressed in the alignment limit,
$ii)$ gluon-gluon Fusion (ggF) depends, instead, parametrically on $t_{\beta}$
and therefore bottom quarks can contribute sizeably to ggF. Direct and
indirect searches, as well as different production mechanisms are hence used
to give bounds on the 2HDM in the $(m_{X},t_{\beta})$ plane (where
$X=H,A,H^{+}$), for a given value of $c_{\alpha-\beta}$. For instance,
searches for scalar particle decays into gauge bosons test the lower parts of
the $(m_{X},t_{\beta})$ plane for type-I and II models [52, 55, 64, 54, 53,
59], while type-II enhanced di-$\tau$ decay set a strong upper bound on
$t_{\beta}$ [51, 60]. For recent analyses on the parameter space of the 2HDM,
see e.g. [74, 75, 76]. In particular, a recent global fit to type-II 2HDM [76]
shows that $m_{H}\gtrsim 500\,$GeV independently of $t_{\beta}$, while the
bound stretches up to $m_{H}\simeq 1\,$TeV for $t_{\beta}\gtrsim 1$.

One can partially use type-I and II 2HDM collider searches to provide
approximate bounds for the non-universal PQ-2HDM models by applying the
present limits on the 33 entries of the Yukawa matrices ($y_{t,b,\tau}$) of
the benchmark M1 and M4 models introduced in Section 2. In this approximation
the type-I and II models searches can be related to the M1 and M4 Yukawa
structures, as summarized in Table 1 for the 3rd generation fermion couplings
to the $H$ scalar field. However, although M1 and M4 models behave only
partially as type-I and II 2HDM, in the following, as a conservative approach,
we will regardless employ bounds from the type-II model, which turn out to be
more stringent.

Model | $y_{t}$ | $y_{b}$ | $y_{\tau}$
---|---|---|---
Type-I | $\dfrac{s_{\alpha}}{s_{\beta}}$ | $\dfrac{s_{\alpha}}{s_{\beta}}$ | $\dfrac{s_{\alpha}}{s_{\beta}}$
Type-II | $\dfrac{s_{\alpha}}{s_{\beta}}$ | $\dfrac{c_{\alpha}}{c_{\beta}}$ | $\dfrac{c_{\alpha}}{c_{\beta}}$
M1 | $\dfrac{s_{\alpha}}{s_{\beta}}$ | $\dfrac{c_{\alpha}}{c_{\beta}}$ | $\dfrac{s_{\alpha}}{s_{\beta}}$
M4 | $\dfrac{c_{\alpha}}{c_{\beta}}$ | $\dfrac{s_{\alpha}}{s_{\beta}}$ | $\dfrac{c_{\alpha}}{c_{\beta}}$
Table 1: Scaling of the $H$ couplings to 3rd generation SM fermions, in the
type-I and II 2HDM compared to the M1 and M4 models. Type-I (II) couplings are
emphasized in blue (red), while grey shaded couplings have no counterpart in
universal 2HDMs. Couplings which are the same in both models are colored in
purple.

###  3.2 Flavour constraints on the non-universal 2HDM

The strongest constraints on the quark flavour violating couplings of the
heavy radial modes of the PQ-2HDM arise from neutral mesons mixing. The
parameter space of a flavourful 2HDM has been previously studied also in
connection with LHC data (see for example [77]). Previous analyses often
relied on Yukawa textures specifically introduced in order to suppress tree-
level flavour changing neutral currents (FCNCs), as e.g. in the case of the
Cheng-Sher ansatz [78]. In this Section, conversely, we do not assume any
particular form for the Yukawa matrices, and perform, instead, an analysis of
meson mixing constraints on the parameter space of a general 2HDM with tree-
level flavour-changing couplings. Our numerical analysis, including QCD
running effects, follows the one in Ref. [77] (see also [79, 80, 81]).
Additionally, we take into account (whenever available) the theoretical error
on the SM predictions for the meson mixing parameters. In the case of $B$
systems theoretical uncertainties are at the level of $\sim 5\%$ [82], while
Kaon physics has not yet reached this level of precision and lattice
calculations provide $K$-form factors with an $\mathcal{O}(30\%)$ error, see
for example Refs. [83, 84].

The details of the analysis are reported in Appendix C. Here, we summarize the
main results in the alignment scenario, which is necessary in order to satisfy
LHC constraints on $c_{\alpha-\beta}$. The bounds from meson mixing can be
recast in this limit as

$\displaystyle\left(\frac{[4.1,\,3.6]\,\text{TeV}}{s_{2\beta}\,m_{H}}\right)^{2}\frac{m_{H}^{2}+m_{A}^{2}}{2m_{A}^{2}}\left|\,\frac{\epsilon_{12}^{d}\epsilon^{d\,*}_{21}}{y_{s}^{2}\,\lambda^{2}}+0.12\frac{m_{H}^{2}-m_{A}^{2}}{m_{H}^{2}+m_{A}^{2}}\frac{(\epsilon_{12}^{d})^{2}+(\epsilon^{d\,*}_{21})^{2}}{y_{s}^{2}\,\lambda^{2}}\right|<1\,,$
(3.1)
$\displaystyle\left(\frac{[1.4,\,1.7]\,\text{TeV}}{s_{2\beta}\,m_{H}}\right)^{2}\frac{m_{H}^{2}+m_{A}^{2}}{2m_{A}^{2}}\left|\,\frac{\epsilon_{13}^{d}\epsilon^{d\,*}_{31}}{y_{b}^{2}\lambda^{6}}+0.17\frac{m_{H}^{2}-m_{A}^{2}}{m_{H}^{2}+m_{A}^{2}}\frac{(\epsilon_{13}^{d})^{2}+(\epsilon^{d\,*}_{31})^{2}}{y_{b}^{2}\lambda^{6}}\right|<1\,,$
(3.2)
$\displaystyle\left(\frac{[1.5,\,1.6]\,\text{TeV}}{s_{2\beta}\,m_{H}}\right)^{2}\frac{m_{H}^{2}+m_{A}^{2}}{2m_{A}^{2}}\left|\,\frac{\epsilon_{23}^{d}\epsilon^{d\,*}_{32}}{y_{b}^{2}\lambda^{4}}+0.17\frac{m_{H}^{2}-m_{A}^{2}}{m_{H}^{2}+m_{A}^{2}}\frac{(\epsilon_{23}^{d})^{2}+(\epsilon^{d\,*}_{32})^{2}}{y_{b}^{2}\lambda^{4}}\right|<1\,,$
(3.3)
$\displaystyle\left(\frac{23\,\text{TeV}}{m_{H}}\right)^{2}\frac{m_{H}^{2}+m_{A}^{2}}{2m_{A}^{2}}\left|\,\frac{\epsilon_{12}^{u}\epsilon^{u\,*}_{21}}{y_{c}^{2}\lambda^{2}}+0.15\frac{m_{H}^{2}-m_{A}^{2}}{m_{H}^{2}+m_{A}^{2}}\frac{(\epsilon_{12}^{u})^{2}+(\epsilon^{u\,*}_{21})^{2}}{y_{c}^{2}\lambda^{2}}\right|<1\,,$
(3.4)

where we have equated the contribution of the Wilson coefficients of the
releavant 4-fermion operators to the mass scale written in the bracket and
rescaled the UV sources of flavour violation $(\epsilon^{u,d})$ with suitable
powers of Yukawas,
$y_{d_{i}(u_{i})}=\frac{\sqrt{2}m_{i}}{vs_{\beta}(c_{\beta})}$, and Cabibbo
angle, $\lambda$. This rescaling corresponds to a CKM-like structure for the
LH rotation matrices, $V_{u_{L}}\sim V_{d_{L}}\sim V_{\rm CKM}$, and
$V_{u_{R}}\sim V_{d_{R}}\sim\mathbbm{1}$ for the RH ones. For example, in Eq.
(3.4) one has
$\epsilon_{12}^{u}=(V_{u_{L}}^{\dagger}Y_{1}^{u}V_{u_{R}})_{12}\sim(V_{u_{L}})_{12}^{*}(Y^{u}_{1})_{22}(V_{u_{R}})_{22}\sim\lambda\frac{\sqrt{2}\,m_{c}}{v\,c_{\beta}}$.
To obtain the numerical values reported in Eqs. (3.1)–(3.4) we have used the
PDG values for $m_{i}(\mu)$, with $\mu=\\{2,\,2.8,\,4.6\\}$ GeV respectively
for $K$, $D$ and $B$ meson oscillations and $\lambda=0.2$. For down-type
flavour transitions, where reliable SM theoretical predictions are available,
the number quoted in the left (right) of the square bracket corresponds to a
negative (positive) contribution from new physics. For
$D^{0}-\overline{D}^{0}$ oscillation we have required, instead, that the new
physics contribution would not overshoot the available experimental bound.

Note that one of the mass terms (e.g. $m_{A}$) can be traded, via Eq. (A.32),
for the scalar potential parameter $\Lambda_{5}$, subject to the
perturbativity constraint $\Lambda_{5}\lesssim 4\pi$. The alignment scenario
with decoupling is readily recovered by taking the $m_{H}^{2}\to m^{2}_{A}$
limit in Eqs. (3.1)–(3.4).

Employing a CKM-like structure for the fermion rotation matrices, we have
studied the implication of the meson mixing bounds on the benchmark models M1
and M4 (see Appendix D for details). In summary, the CKM suppression is
sufficient in order to make these models accessible at the LHC. For M1 the
bound from meson mixing is $m_{H}\gtrsim 200$ GeV, while for M4 is
$m_{H}\gtrsim 1$ TeV. This shows that, despite having FCNCs at tree level,
certain PQ charge arrangements and flavour ansatz can relax the bounds on the
mass scale of the (non-universal) 2HDM, similarly to the case of BGL models
[85] (see also [86]). In Appendix D we have also considered the meson mixing
bounds induced by the SM-like Higgs in the $(c_{\alpha-\beta},t_{\beta})$
plane, in the case in which the heavy Higgs ($H$) is not completely decoupled.
These bounds are found to constrain similar regions of parameter space, as
those arising from the SM-like Higgs LFV decays, considered in Ref. [37].

Regarding instead the leptonic sector of the (non-universal) 2HDM with tree-
level FCNCs, we can mostly rely on previous analyses, such as those in Refs.
[87, 88]. Indeed, the bounds from $l_{i}\to l_{j}\gamma$ [89, 90] have been
only marginally improved [91, 92]; while leptonic decays of the type $l_{i}\to
3l_{j}$ or $l_{i}\to l_{j}l_{k}l_{k}$ are still based on the SINDRUM [93] and
Belle [94] experiments, also taken into account in [88]. On the other hand,
the limits on Higgs LFV decays have been considerably improved in the recent
years by the LHC [95, 96]. In particular, very strong limits on the parameter
space of the (non-universal) 2HDM arise from the $h\to\tau\mu$ and $h\to\tau
e$ processes, while for the case of $\mu$-$e$ transitions $\mu\to e\gamma$
still provides the most relevant constraint. The Higgs LFV branching ratios,
in the notation of Eq. (2.12), read

$\text{BR}(h\to
l_{i}l_{j})=\frac{m_{h}}{16\pi\Gamma_{h}}(|C^{h_{e}}_{ij}|^{2}+|C^{h_{e}}_{ji}|^{2})\,,$
(3.5)

where $\Gamma_{h}$ could receive corrections from new physics (cf. Eq. (E.9)).
By means of Eq. (2.13) we can rewrite the $2\sigma$ bound from $h\to\tau\mu/e$
directly on the $\epsilon^{e}$ coefficients as

$\frac{4.1\,\text{MeV}}{\Gamma_{h}}\frac{|\epsilon^{e}_{i3}|^{2}+|\epsilon^{e}_{3i}|^{2}}{y_{\tau}^{2}}\left(\frac{c_{\alpha-\beta}}{s_{\beta}c_{\beta}}\right)^{2}\lesssim
0.066\,,$ (3.6)

where $i=1,\,2$ (since at $2\sigma$ both the $\mu$ and $e$ channels yield a
similar bound) and $y_{\tau}=\sqrt{2}m_{\tau}/vs_{\beta}$. The Higgs branching
ratio has been normalized to the SM model one, $\Gamma_{h}^{\rm SM}\simeq
4.1\,$MeV.

###  3.3 IR/UV quark flavour connection

As seen in Section 2.3 there is a direct relation between the IR sources of
flavour violation ($C_{f}^{V,A}$) and the UV ones ($\epsilon^{f}_{V,A}$). This
is particularly useful in the case of a positive detection of a new physics
effect, since once the strength of the flavour violation is known in one of
the two sectors, it can be directly translated into the other one. The second
important point is the generality of these correlations. While the flavour
violating axion signals depend on the specific model via the PQ charges and
the fermion mixing matrices $V_{f_{L/R}}$, the correlations derived in Section
2.3 are model-independent, in the sense that they directly relate the flavour
coefficients entering the physical observables.

In the following, we highlight the IR/UV connection in a class of non-
universal PQ models with a 2+1 structure of the PQ charges, like the M1 and M4
models introduced in Section 2.4, exemplifications of scenarios in which
flavour violation in the axion couplings arises either in the LH or RH sector.
The generalization to the most general setup with both LH and RH flavour
violation is then straightforward. Through Eqs. (2.3)–(2.3) we can directly
translate the information on the PQ symmetry breaking scale $f_{a}$, appearing
explicitly in the IR axion couplings to fermions into limits on the heavy 2HDM
scale $m_{H}$ (and $t_{\beta}$) by relating flavour violating axion
observables with meson mixing ones:

$\displaystyle\left(\frac{f_{a}}{10^{11}\,\text{GeV}}\right)^{2}\left(\frac{1\,\text{TeV}}{s_{2\beta}\,m_{H}}\right)^{2}\left(\frac{\text{BR}(K\to\pi
a)}{7.3\cdot 10^{-11}}\right)=\frac{2\,|M^{\text{NP}}_{12}|}{3.5\cdot
10^{-15}\,\text{GeV}}\,,$ (3.7) $\displaystyle\left(\frac{f_{a}}{2.8\cdot
10^{7}\,\text{GeV}}\right)^{2}\left(\frac{1\,\text{TeV}}{s_{2\beta}\,m_{H}}\right)^{2}\left(\frac{\text{BR}(B\to\pi
a)}{2.3\cdot 10^{-5}}\right)=\frac{2\,|M^{\text{NP}}_{12}|}{3.3\cdot
10^{-13}\,\text{GeV}}\,,$ (3.8) $\displaystyle\left(\frac{f_{a}}{8.8\cdot
10^{7}\,\text{GeV}}\right)^{2}\left(\frac{1\,\text{TeV}}{s_{2\beta}\,m_{H}}\right)^{2}\left(\frac{\text{BR}(B\to
Ka)}{7.1\cdot 10^{-6}}\right)=\frac{2\,|M^{\text{NP}}_{12}|}{1.2\cdot
10^{-11}\,\text{GeV}}\,,$ (3.9) $\displaystyle\left(\frac{f_{a}}{1.2\cdot
10^{7}\,\text{GeV}}\right)^{2}\left(\frac{1\,\text{TeV}}{s_{2\beta}\,m_{H}}\right)^{2}\left(\frac{\text{BR}(D\to\pi
a)}{8\cdot 10^{-6}}\right)=\frac{2\,|M^{\text{NP}}_{12}|}{6.7\cdot
10^{-15}\,\text{GeV}}\,.$ (3.10)

These expressions have been written, for simplicity, in the alignment limit
(i.e. $c_{\alpha-\beta}\simeq 0$) assuming $m_{H}^{2}\simeq m_{A}^{2}$.
Flavour violating axion observables and meson mixing parameters have been
normalized to current experimental limits. In particular, the 90% CL limits on
the branching ratios involving the axion filed are taken from Ref. [28].

Figure 1: IR/UV quark flavour connection in the $(m_{a},s_{2\beta}m_{H})$
plane, assuming an axion discovery in a given flavour transition, with solid
(dotted) lines corresponding to current bounds (future sensitivities) on meson
decays. The upper-right region is allowed by current bounds on $K^{0}$ (upper
plot) and $B^{0}_{d}$ (lower plot) meson mixing. LHC bounds (set approximately
at 1 TeV) are displayed in grey, while the brown regions correspond to
astrophysical bounds on the axion decay constant, between $f_{a}\sim 3\times
10^{7}\,$GeV (astrophobic scenario) and $f_{a}\sim 2\times 10^{9}\,$GeV, an
axion with a sizeable coupling to electrons ($C_{e}=1$).

Figure 2: IR/UV quark flavour connection in the $(m_{a},s_{2\beta}m_{H})$
plane, assuming an axion discovery in a given flavour transition, with solid
(dotted) lines corresponding to current bounds (future sensitivities) on meson
decays. The upper-right region is allowed by current bounds on $B^{0}_{s}$
(upper plot) and $D^{0}$ (lower plot) meson mixing. LHC bounds (set
approximately at 1 TeV) are displayed in grey, while the brown regions
correspond to astrophysical bounds on the axion decay constant, between
$f_{a}\sim 3\times 10^{7}\,$GeV (astrophobic scenario) and $f_{a}\sim 2\times
10^{9}\,$GeV, an axion with a sizeable coupling to electrons ($C_{e}=1$).

The connections established in Eqs. (3.7)–(3.10) fully prove their relevance
once an evidence of deviation from the SM will emerge. Assuming, for example,
an axion discovery in a given flavour violating meson decay, it will be then
straightforward to relate it to a possible deviation from SM predictions in
meson mixing observables, induced by the heavy radial modes of the
PQ-2HDM.555Note that meson mixing receives also an IR contribution from the
tree-level exchange of the axion field. However, this contribution is
subleading compared to the UV-induced one, as long as $f_{a}\gg m_{H}$. Fixing
the branching ratio of the flavour violating axion observable in a region
between the current experimental bound and the future expected sensitivity
(the latter is also taken from Ref. [28]), we can translate the current bound
imposed by meson mixing into a relation of the type
$f_{a}<\\#_{f_{a}}s_{2\beta}m_{H}$, with $\\#_{f_{a}}$ a number which can be
read directly from Eqs. (3.7)–(3.10) (or equivalently
$m_{a}>\\#_{m_{a}}s_{2\beta}m_{H}$, using the QCD axion mass relation,
$m_{a}\simeq 5.7\times(10^{6}\ \text{GeV}/f_{a})$ eV). This region corresponds
the white upper-right half-plane in the plots of Figs. 1–2, where we also
superimpose LHC constraints and astrophysical bounds on the axion decay
constant. In the latter case we assume two scenarios corresponding to a
suppressed axion coupling to nucleons and electrons (astrophobic scenario
[21]) and $C_{e}=1$ (corresponding to a sizeable value of the axion coupling
to electrons).

If the axion is discovered in the near future (i.e. in the region between
current bounds and future sensitivities) in $B$\- or $D$-meson transitions, it
is evident from Figs. 1–2 that the interplay of meson mixing and astrophysical
bounds disfavour the possibility of probing the radial modes of the PQ-2HDM at
the LHC (unless the axion is to some extent astrophobic). On the other hand, a
discovery of the axion in $K\to\pi a$ leaves open plenty of parameter space
for a light PQ-2HDM in the LHC range, since astrophysical constraints are not
relevant for large $f_{a}$.

A similar discussion applies if a SM deviation will be detected in meson
mixing. In such a case Figs. 1–2 would be filled in the upper-right plane, and
one could draw well-definite conclusions regarding the possibility to observe
an axion within the reach of current or future meson decay experiments. This
scenario is however more difficult to be realized in practice, since it relies
on a precise SM computation for the meson mixing parameters.

###  3.4 IR/UV lepton flavour connection

In the charged-lepton sector it is possible to exploit Eq. (2.24) in order to
connect LFV axion processes like $l_{i}\to l_{j}\,a$ with other LFV
observables such as $l_{i}\to 3l_{j}$, $l_{i}\to l_{j}l_{k}l_{k}$ and
$l_{i}\to l_{j}\gamma$. However, as shown in Appendix E, the relation between
this two classes of LFV observables is not straightforward as in the quark
case considered in Section 3.3, because of the appearance of several couplings
which break the 1-to-1 flavour correspondence. Hence, a global analysis would
be necessary in this case.

Another possibility is to consider LFV Higgs decays $h\to l_{i}l_{j}$, which
actually turns out to be in direct connection with LFV axion processes, thus
providing a useful simplification. Moreover, there is a recent hint from ATLAS
[96] about a possible deviation from the SM in the Higgs branching ratios
$\text{BR}(h\to\tau e)=0.09\pm 0.06\,\%$ and
$\text{Br}(h\to\tau\mu)=0.11^{+0.05}_{-0.04}\,\%$. These provide indeed a
useful benchmark in order to highlight the IR/UV lepton flavour connection.
Note, also, that deviations at the 0.1$\%$ level in $\text{BR}(h\to\tau
e/\mu)$ within a non-universal 2HDM are still allowed by other LFV observables
(see e.g. [97]). Employing Eq. (3.5) and the IR/UV connection formulae in
Section 2.3, we obtain the following expression for the LFV Higgs decays in
terms of axion couplings

$\text{BR}(h\to
l_{i}l_{j})\simeq\frac{m_{h}}{16\pi\Gamma_{h}}\left(\frac{c_{\alpha-\beta}}{s_{\beta}c_{\beta}}\right)^{2}\frac{m_{l_{i}}^{2}}{v^{2}}|(C^{L,R}_{e})_{ij}|^{2}\,,$
(3.11)

which applies to models with 2+1 PQ charges (the extension to general non-
universal DFSZ models, with both LH and RH axion flavour violation, is
provided in Eq. (E)). Given a LFV process involving the axion field [29]

$\text{BR}(l_{i}\to
l_{j}\,a)\simeq\frac{m_{l_{i}}^{3}}{16\pi\Gamma_{l_{i}}}\frac{|(C^{L,R}_{e})_{ij}|^{2}}{2f_{a}^{2}}\,,$
(3.12)

we can provide a direct relation for the two branching ratios, that reads

$\text{BR}(h\to l_{i}l_{j})\simeq\text{BR}(l_{i}\to
a\,l_{j})\frac{2m_{h}}{m_{l_{i}}}\frac{\Gamma_{l_{i}}}{\Gamma_{h}}\frac{f_{a}^{2}}{v^{2}}\left(\frac{c_{\alpha-\beta}}{c_{\beta}s_{\beta}}\right)^{2}\,.$
(3.13)

Employing the data in Table 2, we display in Fig. 3 the correlation provided
by Eq. (3.13) in the $(m_{a},c_{\alpha-\beta}/s_{2\beta})$ plane. Assuming as
a benchmark the hint from ATLAS on $h\to\tau\mu/e$ decays, we impose the
bounds and future sensitivities from the axion processes $\tau\to ea$ (top
panel) and $\tau\to\mu a$ (bottom panel), as well as reference astrophysical
limits on the axion mass and LHC constraints on $c_{\alpha-\beta}$ for typical
model parameters. The combination of all bounds corners a region of the axion
mass parameter space that will be targeted by future axion searches (such as
those with helioscopes) and that is characterized by a remarkable
complementarity between different types of observables, ranging from flavour
physics and LHC direct searches to axion astrophysics.

$i,\,j$ | BR($h\to l_{i}l_{j}$) | Current BR($l_{i}\to a\,l_{j}$) | Expected BR($l_{i}\to al_{j}$)
---|---|---|---
$\tau,\,e$ |  $0.09\pm 0.06\,\%$ [96] |  $9.4\times 10^{-4}$ [98] |  $8.3\times 10^{-6}$ [29]
$\tau,\,\mu$ |  $0.11^{+0.05}_{-0.04}\,\%$ [96] |  $5.9\times 10^{-4}$ [98] |  $2.0\times 10^{-5}$ [29]
Table 2: Input values used in the LFV analysis in Fig. 3. Lepton decays
correspond to $95\%$ CL bounds, while the Higgs branching ratios are
measurements at $1\sigma$.

Figure 3: IR/UV lepton flavour connection in the
$(m_{a},c_{\alpha-\beta}/s_{2\beta})$ plane, assuming the hint on LFV Higgs
decays $h\to\tau\mu/e$ and imposing the bound (colored purple region being
excluded) from $\tau\to ea$ (top panel) and $\tau\to\mu a$ (bottom panel).
Solid (dotted) lines correspond to current bounds (future sensitivities) from
$\tau$ decays, while dashed lines refer to the 1$\sigma$ region for the
$h\to\tau\mu/a$ hint. For input data see Table 2. LHC bounds on
$c_{\alpha-\beta}$ (for $t_{\beta}=1$) are displayed in grey as a reference in
the case of the type-I and -II 2HDM, while the brown regions correspond to
astrophysical bounds on the axion decay constant, between $f_{a}\sim 3\times
10^{7}\,$GeV (astrophobic scenario) and $f_{a}\sim 2\times 10^{9}\,$GeV, for
an axion with a sizeable coupling to electrons ($C_{e}=1$).

##  4 Conclusions

In this paper we have studied the flavour phenomenology of non-universal DFSZ
axion models, focusing on the interplay between IR and UV sources of flavour
violation. These consist respectively in the flavour violating axion couplings
to fermions, whose chiral counterparts are denoted as $C_{f}^{L,R}$ (see Eq.
(2.8)) and the flavour violating couplings of the heavy radial modes of the
PQ-2HDM, which are expressed in terms of the $\epsilon^{f}$ parameters
(defined in Eq. (2.17)). The key relations, providing a direct link between
the $\epsilon^{f}$ and $C_{f}^{L,R}$ coefficients, are given by Eqs.
(2.22)–(2.24). These remarkably simple equations, that lock the flavour
violating pattern of the axion field to that of the heavy radial modes, hold
in general for any (non-universal) DFSZ models, with a renormalizable 2HDM
embedded. In fact, in order to derive Eqs. (2.22)–(2.24) we have exploited the
${\rm U}(1)_{\rm PQ}$ invariance of the PQ-2HDM Yukawa sector made explicit in
Eqs. (2.2)–(2.4).

With the help of Eqs. (2.22)–(2.24) we have hence investigated the flavour
phenomenology both in the quark and charged-lepton sectors, also in connection
with LHC limits on the PQ-2HDM and astrophysical constraints on axion
couplings. In the quark sector we focused on meson mixing observables, which
yield the strongest constraints on the flavour violating couplings of the non-
universal PQ-2HDM radial modes. We have also derived a set of useful formulae
(see Eqs. (3.1)–(3.4)), from which present bounds on meson mixing mass
differences can be translated to tree-level flavour violating couplings of a
general 2HDM in the alignment limit, regardless of the presence of the PQ
symmetry.

In order to further exemplify the IR/UV flavour connection we have considered
two explicit examples of 2+1 models, labelled M1 and M4, which are prototypes
of models with quark flavour violation in axion couplings, respectively in the
LH and in the RH down-quark sector. Assuming an axion discovery in flavour
violating meson decays, and imposing in turn the bounds from meson mixing
experiments we have obatined the results shown in Figs. 1–2. The main message
that can be learnt with such exercise is that if an axion is discovered in the
near future in $B$\- or $D$-meson transitions, then the interplay of meson
mixing and astrophysical bounds disfavours the possibility of probing the
radial modes of the PQ-2HDM at the LHC. On the other hand, an axion discovery
in $K\to\pi a$ leaves open plenty of parameter space for a light PQ-2HDM in
the LHC range.

A similar IR/UV connection can be envisaged as well for the charged-lepton
sector. However, LFV observables, such as for example $l_{i}\to l_{j}\gamma$,
turn out to depend on several axion couplings, making a direct connection with
LFV axion processes not straightforward. To this end, we have considered
instead SM-like Higgs decays $h\to l_{i}l_{j}$, which actually provide a
1-to-1 flavour correspondence with LFV axion processes such as $l_{i}\to
l_{j}a$. Given the fact that ATLAS [96] currently hints to a mild deviation
from the SM in the Higgs decays $h\to\tau e/\mu$, we have taken this as a
benchmark scenario and imposed the bounds from axion LFV processes, exploiting
again the relation between IR and UV sources of flavour violation. The results
of this analysis are displayed in Fig. 3, which shows again a strong
complementarity between flavour observables, LHC direct searches and axion
astrophysics.

## Acknowledgements

We thank Mario Fernández Navarro, Belén Gavela, Luca Merlo and Robert Ziegler
for useful discussions and comments. This work received funding from the
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement N∘ 860881-HIDDeN, grant agreement N∘
101086085–ASYMMETRY and by the INFN Iniziative Specifica APINE. The work of
LDL is also supported by the project “CPV-Axion” under the Supporting TAlent
in ReSearch@University of Padova (STARS@UNIPD).

##  Appendix A PQ-2HDM scalar potential

In this Appendix we provide some details on the scalar potential, with a
special focus on the PQ-2HDM spectrum in the PQ broken phase (see also [99,
100]). The DFSZ scalar potential, in the presence of the cubic invariant
$H_{2}^{\dagger}H_{1}\phi$, can be written as

$\displaystyle V\left(H_{1},H_{2},\phi\right)$
$\displaystyle=\mu_{1}^{2}|H_{1}|^{2}+\mu_{2}^{2}|H_{2}|^{2}+\frac{\lambda_{1}}{2}|H_{1}|^{4}+\frac{\lambda_{2}}{2}|H_{2}|^{4}+\lambda_{3}|H_{1}|^{2}|H_{2}|^{2}$
$\displaystyle+\lambda_{4}(H_{1}^{\dagger}H_{2})(H_{2}^{\dagger}H_{1})+\frac{\lambda^{\phi}_{1}}{2}|\phi|^{2}|H_{1}|^{2}+\frac{\lambda^{\phi}_{2}}{2}|\phi|^{2}|H_{2}|^{2}+\lambda^{\phi}_{3}\left(|\phi|^{2}-\frac{v^{2}_{\phi}}{2}\right)^{2}$
$\displaystyle-\left(\mu_{\phi}H^{\dagger}_{2}H_{1}\phi+\text{h.c.}\right)\,,$
(A.1)

where the Higgs doublets are decomposed as

$H_{1}=\begin{pmatrix}\frac{1}{\sqrt{2}}(v_{1}+h^{0}_{1}+i\eta_{1})\\\
h^{-}_{1}\end{pmatrix}\,,\quad
H_{2}=\begin{pmatrix}\frac{1}{\sqrt{2}}(v_{2}+h^{0}_{2}+i\eta_{2})\\\
h^{-}_{2}\end{pmatrix}\,.\quad$ (A.2)

In the PQ-broken phase, with $\langle\phi\rangle=v_{\phi}/\sqrt{2}$, we can
integrate out the radial mode of the $\phi$ field and obtain the leading-order
PQ-2HDM potential

$\displaystyle V\left(H_{1},H_{2}\right)$
$\displaystyle=m_{11}^{2}|H_{1}|^{2}+m_{22}^{2}|H_{2}|^{2}-\left(m_{12}H^{\dagger}_{1}H_{2}+\text{h.c.}\right)$
$\displaystyle+\frac{\lambda_{1}}{2}|H_{1}|^{4}+\frac{\lambda_{2}}{2}|H_{2}|^{4}+\lambda_{3}|H_{1}|^{2}|H_{2}|^{2}+\lambda_{4}(H_{1}^{\dagger}H_{2})(H_{2}^{\dagger}H_{1})\,,$
(A.3)

where

$m_{11}^{2}=\mu^{2}_{1}+\frac{\lambda_{1}^{\phi}}{4}v^{2}_{\phi}\,,\qquad
m_{22}^{2}=\mu^{2}_{2}+\frac{\lambda_{2}^{\phi}}{4}v^{2}_{\phi}\,,\qquad
m_{12}^{2}=\mu^{*}_{\phi}\frac{v_{\phi}}{\sqrt{2}}\,.$ (A.4)

To implement the $v_{\phi}\gg v$ hierarchy, we impose the “ultra-weak” scaling
$\lambda_{1,2}^{\phi}\sim v^{2}/v^{2}_{\phi}$ and $\mu_{\phi}\sim
v^{2}/v_{\phi}$, which also guarantees that the hierarchy between the PQ and
the electroweak breaking scales remains radiatively stable thanks to the
emergence of an extended Poincaré symmetry of the action [101, 102]. This is
reflected by the fact that $\lambda_{1,2}^{\phi}\to 0$ and $\mu_{\phi}\to 0$
correspond to fixed points of the renormalization group evolution (see e.g.
App. B in [103]).

To match with the more standard 2HDM notation (see e.g. [48]) we consider the
rotated doublets

$\begin{pmatrix}\Phi_{1}\\\ \Phi_{2}\\\
\end{pmatrix}=\begin{pmatrix}c_{\beta}&s_{\beta}\\\
-s_{\beta}&c_{\beta}\end{pmatrix}\begin{pmatrix}\tilde{H}_{1}\\\
\tilde{H}_{2}\\\ \end{pmatrix}\,,$ (A.5)

so that $\Phi_{1}$ picks up the whole electroweak VEV

$\Phi_{1}=\begin{pmatrix}G^{+}\\\
-\frac{1}{\sqrt{2}}(v+c_{\beta}h^{0}_{1}+s_{\beta}h^{0}_{2}-iG^{0})\end{pmatrix}\,,\quad\Phi_{2}=\begin{pmatrix}H^{+}\\\
-\frac{1}{\sqrt{2}}(-s_{\beta}h^{0}_{1}+c_{\beta}h^{0}_{2}-iA)\end{pmatrix}\,.$
(A.6)

Here we have defined the would-be Goldstone bosons

$G^{+}=c_{\beta}h^{+}_{1}+s_{\beta}h^{+}_{2}\,,\quad
G^{0}=c_{\beta}\eta_{1}+s_{\beta}\eta_{2}\,,$ (A.7)

and the scalar radial modes

$H^{+}=-s_{\beta}h^{+}_{1}+c_{\beta}h^{+}_{2}\,,\quad
A^{0}=-s_{\beta}\eta_{1}+c_{\beta}\eta_{2}\,.$ (A.8)

In the rotated basis the scalar potential reads

$\displaystyle V\left(\Phi_{1},\,\Phi_{2}\right)$
$\displaystyle=M_{11}^{2}\,\Phi_{1}^{\dagger}\Phi_{1}+M_{22}^{2}\,\Phi_{2}^{\dagger}\Phi_{2}-\left(M_{12}\Phi_{1}^{\dagger}\Phi_{2}+\text{h.c.}\right)$
$\displaystyle+\frac{\Lambda_{1}}{2}(\Phi_{1}^{\dagger}\Phi_{1})^{2}+\frac{\Lambda_{2}}{2}(\Phi_{2}^{\dagger}\Phi_{2})^{2}+\Lambda_{3}(\Phi_{1}^{\dagger}\Phi_{1})(\Phi_{2}^{\dagger}\Phi_{2})+\Lambda_{4}(\Phi_{1}^{\dagger}\Phi_{2})(\Phi_{2}^{\dagger}\Phi_{1})$
$\displaystyle+\left\\{\frac{1}{2}\Lambda_{5}(\Phi_{1}^{\dagger}\Phi_{2})^{2}+\left[\Lambda_{6}(\Phi_{1}^{\dagger}\Phi_{1})+\Lambda_{7}(\Phi_{2}^{\dagger}\Phi_{2})\right]\Phi_{1}^{\dagger}\Phi_{2}+\text{h.c.}\right\\}\,,$
(A.9)

with the new parameters given by

$\displaystyle M_{11}^{2}$
$\displaystyle=m_{11}^{2}c_{\beta}^{2}+m_{22}^{2}s_{\beta}^{2}-\mbox{Re}\,{m_{12}^{2}}\,s_{2\beta}\,,$
(A.10) $\displaystyle M_{22}^{2}$
$\displaystyle=m_{11}^{2}s_{\beta}^{2}+m_{22}^{2}c_{\beta}^{2}-\mbox{Re}\,{m_{12}^{2}}\,s_{2\beta}\,,$
(A.11) $\displaystyle M_{12}^{2}$
$\displaystyle=\frac{1}{2}\left(m_{11}^{2}-m_{22}^{2}\right)s_{2\beta}+\mbox{Re}\,{m_{12}^{2}}\,c_{2\beta}+i\,\mbox{Im}\,{m_{12}^{2}}\,,$
(A.12) $\displaystyle\Lambda_{1}$
$\displaystyle=\lambda_{1}c^{4}_{\beta}+\lambda_{2}s^{4}_{\beta}+\frac{1}{2}(\lambda_{3}+\lambda_{4})s^{2}_{2\beta}\,,$
(A.13) $\displaystyle\Lambda_{2}$
$\displaystyle=\lambda_{1}s^{4}_{\beta}+\lambda_{2}c^{4}_{\beta}+\frac{1}{2}(\lambda_{3}+\lambda_{4})s^{2}_{2\beta}\,,$
(A.14) $\displaystyle\Lambda_{3}$
$\displaystyle=\frac{1}{4}s^{2}_{2\beta}\left[\lambda_{1}+\lambda_{2}-2(\lambda_{3}+\lambda_{4})\right]+\lambda_{3}\,,$
(A.15) $\displaystyle\Lambda_{4}$
$\displaystyle=\frac{1}{4}s^{2}_{2\beta}\left[\lambda_{1}+\lambda_{2}-2(\lambda_{3}+\lambda_{4})\right]+\lambda_{4}\,,$
(A.16) $\displaystyle\Lambda_{5}$
$\displaystyle=\frac{1}{4}s^{2}_{2\beta}\left[\lambda_{1}+\lambda_{2}-2(\lambda_{3}+\lambda_{4})\right]\,,$
(A.17) $\displaystyle\Lambda_{6}$
$\displaystyle=-\frac{1}{2}s_{2\beta}\left[\lambda_{1}c_{\beta}^{2}-\lambda_{2}s_{\beta}^{2}-c_{2\beta}(\lambda_{3}+\lambda_{4})\right]\,,$
(A.18) $\displaystyle\Lambda_{7}$
$\displaystyle=-\frac{1}{2}s_{2\beta}\left[\lambda_{1}s_{\beta}^{2}-\lambda_{2}c_{\beta}^{2}+c_{2\beta}(\lambda_{3}+\lambda_{4})\right]\,.$
(A.19)

The stationary conditions read

$\displaystyle\left.\frac{\partial
V}{\partial\Phi_{1}}\right|_{\langle\Phi_{2}\rangle=0,\,\langle\Phi_{1}\rangle=\frac{v}{\sqrt{2}}}=0\,,\quad$
$\displaystyle\Longrightarrow\quad
M_{11}^{2}\frac{v}{\sqrt{2}}+\Lambda_{1}\frac{v^{3}}{2\sqrt{2}}=0\,,$ (A.20)
$\displaystyle\left.\frac{\partial
V}{\partial\Phi_{2}}\right|_{\langle\Phi_{2}\rangle=0,\,\langle\Phi_{1}\rangle=\frac{v}{\sqrt{2}}}=0\,,\quad$
$\displaystyle\Longrightarrow\quad-
M_{12}^{2}\frac{v}{\sqrt{2}}+\Lambda_{6}\frac{v^{3}}{2\sqrt{2}}=0\,,$ (A.21)

from which we can substitute $M^{2}_{11}$ and $M^{2}_{12}$ into the scalar
potential. The squared mass terms of the CP-odd and charged scalars can be
read directly from the Lagrangian

$\displaystyle m_{H^{\pm}}^{2}$
$\displaystyle=M^{2}_{22}+\frac{1}{2}v^{2}\Lambda_{3}\,,$ (A.22)
$\displaystyle m_{A}^{2}$
$\displaystyle=m^{2}_{H^{\pm}}+\frac{1}{2}v^{2}(\Lambda_{4}-\Lambda_{5})=M^{2}_{22}+\frac{1}{2}v^{2}(\Lambda_{3}+\Lambda_{4}-\Lambda_{5})\,,$
(A.23)

while to obtain those of the CP-even Higgses we perform the rotation

$\begin{pmatrix}h^{0}_{2}\\\ h^{0}_{1}\\\
\end{pmatrix}=\begin{pmatrix}c_{\alpha}&s_{\alpha}\\\
-s_{\alpha}&c_{\alpha}\end{pmatrix}\begin{pmatrix}h\\\ H\\\ \end{pmatrix}\,.$
(A.24)

This is equivalent to go into the basis where only one Higgs gets the whole
electroweak VEV

$\begin{pmatrix}m_{H}^{2}&&0\\\
0&&m_{h}^{2}\end{pmatrix}=\begin{pmatrix}c_{\alpha-\beta}&&s_{\alpha-\beta}\\\
-s_{\alpha-\beta}&&c_{\alpha-\beta}\end{pmatrix}\begin{pmatrix}\Lambda_{1}v^{2}&&\Lambda_{6}v^{2}\\\
\Lambda_{6}v^{2}&&m_{A}^{2}+\Lambda_{5}v^{2}\end{pmatrix}\begin{pmatrix}c_{\alpha-\beta}&&-s_{\alpha-\beta}\\\
s_{\alpha-\beta}&&c_{\alpha-\beta}\end{pmatrix}\,,$ (A.25)

with eigenvalues

$m_{h,\,H}^{2}=\frac{1}{2}\left(m_{A}^{2}+v^{2}(\Lambda_{1}+\Lambda_{5})\mp\sqrt{\left[m_{A}^{2}+(\Lambda_{5}-\Lambda_{1})v^{2}\right]^{2}+4v^{4}\Lambda_{6}^{2}}\right)\,,$
(A.26)

and mixing angle given by [48]

$s_{2(\alpha-\beta)}=\frac{2\Lambda_{6}v^{2}}{m_{H}^{2}-m_{h}^{2}},\,\,\,\,\,\,\,\,c_{2(\alpha-\beta)}=\frac{(\Lambda_{1}-\Lambda_{5})v^{2}-m_{A}^{2}}{m_{H}^{2}-m_{h}^{2}}\,.$
(A.27)

A useful quantity (appearing into Higgs couplings) is

$c^{2}_{\alpha-\beta}=\frac{\Lambda_{1}v^{2}-m_{h}^{2}}{m_{H}^{2}-m_{h}^{2}}=\frac{\Lambda^{2}_{6}v^{4}}{(m_{H}^{2}-m_{h}^{2})(m_{H}^{2}-\Lambda_{1}v^{2})}\,,$
(A.28)

while some of the previous equations can be rearranged as

$\displaystyle\Lambda_{1}v^{2}$
$\displaystyle=m_{h}^{2}s_{\alpha-\beta}^{2}+m_{H}^{2}c_{\alpha-\beta}^{2}\,,$
(A.29) $\displaystyle\Lambda_{6}v^{2}$
$\displaystyle=(m_{h}^{2}-m_{H}^{2})\,s_{\alpha-\beta}c_{\alpha-\beta}\,,$
(A.30) $\displaystyle m_{A}^{2}+\Lambda_{5}v^{2}$
$\displaystyle=m_{H}^{2}s_{\alpha-\beta}^{2}+m_{h}^{2}c_{\alpha-\beta}^{2}\,.$
(A.31)

Finally, Eq. (2.11) provides the embedding of the physical mass eigenstates in
the original basis of the PQ-Higgs doublets.

Note that $s_{\alpha-\beta}=0$ corresponds to $\Lambda_{6}=0$ . In this
scenario the field $H$ is aligned to the electroweak VEV, making $H$ to be SM-
like with $m^{2}_{H}=\Lambda_{1}v^{2}$ [49]. In the opposite case where
$c_{\alpha-\beta}=0$ we have instead $m_{h}^{2}=\Lambda_{1}v^{2}$,
corresponding to a SM-like $h$. This can be obtained by taking $m_{H}^{2}\gg
v^{2}$, which from Eq. (A.28) that corresponds to $c_{\alpha-\beta}\sim
1/m_{H}^{2}$. This scenario is called alignment with decoupling [50]. From Eq.
(A.31) one also obtains

$m_{A}^{2}+\Lambda_{5}v^{2}\simeq m_{H}^{2}\,.$ (A.32)

Hence, in the decoupling scenario the scalar masses become degenerate,
$m_{H}^{2}\simeq m_{A}^{2}\simeq m_{H^{\pm}}^{2}$. This is instead not the
case for the alignment without decoupling, where the scalar masses remain at
the electroweak scale and, in general, are non-degenerate. Such a case
requires $\Lambda_{6}\to 0$ and it can lead to a non-SM-like Higgs either
above or below $125\,$GeV, see [49, 50] respectively.

##  Appendix B 4-fermion operators

We list here the 4-fermion operators which arise in Eq. (2.18), upon
integrating out at tree level the radial modes of the PQ-2HDM. In the quark
sector we obtain

$\displaystyle\mathcal{L}^{\rm 4\text{-}quark}_{\rm EFT}$
$\displaystyle=(C_{ijkl}^{\overline{u}_{L}u_{R}\overline{u}_{L}u_{R}}\overline{u}_{Li}u_{Rj}\overline{u}_{Lk}u_{Rl}+\text{h.c.})+(C_{ijkl}^{\overline{d}_{L}d_{R}\overline{d}_{L}d_{R}}\overline{d}_{Li}d_{Rj}\overline{d}_{Lk}d_{Rl}+\text{h.c.})$
$\displaystyle+(C_{ijkl}^{\overline{u}_{L}u_{R}\overline{d}_{L}d_{R}}\overline{u}_{Li}u_{Rj}\overline{d}_{Lk}d_{Rl}+\text{h.c.})+(C_{ijkl}^{\overline{u}_{L}u_{R}\overline{d}_{R}d_{L}}\overline{u}_{Li}u_{Rj}\overline{d}_{Rk}d_{Ll}+\text{h.c.})$
$\displaystyle+C_{ijkl}^{\overline{u}_{L}u_{R}\overline{u}_{R}u_{L}}\overline{u}_{Li}u_{Rj}\overline{u}_{Rk}u_{Ll}+C_{ijkl}^{\overline{d}_{L}d_{R}\overline{d}_{R}d_{L}}\overline{d}_{Li}d_{Rj}\overline{d}_{Rk}d_{Ll}$
$\displaystyle+(C_{ijkl}^{\overline{u}_{L}d_{R}\overline{d}_{L}u_{R}}\overline{u}_{Li}d_{Rj}\overline{d}_{Lk}u_{Rl}+\text{h.c.})+C_{ijkl}^{\overline{u}_{L}d_{R}\overline{d}_{R}u_{L}}\overline{u}_{Li}d_{Rj}\overline{d}_{Rk}u_{Ll}$
$\displaystyle+C_{ijkl}^{\overline{d}_{L}u_{R}\overline{u}_{R}d_{L}}\overline{d}_{Li}u_{Rj}\overline{u}_{Rk}d_{Ll}\,,$
(B.1)

where the 4-quark Wilson coefficients are found to be

$\displaystyle C^{\overline{u}_{L}u_{R}\overline{u}_{L}u_{R}}_{ijkl}$
$\displaystyle=\frac{1}{2m_{H}^{2}}C^{H_{u}}_{ij}C^{H_{u}}_{kl}+\frac{1}{2m_{A}^{2}}C^{A_{u}}_{ij}C^{A_{u}}_{kl}+\frac{1}{2m_{h}^{2}}C^{h_{u}}_{ij}C^{h_{u}}_{kl}\,,$
(B.2) $\displaystyle C^{\overline{u}_{L}u_{R}\overline{d}_{L}d_{R}}_{ijkl}$
$\displaystyle=\frac{1}{2m_{H}^{2}}C^{H_{u}}_{ij}C^{H_{d}}_{kl}+\frac{1}{2m_{A}^{2}}C^{A_{u}}_{ij}C^{A_{d}}_{kl}+\frac{1}{2m_{h}^{2}}C^{h_{u}}_{ij}C^{h_{d}}_{kl}\,,$
(B.3) $\displaystyle C^{\overline{d}_{L}d_{R}\overline{d}_{L}d_{R}}_{ijkl}$
$\displaystyle=\frac{1}{2m_{H}^{2}}C^{H_{d}}_{ij}C^{H_{d}}_{kl}+\frac{1}{2m_{A}^{2}}C^{A_{d}}_{ij}C^{A_{d}}_{kl}+\frac{1}{2m_{h}^{2}}C^{h_{d}}_{ij}C^{h_{d}}_{kl}\,,$
(B.4) $\displaystyle C^{\overline{u}_{L}u_{R}\overline{d}_{R}d_{L}}_{ijkl}$
$\displaystyle=\frac{1}{2m_{H}^{2}}C^{H_{u}}_{ij}(C^{H_{d}}_{lk})^{*}+\frac{1}{2m_{A}^{2}}C^{A_{u}}_{ij}(C^{A_{d}}_{lk})^{*}+\frac{1}{2m_{h}^{2}}C^{h_{u}}_{ij}(C^{h_{d}}_{lk})^{*}\,,$
(B.5) $\displaystyle C^{\overline{u}_{L}u_{R}\overline{u}_{R}u_{L}}_{ijkl}$
$\displaystyle=\frac{1}{m_{H}^{2}}C^{H_{u}}_{ij}(C^{H_{u}}_{lk})^{*}+\frac{1}{m_{A}^{2}}C^{A_{u}}_{ij}(C^{A_{u}}_{lk})^{*}+\frac{1}{m_{h}^{2}}C^{h_{u}}_{ij}(C^{h_{u}}_{lk})^{*}\,,$
(B.6) $\displaystyle C^{\overline{d}_{L}d_{R}\overline{d}_{R}d_{L}}_{ijkl}$
$\displaystyle=\frac{1}{m_{H}^{2}}C^{H_{d}}_{ij}(C^{H_{d}}_{lk})^{*}+\frac{1}{m_{A}^{2}}C^{A_{d}}_{ij}(C^{A_{d}}_{lk})^{*}+\frac{1}{m_{h}^{2}}C^{h_{d}}_{ij}(C^{h_{d}}_{lk})^{*}\,,$
(B.7) $\displaystyle C^{\overline{u}_{L}d_{R}\overline{d}_{L}u_{R}}_{ijkl}$
$\displaystyle=\frac{2}{m_{H^{\pm}}^{2}}V_{nj}V^{*}_{mk}C^{H^{+}_{d}}_{in}C^{H^{-}_{u}}_{ml}\,,$
(B.8) $\displaystyle C^{\overline{u}_{L}d_{R}\overline{d}_{R}u_{L}}_{ijkl}$
$\displaystyle=\frac{2}{m_{H^{\pm}}^{2}}V_{nj}V_{mk}^{*}C^{H^{+}_{d}}_{in}(C^{H^{+}_{d}}_{lm})^{*}\,,$
(B.9) $\displaystyle C^{\overline{d}_{L}u_{R}\overline{u}_{R}d_{L}}_{ijkl}$
$\displaystyle=\frac{2}{m_{H^{\pm}}^{2}}V^{*}_{ni}V_{ml}C^{H^{-}_{u}}_{nj}(C^{H^{-}_{u}}_{mk})^{*}\,.$
(B.10)

For the semi-leptonic operators we find

$\displaystyle\mathcal{L}^{\rm semi\text{-}lept}_{\rm EFT}$
$\displaystyle=(C_{ijkl}^{\overline{q}_{L}q_{R}\overline{e}_{L}e_{R}}\overline{q}_{Li}q_{Rj}\overline{e}_{Lk}e_{Rl}+\text{h.c.})+(C_{ijkl}^{\overline{q}_{L}q_{R}\overline{e}_{R}e_{L}}\overline{q}_{Li}q_{Rj}\overline{e}_{Rk}e_{Ll}+\text{h.c.})$
$\displaystyle+(C_{ijkl}^{\overline{\nu}_{L}e_{R}\overline{d}_{R}u_{L}}\overline{\nu}_{Li}e_{Rj}\overline{d}_{Rk}u_{Ll}+\text{h.c.})+(C_{ijkl}^{\overline{\nu}_{L}e_{R}\overline{d}_{R}u_{L}}\overline{\nu}_{Li}e_{Rj}\overline{d}_{Rk}u_{Ll}+\text{h.c.})\,,$
(B.11)

with Wilson coefficients

$\displaystyle C^{\overline{u}_{L}u_{R}\overline{e}_{L}e_{R}}_{ijkl}$
$\displaystyle=\frac{1}{m_{H}^{2}}C^{H_{u}}_{ij}C^{H_{e}}_{kl}+\frac{1}{m_{A}^{2}}C^{A_{u}}_{ij}C^{A_{e}}_{kl}+\frac{1}{m_{h}^{2}}C^{h_{u}}_{ij}C^{h_{e}}_{kl}\,,$
(B.12) $\displaystyle C^{\overline{d}_{L}d_{R}\overline{e}_{L}e_{R}}_{ijkl}$
$\displaystyle=\frac{1}{m_{H}^{2}}C^{H_{d}}_{ij}C^{H_{e}}_{kl}+\frac{1}{m_{A}^{2}}C^{A_{d}}_{ij}C^{A_{e}}_{kl}+\frac{1}{m_{h}^{2}}C^{h_{d}}_{ij}C^{h_{e}}_{kl}\,,$
(B.13) $\displaystyle C^{\overline{u}_{L}u_{R}\overline{e}_{R}e_{L}}_{ijkl}$
$\displaystyle=\frac{1}{m_{H}^{2}}C^{H_{u}}_{ij}(C^{H_{e}}_{lk})^{*}+\frac{1}{m_{A}^{2}}C^{A_{u}}_{ij}(C^{A_{e}}_{lk})^{*}+\frac{1}{m_{h}^{2}}C^{h_{u}}_{ij}(C^{h_{e}}_{lk})^{*}\,,$
(B.14) $\displaystyle C^{\overline{d}_{L}d_{R}\overline{e}_{R}e_{L}}_{ijkl}$
$\displaystyle=\frac{1}{m_{H}^{2}}C^{H_{d}}_{ij}(C^{H_{e}}_{lk})^{*}+\frac{1}{m_{A}^{2}}C^{A_{d}}_{ij}(C^{A_{e}}_{lk})^{*}+\frac{1}{m_{h}^{2}}C^{h_{d}}_{ij}(C^{h_{e}}_{lk})^{*}\,,$
(B.15) $\displaystyle C^{\overline{\nu}_{L}e_{R}\overline{d}_{L}u_{R}}_{ijkl}$
$\displaystyle=\frac{2}{m_{H^{\pm}}^{2}}U^{*}_{ni}V_{mk}^{*}C^{H^{+}_{e}}_{nj}C^{H^{-}_{u}}_{ml},$
(B.16) $\displaystyle C^{\overline{\nu}_{L}e_{R}\overline{d}_{R}u_{L}}_{ijkl}$
$\displaystyle=\frac{2}{m_{H^{\pm}}^{2}}U^{*}_{ni}V_{mk}^{*}C^{H^{-}_{e}}_{nj}(C^{H^{+}_{d}}_{lm})^{*}\,.$
(B.17)

Finally, the purely leptonic effective operators are

$\displaystyle\mathcal{L}^{\rm 4\text{-}lept}_{\rm EFT}$
$\displaystyle=(C_{ijkl}^{\overline{e}_{L}e_{R}\overline{e}_{L}e_{R}}\overline{e}_{Li}e_{Rj}\overline{e}_{Lk}e_{Rl}+\text{h.c.})+C_{ijkl}^{\overline{e}_{L}e_{R}\overline{e}_{R}e_{L}}\overline{e}_{Li}e_{Rj}\overline{e}_{Rk}e_{Ll}$
$\displaystyle+C_{ijkl}^{\overline{\nu}_{L}e_{R}\overline{e}_{R}\nu_{L}}\overline{\nu}_{Li}e_{Rj}\overline{e}_{Rk}\nu_{Ll}\,,$
(B.18)

with Wilson coefficients

$\displaystyle C^{\overline{e}_{L}e_{R}\overline{e}_{L}e_{R}}_{ijkl}$
$\displaystyle=\frac{1}{2m_{H}^{2}}C^{H_{e}}_{ij}C^{H_{e}}_{kl}+\frac{1}{2m_{A}^{2}}C^{A_{e}}_{ij}C^{A_{e}}_{kl}+\frac{1}{2m_{h}^{2}}C^{h_{e}}_{ij}C^{h_{e}}_{kl}\,,$
(B.19) $\displaystyle C^{\overline{e}_{L}e_{R}\overline{e}_{R}e_{L}}_{ijkl}$
$\displaystyle=\frac{1}{m_{H}^{2}}C^{H_{e}}_{ij}(C^{H_{e}}_{lk})^{*}+\frac{1}{m_{A}^{2}}C^{A_{e}}_{ij}(C^{A_{e}}_{lk})^{*}+\frac{1}{m_{h}^{2}}C^{h_{e}}_{ij}(C^{h_{e}}_{lk})^{*}\,,$
(B.20) $\displaystyle
C^{\overline{\nu}_{L}e_{R}\overline{e}_{R}\nu_{L}}_{ijkl}$
$\displaystyle=\frac{2}{m_{H^{\pm}}^{2}}U^{*}_{ni}U_{ml}C^{H^{+}_{e}}_{nj}(C^{H^{+}_{e}}_{mk})^{*}\,.$
(B.21)

##  Appendix C Meson mixing bounds

The neutral current 4-quark operators in Eq. (B.1) contribute to $\Delta F=2$
transitions, which are constrained by meson mixing observables. The effective
$\Delta F=2$ Hamiltonian is given by

$\mathcal{H}^{\Delta
F=2}_{\text{eff}}=\sum_{i}^{5}C_{i}Q_{i}+\sum_{i}^{5}\tilde{C}_{i}\tilde{Q}_{i}\,,$
(C.1)

with

$\displaystyle Q_{1}$
$\displaystyle=\overline{q}^{\alpha}_{L}\gamma^{\mu}q^{\prime\alpha}_{L}\overline{q}^{\beta}_{L}\gamma_{\mu}q^{\prime\beta}_{L}\,,$
(C.2) $\displaystyle Q_{2}$
$\displaystyle=\overline{q}^{\alpha}_{R}q^{\prime\alpha}_{L}\overline{q}^{\beta}_{R}q^{\prime\beta}_{L}\,,$
(C.3) $\displaystyle Q_{3}$
$\displaystyle=\overline{q}^{\alpha}_{R}q^{\prime\beta}_{L}\overline{q}^{\beta}_{R}q^{\prime\alpha}_{L}\,,$
(C.4) $\displaystyle Q_{4}$
$\displaystyle=\overline{q}^{\alpha}_{R}q^{\prime\alpha}_{L}\overline{q}^{\beta}_{L}q^{\prime\beta}_{R}\,,$
(C.5) $\displaystyle Q_{5}$
$\displaystyle=\overline{q}^{\alpha}_{R}q^{\prime\beta}_{L}\overline{q}^{\beta}_{L}q^{\prime\alpha}_{R}\,.$
(C.6)

Here, $\alpha$ and $\beta$ are color indices, while $q$ and $q^{\prime}$ refer
to two different quark flavours in the mesonic system. In the 2HDM case, that
is relevant for our analysis, only the operators $Q_{2}$ and $Q_{4}$ are
generated at tree level. However, it is well-known that QCD running effects
induce other operators, and give a sizable contribution to those already
present. Running effects will be taken into account by using the analytic
formulas provided in Refs. [79, 80, 81], respectively for $K$, $B$ and $D$
meson mixing. The master formula encoding the running of the Wilson
coefficients is given by

$C_{r}(\mu)=\sum_{i}\sum_{s}\left(b_{i}^{(r,s)}+\eta\,c^{(r,s)}_{i}\right)\eta^{a_{i}}\,C_{i}(M_{H})\,,$
(C.7)

where, $a,\,b,\,c$ are so-called “magic numbers”, which depend on the process
at hand and are taken from the above-mentioned references. The running of the
strong coupling is taken care by the factor
$\eta=\alpha_{s}(m_{H})/\alpha_{s}(M_{Z})$, where we take as an imput
$\alpha_{s}(M_{Z})$ and compute the QCD running at 2-loops up to the scale of
the heavy NP, $m_{H}$, by means of the code RunDec [104]. The low-energy scale
$\mu$ is set to $2$, $4.6$ and $2.8$ GeV respectively for the $K$, $B$ and $D$
meson systems, in order to match with the scale of the corresponding hadronic
matrix elements. The numerical inputs for the meson mixing analysis are given
in Table 3, where we defined the QCD running coefficients
$\eta_{i}=C_{i}(\mu)/C_{i}(m_{H})$ and $\eta_{ij}=C_{i}(\mu)/C_{j}(m_{H})$,
and we have taken as a reference value $m_{H}=1$ TeV.

| $K^{0}-\overline{K}^{0}$ | $B^{0}-\overline{B}^{0}$ | $B^{0}_{s}-\overline{B}^{0}_{s}$ | $D-\overline{D}$
---|---|---|---|---
$\Delta m\,$[GeV]$\times 10^{-15}$ | $3.484\pm 0.006$ | $333.4\pm 1.3$ | $11693\pm 8$ | $6.7\pm 0.8$
$M_{M}\,[\text{GeV}]$ | 0.497 | 5.279 | 5.367 | 1.864
$f_{M}\,[\text{GeV}]$ | 0.156 | 0.192 | 0.228 | 0.209
$B_{2}$ | 0.50 | 0.77 | 0.82 | 0.66
$B_{4}$ | 0.92 | 1.08 | 1.03 | 0.84
$\eta_{2}$ | 2.50 | 1.99 | 1.99 | 2.30
$\eta_{32}$ | $-1.70\times 10^{-3}$ | $-0.024$ | $-0.024$ | $-0.013$
$\eta_{4}$ | 4.78 | 3.10 | 3.10 | 3.97
$\eta_{45}$ | 0.18 | 0.086 | 0.086 | 0.14
Table 3: Numerical inputs for the analysis of $\Delta F=2$ observables,
including experimental values for meson-antimeson mass differences and meson
masses [105], meson decay constants [83], bag parameters for
$K^{0}-\overline{K}^{0}$ [83], $B_{d}^{0}-\overline{B}_{d}^{0}$ and
$B_{s}^{0}-\overline{B}_{s}^{0}$ [106], $D^{0}-\overline{D}^{0}$ [107], as
well as QCD running coefficients (see text).

Taking into account the renormalization group evolution of the Wilson
coefficient for $\Delta F=2$ transition, the meson mixing amplitudes read [79]

$\displaystyle M_{ij}=\frac{\bra{X}\mathcal{H}^{\Delta
F=2}_{\text{eff}}\ket{\overline{X}}}{2M_{X}}$
$\displaystyle=M_{X}f_{X}^{2}\left[-\frac{5}{24}\frac{M_{X}^{2}}{\left(m^{q}_{i}+m^{q}_{j}\right)^{2}}\eta_{2}B_{2}(\mu)\left((C_{2})_{ij}+(C_{2}^{*})_{ji}\right)\right.$
$\displaystyle\left.+\left(\frac{1}{24}+\frac{1}{4}\frac{M_{X}^{2}}{\left(m^{q}_{i}+m^{q}_{j}\right)^{2}}\right)\eta_{4}B_{4}(\mu)(C_{4})_{ij}\right]\,,$
(C.8)

where $B_{i}$ are the so-called bag parameters which are calculated in lattice
QCD (see Table 3 for the employed numerical values). Including also the
contribution of NP, the meson-antimeson mass difference is given by
(neglecting CP-violating phases)

$\Delta m=2|M_{12}^{\text{SM}}+M_{12}^{\text{NP}}|\,.$ (C.9)

Taking into account that the new physics contribution can be either positive
or negative, and that the SM uncertainty is much larger than the experimental
one (which is henceforth neglected), we conservatively consider the $2\sigma$
bounds:

$\displaystyle\Delta m^{exp}=\Delta m^{\text{SM}}-2\sigma_{\Delta
m}^{\text{SM}}+2|M_{12}^{\text{NP}}|\,,\quad\text{for}\quad M_{12}>0\,,$
(C.10) $\displaystyle\Delta m^{exp}=\Delta m^{\text{SM}}+2\sigma_{\Delta
m}^{\text{SM}}-2|M_{12}^{\text{NP}}|\,,\quad\text{for}\quad M_{12}<0\,.$
(C.11)

Regarding theoretical predictions, for the $B_{d}^{0}-\overline{B}_{d}^{0}$
and $B_{s}^{0}-\overline{B}_{s}^{0}$ mixing we will employ the values $\Delta
m_{d}=3.51^{+0.15}_{-0.24}\times 10^{-13}\,$GeV and $\Delta
m_{s}=1.211^{+0.046}_{-0.079}\times 10^{-11}\,$GeV from Ref. [82], which are
obtained by averaging lattice QCD and sum rules results. For
$K^{0}-\overline{K}^{0}$ mixing we take $\Delta m_{K}=(3.3\pm 1.2)\times
10^{-15}\,\text{GeV}$, that is the FLAG21 value [83] (see also [84]). For
$D^{0}-\overline{D}^{0}$ mixing a reliable SM prediction is still lacking, and
hence we will just impose that new physics does not overshoot the experimental
value at $2\sigma$. Taking all of the above into account and applying Eqs.
(C.10)–(C.11) we obtain

$\displaystyle\left|-0.16\,\left((C^{d}_{2})_{12}+(C^{d}_{2})^{*}_{21}\right)+0.69\,(C^{d}_{4})_{12}\right|=\left[2.0,\,2.5\right]\cdot
10^{-15}\,\text{GeV}\,,$ (C.12)
$\displaystyle\left|-0.20\,\left((C^{d}_{2})_{13}+(C^{d}_{2})^{*}_{31}\right)+0.57\,(C^{d}_{4})_{13}\right|=\left[4.6,\,3.0\right]\cdot
10^{-14}\,\text{GeV}\,,$ (C.13)
$\displaystyle\left|-0.28\,\left((C^{d}_{2})_{23}+(C^{d}_{2})^{*}_{32}\right)+0.82\,(C^{d}_{4})_{23}\right|=\left[1.3,\,1.2\right]\cdot
10^{-12}\,\text{GeV}\,,$ (C.14)

where the number to the left is for the case in which the new physics
contribution is negative, while the number on the right is for the case in
which it is positive. For $D$ mixing instead the bound reads

$\left|-0.16\,\left((C^{u}_{2})_{12}+(C^{u}_{2})^{*}_{21}\right)+0.53\,(C^{u}_{4})_{12}\right|=8.3\cdot
10^{-15}\,\text{GeV}\,.$ (C.15)

In this Appendix we have fixed as a reference the running of the $d$($u$)-type
sector to $1(10)$ TeV, while in the numerical analysis in Section 3 we have
properly taken into account running effects as a function of $m_{H}$. Note
that employing model-independent constraints on individual Wilson coefficients
would lead to an incorrect bound for the 2HDM case, since a specific
combination of Wilson coefficients is actually constrained.

The operators associated to the Wilson coefficients $C_{2}^{f}$ and
$C_{4}^{f}$ are obtained from Eq. (B.1), with the conditions $i=k$, $j=l$ and
$i\neq j$. Using this information, one can rewrite the Wilson coefficients in
terms of $\epsilon_{ij}^{f}$ as follows

$\displaystyle C_{2}^{d}$
$\displaystyle=\frac{(\epsilon_{ij}^{d})^{2}}{4c_{\beta}^{2}}\left(\frac{m_{A}^{2}-m_{H}^{2}}{m_{H}^{2}m_{A}^{2}}+c_{\alpha-\beta}^{2}\frac{m_{H}^{2}-m_{h}^{2}}{m_{H}^{2}m_{h}^{2}}\right)\,,$
(C.16) $\displaystyle C_{4}^{d}$
$\displaystyle=\frac{\epsilon_{ij}^{d\,}\epsilon_{ji}^{f\,*}}{2c_{\beta}^{2}}\left(\frac{m_{A}^{2}+m_{H}^{2}}{m_{H}^{2}m_{A}^{2}}+c_{\alpha-\beta}^{2}\frac{m_{H}^{2}-m_{h}^{2}}{m_{H}^{2}m_{h}^{2}}\right)\,,$
(C.17)

and similarly for the up-quark sector, which is obtained by replacing
$c_{\beta}\to s_{\beta}$ in the above formulae. Since $c_{\alpha-\beta}$ is
bounded by Higgs observables, the second term in Eqs. (C.16)–(C.17) is
somewhat suppressed. In fact, this becomes obvious in the case of alignment,
$c_{\alpha-\beta}\simeq 0$, with and without decoupling. In the alignment
scenarios the Wilson coefficients can be further simplified as

$\displaystyle C_{2}^{f}=0\,$ $\displaystyle,\quad
C_{4}^{f}=\frac{\epsilon_{ij}^{f}\epsilon_{ji}^{f\,*}}{c_{\beta}^{2}\,m_{H}^{2}}\,,\quad\quad\quad\quad\quad\
\text{(with decoupling)}$ (C.18) $\displaystyle
C_{2}^{f}=\frac{\epsilon_{ij}^{f\,2}}{4c_{\beta}^{2}}\frac{m_{A}^{2}-m_{H}^{2}}{m_{H}^{2}m_{A}^{2}}\,$
$\displaystyle,\quad
C_{4}^{f}=\frac{\epsilon_{ij}^{f\,}\epsilon_{ji}^{f\,*}}{2c_{\beta}^{2}}\frac{m_{A}^{2}+m_{H}^{2}}{m_{H}^{2}m_{A}^{2}}\,,\quad\text{(without
decoupling).}$ (C.19)

Plugging in these expressions in Eqs. (C.12)–(C.15) and employing the
definitions in Eqs. (2.3)–(2.3) it is straightforward to obtain the meson
mixing bounds for the alignment scenarios, see Eqs. (3.1)–(3.4) .

##  Appendix D Model-dependent bounds on the PQ-2HDM

Figure 4: Left (right) panel plots show the meson mixing bounds for the M1
(M4) model, with the solid (dashed) line corresponding to a positive
(negative) new physics contribution. The top panels correspond to
$c_{\alpha-\beta}=0$ and the bottom ones to $c_{\alpha-\beta}=0.1$. The black
vertical line in the bottom panels indicates the value of $m_{H}$ for which
the condition $c_{\alpha-\beta}=0.1$ saturates the perturbativity of the
$\Lambda_{1}$ parameter, i.e. $\Lambda_{1}\simeq 4\pi$. LHC direct searches
limits are displayed as well, taking as a reference those for the case of the
type-II 2HDM (grey shaded region). Figure 5: Left (right) panel plots show the
meson mixing bounds due to the SM-like Higgs for the M1 (M4) model, with the
solid (dashed) line corresponding to a positive (negative) new physics
contribution. The top panels correspond to a decoupled $m_{H}$ and the bottom
ones to $m_{H}=1$ TeV.

In this Appendix we analyze other aspects of the PQ-2HDM which do not rely on
the IR/UV connection discussed in this paper, but we assume a certain
structure for the PQ charges and the mixing matrices. A possible ansatz is
given for instance by $V_{d_{L,R}}=V_{\text{CKM}}$ and
$V_{u_{L,R}}=\mathbf{1}$. Considering the M1 and M4 models introduced in
Section 2, we can see from Fig. 4 that the bounds become much weaker once we
introduce this CKM-like scenario. In the M1 model they are around $200\,$GeV.
In fact, in this case the strongest bound does not come from $K$-mixing, since
the PQ charge structure is such that an $s$-$d$ transition is obtained after
two CKM rotations of strength $\sim\lambda^{3}\times\lambda^{2}$ and hence it
is somewhat protected by the 2+1 structure of the PQ charges. In the M4 model
instead, due to the PQ charge structure, the $K$-mixing leads to the strongest
bound, being the $s$-$d$ transition of order $\lambda$.

Note that in the meson mixing bounds (cf. Eqs. (3.1)–(3.4)) there could be a
cancellation for $c_{\alpha-\beta}\neq 0$ (or $m_{H}\neq m_{A}$) which
correspond to the bottom plots in Fig. 4. Hence, meson mixing does not rule
out a relevant portion of the parameter space. For the M1 the model this is of
the order of $200\,$GeV, while M4 is still in the reach of future LHC
searches. Also, in the eventual discovery of a heavy Higgs, if we observe as
well flavour mixing, the observed channel could help in reconstructing the
flavour structure of the 2HDM.

For completeness, and for comparison with Ref. [37], we also show in Fig. 5
the meson mixing bounds in the $(c_{\alpha-\beta},t_{\beta})$ plane, for a
fixed value of $m_{H}$. In particular, the top panel plots correspond to a
mass parameter that is large enough to be considered in the decoupling regime,
while the bottom ones are for $m_{H}=1\,$TeV. Here, we display the bound
imposed by the SM-like Higgs component contributing to meson mixing. Note that
it is not possible to decouple the heavier Higgs without sending to zero the
mixing angle of the CP-even states. If we fix the mass of the second Higgs at
$1\,$TeV, as in the bottom row of Fig. 5, then a possible cancellation in the
meson mixing formula (cf. Eqs. (3.1)–(3.4)) is at play. This is an interesting
parameter space region, since it can be potentially probed by Higgs flavour
violating decays, such as $h\to bs$. As a final comment, we note that meson
mixing constraints in the $(c_{\alpha-\beta},t_{\beta})$ plane are comparable
to those obtained from lepton flavour violating observables (see e.g. [37]).

##  Appendix E Lepton flavour violating observables

In this Appendix we provide a compendium of formulae for LFV processes in the
parametrization that allows us to make a direct connection between the IR and
UV sources of flavour violation, including the leptonic decays $l_{i}\to
3l_{j}$, $l_{i}\to l_{j}l_{k}l_{k}$ and $l_{i}\to\gamma l_{j}$ as well as the
LFV Higgs decays $h\to l_{i}l_{j}$. Employing the IR/UV relation in Eq.
(2.24), the contribution of the heavy radial modes of the PQ-2HDM to the
leptonic decays, $l_{i}\to 3l_{j}$, can be written as [87]

$\displaystyle\text{BR}(l_{i}\to 3l_{j})$
$\displaystyle=\frac{m_{\tau}^{5}}{96\pi^{3}\Gamma_{\tau}\,m_{H}^{4}}\left(\frac{m_{j}}{s_{2\beta}v}\right)^{4}(C^{A}_{e})_{jj}^{2}$
$\displaystyle\times\left[\left|(C^{R}_{e})_{ji}+\frac{m_{\tau}}{m_{j}}(C^{L}_{e})_{ji}\right|^{2}+\left|(C^{L}_{e})_{ji}+\frac{m_{\tau}}{m_{e}}(C^{R}_{e})_{ji}\right|^{2}\right]\,,$
(E.1)

where we have neglected the mass of the final state lepton. This equation
contains a combination of diagonal and off-diagonal components, thus making
the IR/UV connection a bit more involved.666In this case one could fix the
diagonal axion coupling to electrons to the value suggested by the cooling
hints [45], and then impose the LFV bounds in analogy to what done in Section
3.3 for the quark sector. Again, here we consider $f_{a}\gg m_{H}$, so that
the direct axion contribution to this observable is negligible.

Within 2+1 PQ models (featuring axion flavour violation either in the LH or RH
sector) the expression above gets further simplified into

$\text{BR}(l_{i}\to
3l_{j})=\frac{m_{f_{i}}^{5}}{96\pi^{3}\Gamma_{i}\,m_{H}^{4}}\left(\frac{\sqrt{m_{f_{j}}m_{f_{i}}}}{s_{2\beta}v}\right)^{4}(C^{A}_{e})_{jj}^{2}\left|(C^{R,L}_{e})_{ji}\right|^{2}\,.$
(E.2)

For $l_{i}\to l_{j}\overline{l}_{k}l_{k}$ decays we obtain instead (always in
the case of 2+1 PQ models)

$\displaystyle\text{BR}(l_{i}\to l_{j}\overline{l}_{k}l_{k})$
$\displaystyle=\frac{m_{f_{i}}^{5}}{96\pi^{3}\Gamma_{i}\,m_{H}^{4}}\left(\frac{m_{f_{i}}}{s_{2\beta}v}\right)^{4}\left[(C^{A}_{e})_{kk}^{2}\left|(C^{R,L}_{e})_{ji}\right|^{2}\frac{m_{k}^{2}}{m_{i}^{2}}\left(1+\frac{m_{j}^{2}}{m_{i}^{2}}\right)\right.$
$\displaystyle+\left.\left|(C^{R,L}_{e})_{jk}\right|^{2}\left|(C^{R,L}_{e})_{ki}\right|^{2}\left(\frac{m_{k}^{4}}{m_{i}^{4}}+\frac{m_{j}^{2}}{m_{i}^{2}}\right)\right]\,.$
(E.3)

Following [87] we can also write the one-loop contribution of the PQ-2HDM to
the radiative LFV decays as

$\text{BR}(l_{i}\to
l_{j}\gamma)=\frac{m_{l_{i}}^{5}}{4\pi\Gamma_{l_{i}}}\left(|C_{R}^{l_{j},l_{i}}|^{2}+|C_{L}^{l_{j},l_{i}}|^{2}\right)\,,$
(E.4)

in terms of the coefficients $C_{L(R)}^{l_{j},l_{i}}=C_{L(R),+}+C_{L(R),0}$,
denoting the sum of the charged and neutral contributions. For simplicity we
will give them here for the case of 2+1 PQ models. The charged contribution
reads

$\displaystyle
C_{R,+}^{l_{j},l_{i}}=\frac{e}{384\pi^{2}m_{H^{\pm}}^{2}}\left|C^{H^{+}_{e}}\right|^{2}_{ji}\,,\qquad
C_{L,+}^{l_{j},l_{i}}=\frac{e}{384\pi^{2}m_{H^{\pm}}^{2}}\frac{m_{l_{j}}}{m_{l_{i}}}\left|C^{H^{+}_{e}}\right|^{2}_{ji}\,,$
(E.5)

which we did not expand in terms of the axion couplings, since it includes
both diagonal and off-diagonal couplings. The neutral contribution reads
instead (working in the decoupling limit, $m_{H}^{2}\simeq m_{A}^{2}$)

$\displaystyle C_{L(R),0}^{l_{j},l_{i}}$
$\displaystyle=\frac{-e}{96\pi^{2}\,m_{H}^{2}}\left(\frac{\sqrt{m_{e_{i}}m_{e_{j}}}}{vc_{\beta}s_{\beta}}\right)^{2}\Bigg{[}2(C^{A}_{e})_{ii}(C^{L(R)}_{e})_{ji}+\left(1+\frac{m_{e_{j}}^{2}}{m_{e_{i}}^{2}}\right)(C^{A}_{e})_{jj}(C^{L(R)}_{e})_{ji}$
$\displaystyle-\left(1+\frac{m_{e_{k}}^{2}}{m_{e_{i}}^{2}}\right)(C^{L(R)}_{e})_{jk}(C^{L\,(R)}_{e})^{*}_{ik}\Bigg{]}\,,$
(E.6) $\displaystyle C_{R(L),0}^{l_{j},l_{i}}$
$\displaystyle=\frac{e}{96\pi^{2}\,m_{H}^{2}}\left(\frac{\sqrt{m_{e_{i}}m_{e_{j}}}}{vc_{\beta}s_{\beta}}\right)^{2}\frac{m_{e_{j}}}{m_{e_{i}}}\Bigg{[}2(C^{A}_{e})_{jj}(C^{L(R)}_{e})_{ji}+\left(1+\frac{m_{e_{i}}^{2}}{m_{e_{j}}^{2}}\right)(C^{A}_{e})_{ii}(C^{L(R)}_{e})_{ji}$
$\displaystyle-\left(1+\frac{m_{e_{k}}^{2}}{m_{e_{j}}^{2}}\right)(C^{L(R)}_{e})_{jk}(C^{L(R)}_{e})^{*}_{ik}\Bigg{]}\,,$
(E.7)

where $k\neq i,\,j$ is not summed over. Note that due to the structure of the
loop contribution, it is not possible to obtain a straightforward connection
with LFV axion observables. However, the bounds from radiative LFV decays such
as $\mu\to e\gamma$ are very relevant for constraining the parameter space and
they should be eventually taken into acccount in a global analysis.

We finally consider LFV Higgs decays, whose connection with LFV axion
processes was discussed in Section 3.4. In this Appendix we provide full
formulae, including the effects of a modified Higgs width and show that also
flavour conserving Higgs decays can be related to axion physics. The complete
formula for the LFV Higgs decays can be written as (in terms of LFV axion
couplings, exploiting Eq. (2.24))

$\displaystyle\text{BR}(h\to l_{i}l_{j})$
$\displaystyle=\frac{m_{h}}{16\pi\Gamma_{h}}\left(\frac{c_{\alpha-\beta}}{vs_{\beta}c_{\beta}}\right)^{2}\Big{[}\left(m_{l_{i}}^{2}+m_{l_{j}}^{2}\right)\left(|(C^{L}_{e})_{ij}|^{2}+|(C^{R}_{e})_{ij}|^{2}\right)$
$\displaystyle-4m_{l_{i}}m_{l_{j}}\mbox{Re}\,{\left((C^{L}_{e})_{ij}(C^{R}_{e})^{*}_{ij}\right)}\Big{]}\,,$
(E.8)

where the Higgs width is modified as [37]

$\Gamma_{h}\simeq\frac{k_{h}^{2}}{1-B^{\text{BSM}}_{h}}\Gamma_{\text{SM}}\,,$
(E.9)

with $B^{\text{BSM}}_{h}$ denoting the branching ratio of the beyond-the-SM
contributions to the Higgs decay width, $\Gamma^{\text{SM}}_{h}\simeq
4.1\,$MeV and $k_{h}^{2}$ parametrizing the effect due to Higgs couplings
modifications. The latter is defined as [37]

$\displaystyle k_{h}^{2}$ $\displaystyle\simeq
0.58k_{b}^{2}+0.22k_{W}^{2}+0.08k_{g}^{2}+0.06k_{\tau}^{2}+0.03k_{Z}^{2}+0.03k_{c}^{2}+2.3\times
10^{-3}k_{\gamma}^{2}$ $\displaystyle+1.5\times 10^{-3}k_{Z\gamma}^{2}+4\times
10^{-4}k_{s}^{2}+2.2\times 10^{-4}k_{\mu}^{2}\,,$ (E.10)

in terms of the $k$-parameters

$k_{f_{i}}\equiv\frac{\sqrt{2}m_{f}}{v}C^{h_{f}}_{ii}=-s_{\alpha-\beta}+\frac{c_{\alpha-\beta}}{c_{\beta}s_{\beta}}(C^{A}_{f})_{ii}\,,$
(E.11)

where in the last step we have exploited again Eq. (2.24). In particular, the
modified Higgs interactions with gauge couplings read

$\displaystyle k_{W}=k_{Z}=s_{\alpha-\beta}\,,$ (E.12) $\displaystyle
k^{2}_{Z\gamma}\simeq 0.00348k_{t}^{2}+1.121k_{W}^{2}-0.1249k_{t}k_{W}\,,$
(E.13) $\displaystyle k_{g}^{2}\simeq
1.04k_{t}^{2}+0.002k_{b}^{2}-0.04k_{b}k_{t}\,,$ (E.14) $\displaystyle
k_{\gamma}^{2}\simeq 1.59k_{W}^{2}+0.07k_{t}^{2}-0.67k_{W}k_{t}\,.$ (E.15)

Hence, the above expressions show that deviations in Higgs couplings can also
be related to diagonal axion couplings.

## References

  * [1] R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440–1443.
  * [2] R. D. Peccei and H. R. Quinn, “Constraints Imposed by CP Conservation in the Presence of Instantons,” Phys. Rev. D16 (1977) 1791–1797.
  * [3] S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40 (1978) 223–226.
  * [4] F. Wilczek, “Problem of Strong p and t Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * [5] M. Dine and W. Fischler, “The Not So Harmless Axion,” Phys. Lett. B 120 (1983) 137–141.
  * [6] L. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. B 120 (1983) 133–136.
  * [7] J. Preskill, M. B. Wise, and F. Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. B 120 (1983) 127–132.
  * [8] A. R. Zhitnitsky, “On Possible Suppression of the Axion Hadron Interactions. (In Russian),” Sov. J. Nucl. Phys. 31 (1980) 260. [Yad. Fiz.31,497(1980)].
  * [9] M. Dine, W. Fischler, and M. Srednicki, “A Simple Solution to the Strong CP Problem with a Harmless Axion,” Phys. Lett. B104 (1981) 199–202.
  * [10] J. E. Kim, “Weak Interaction Singlet and Strong CP Invariance,” Phys. Rev. Lett. 43 (1979) 103.
  * [11] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, “Can Confinement Ensure Natural CP Invariance of Strong Interactions?,” Nucl. Phys. B166 (1980) 493.
  * [12] R. D. Peccei, T. T. Wu, and T. Yanagida, “A Viable Axion Model,” Phys. Lett. B 172 (1986) 435–440.
  * [13] L. M. Krauss and F. Wilczek, “A Shortlived Axion Variant,” Phys. Lett. B 173 (1986) 189–192.
  * [14] L. M. Krauss and D. J. Nash, “A Viable Weak Iinteraction Axion?,” Phys. Lett. B 202 (1988) 560–567.
  * [15] M. Clemente, E. Berdermann, P. Kienle, H. Tsertos, W. Wagner, C. Kozhuharov, F. Bosch, and W. Konig, “Narrow positron lines from U - U and U - Th collisions,” Phys. Lett. B 137 (1984) 41–46.
  * [16] J. Schweppe et al., “Observation of a peak structure in positron spectra from U + Cm collisions,” Phys. Rev. Lett. 51 (1983) 2261–2264.
  * [17] Y. Ema, K. Hamaguchi, T. Moroi, and K. Nakayama, “Flaxion: a minimal extension to solve puzzles in the standard model,” JHEP 01 (2017) 096, arXiv:1612.05492 [hep-ph].
  * [18] L. Calibbi, F. Goertz, D. Redigolo, R. Ziegler, and J. Zupan, “Minimal axion model from flavor,” Phys. Rev. D 95 no. 9, (2017) 095009, arXiv:1612.08040 [hep-ph].
  * [19] F. Arias-Aragon and L. Merlo, “The Minimal Flavour Violating Axion,” JHEP 10 (2017) 168, arXiv:1709.07039 [hep-ph]. [Erratum: JHEP 11, 152 (2019)].
  * [20] F. Björkeroth, L. Di Luzio, F. Mescia, and E. Nardi, “$U(1)$ flavour symmetries as Peccei-Quinn symmetries,” JHEP 02 (2019) 133, arXiv:1811.09637 [hep-ph].
  * [21] L. Di Luzio, F. Mescia, E. Nardi, P. Panci, and R. Ziegler, “Astrophobic Axions,” Phys. Rev. Lett. 120 no. 26, (2018) 261803, arXiv:1712.04940 [hep-ph].
  * [22] F. Bjorkeroth, L. Di Luzio, F. Mescia, E. Nardi, P. Panci, and R. Ziegler, “Axion-electron decoupling in nucleophobic axion models,” Phys. Rev. D 101 no. 3, (2020) 035027, arXiv:1907.06575 [hep-ph].
  * [23] L. Darmé, L. Di Luzio, M. Giannotti, and E. Nardi, “Selective enhancement of the QCD axion couplings,” Phys. Rev. D 103 no. 1, (2021) 015034, arXiv:2010.15846 [hep-ph].
  * [24] L. Di Luzio, F. Mescia, E. Nardi, and S. Okawa, “Renormalization group effects in astrophobic axion models,” Phys. Rev. D 106 no. 5, (2022) 055016, arXiv:2205.15326 [hep-ph].
  * [25] M. Badziak and K. Harigaya, “Naturally astrophobic QCD axion,” arXiv:2301.09647 [hep-ph].
  * [26] F. Takahashi and W. Yin, “Hadrophobic Axion from GUT,” arXiv:2301.10757 [hep-ph].
  * [27] C. Cornella, P. Paradisi, and O. Sumensari, “Hunting for ALPs with Lepton Flavor Violation,” JHEP 01 (2020) 158, arXiv:1911.06279 [hep-ph].
  * [28] 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 no. 1, (2020) 015023, arXiv:2002.04623 [hep-ph].
  * [29] L. Calibbi, D. Redigolo, R. Ziegler, and J. Zupan, “Looking forward to lepton-flavor-violating ALPs,” JHEP 09 (2021) 173, arXiv:2006.04795 [hep-ph].
  * [30] M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, “Flavor probes of axion-like particles,” JHEP 09 (2022) 056, arXiv:2110.10698 [hep-ph].
  * [31] A. W. M. Guerrera and S. Rigolin, “Revisiting $K\rightarrow\pi a$ decays,” Eur. Phys. J. C 82 no. 3, (2022) 192, arXiv:2106.05910 [hep-ph].
  * [32] J. A. Gallo, A. W. M. Guerrera, S. Peñaranda, and S. Rigolin, “Leptonic meson decays into invisible ALP,” Nucl. Phys. B 979 (2022) 115791, arXiv:2111.02536 [hep-ph].
  * [33] Y. Jho, S. Knapen, and D. Redigolo, “Lepton-flavor violating axions at MEG II,” JHEP 10 (2022) 029, arXiv:2203.11222 [hep-ph].
  * [34] A. W. M. Guerrera and S. Rigolin, “ALP production in weak mesonic decays,” Fortsch. Phys. 2023 (11, 2022) 2200192, arXiv:2211.08343 [hep-ph].
  * [35] I. G. Irastorza and J. Redondo, “New experimental approaches in the search for axion-like particles,” Prog. Part. Nucl. Phys. 102 (2018) 89–159, arXiv:1801.08127 [hep-ph].
  * [36] P. Sikivie, “Invisible Axion Search Methods,” Rev. Mod. Phys. 93 no. 1, (2021) 015004, arXiv:2003.02206 [hep-ph].
  * [37] M. Badziak, G. Grilli di Cortona, M. Tabet, and R. Ziegler, “Flavor-violating Higgs decays and stellar cooling anomalies in axion models,” JHEP 10 (2021) 181, arXiv:2107.09708 [hep-ph].
  * [38] L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, “The landscape of QCD axion models,” Phys. Rept. 870 (2020) 1–117, arXiv:2003.01100 [hep-ph].
  * [39] L. Di Luzio, F. Mescia, and E. Nardi, “Redefining the Axion Window,” Phys. Rev. Lett. 118 no. 3, (2017) 031801, arXiv:1610.07593 [hep-ph].
  * [40] L. Di Luzio, F. Mescia, and E. Nardi, “Window for preferred axion models,” Phys. Rev. D 96 no. 7, (2017) 075003, arXiv:1705.05370 [hep-ph].
  * [41] K. Choi, S. H. Im, C. B. Park, and S. Yun, “Minimal Flavor Violation with Axion-like Particles,” JHEP 11 (2017) 070, arXiv:1708.00021 [hep-ph].
  * [42] M. Chala, G. Guedes, M. Ramos, and J. Santiago, “Running in the ALPs,” Eur. Phys. J. C 81 no. 2, (2021) 181, arXiv:2012.09017 [hep-ph].
  * [43] M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, “The Low-Energy Effective Theory of Axions and ALPs,” JHEP 04 (2021) 063, arXiv:2012.12272 [hep-ph].
  * [44] K. Choi, S. H. Im, H. J. Kim, and H. Seong, “Precision axion physics with running axion couplings,” JHEP 08 (2021) 058, arXiv:2106.05816 [hep-ph].
  * [45] L. Di Luzio, M. Fedele, M. Giannotti, F. Mescia, and E. Nardi, “Stellar evolution confronts axion models,” JCAP 02 (2022) 035, arXiv:2109.10368 [hep-ph].
  * [46] ATLAS Collaboration, “Combined measurements of Higgs boson production and decay using up to $139$ fb-1 of proton-proton collision data at $\sqrt{s}=13$ TeV collected with the ATLAS experiment,”.
  * [47] L. J. Hall and M. B. Wise, “Flavor Changing Higgs Boson Couplings,” Nucl. Phys. B 187 (1981) 397–408.
  * [48] J. F. Gunion and H. E. Haber, “The CP conserving two Higgs doublet model: The Approach to the decoupling limit,” Phys. Rev. D 67 (2003) 075019, arXiv:hep-ph/0207010.
  * [49] J. Bernon, J. F. Gunion, H. E. Haber, Y. Jiang, and S. Kraml, “Scrutinizing the alignment limit in two-Higgs-doublet models: mh=125 GeV,” Phys. Rev. D 92 no. 7, (2015) 075004, arXiv:1507.00933 [hep-ph].
  * [50] J. Bernon, J. F. Gunion, H. E. Haber, Y. Jiang, and S. Kraml, “Scrutinizing the alignment limit in two-Higgs-doublet models. II. mH=125 GeV,” Phys. Rev. D 93 no. 3, (2016) 035027, arXiv:1511.03682 [hep-ph].
  * [51] ATLAS Collaboration, M. Aaboud et al., “Search for additional heavy neutral Higgs and gauge bosons in the ditau final state produced in 36 fb-1 of pp collisions at $\sqrt{s}=13$ TeV with the ATLAS detector,” JHEP 01 (2018) 055, arXiv:1709.07242 [hep-ex].
  * [52] ATLAS Collaboration, M. Aaboud et al., “Search for $WW/WZ$ resonance production in $\ell\nu qq$ final states in $pp$ collisions at $\sqrt{s}=$ 13 TeV with the ATLAS detector,” JHEP 03 (2018) 042, arXiv:1710.07235 [hep-ex].
  * [53] ATLAS Collaboration, M. Aaboud et al., “Searches for heavy $ZZ$ and $ZW$ resonances in the $\ell\ell qq$ and $\nu\nu qq$ final states in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector,” JHEP 03 (2018) 009, arXiv:1708.09638 [hep-ex].
  * [54] ATLAS Collaboration, M. Aaboud et al., “Search for heavy ZZ resonances in the $\ell^{+}\ell^{-}\ell^{+}\ell^{-}$ and $\ell^{+}\ell^{-}\nu\overline{\nu}$ final states using proton–proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector,” Eur. Phys. J. C 78 no. 4, (2018) 293, arXiv:1712.06386 [hep-ex].
  * [55] ATLAS Collaboration, M. Aaboud et al., “Search for heavy resonances decaying into $WW$ in the $e\nu\mu\nu$ final state in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector,” Eur. Phys. J. C 78 no. 1, (2018) 24, arXiv:1710.01123 [hep-ex].
  * [56] ATLAS Collaboration, M. Aaboud et al., “Search for heavy resonances decaying into a $W$ or $Z$ boson and a Higgs boson in final states with leptons and $b$-jets in 36 fb-1 of $\sqrt{s}=13$ TeV $pp$ collisions with the ATLAS detector,” JHEP 03 (2018) 174, arXiv:1712.06518 [hep-ex]. [Erratum: JHEP 11, 051 (2018)].
  * [57] ATLAS Collaboration, M. Aaboud et al., “Search for a heavy Higgs boson decaying into a $Z$ boson and another heavy Higgs boson in the $\ell\ell bb$ final state in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector,” Phys. Lett. B 783 (2018) 392–414, arXiv:1804.01126 [hep-ex].
  * [58] ATLAS Collaboration, G. Aad et al., “Reconstruction and identification of boosted di-$\tau$ systems in a search for Higgs boson pairs using 13 TeV proton-proton collision data in ATLAS,” JHEP 11 (2020) 163, arXiv:2007.14811 [hep-ex].
  * [59] ATLAS Collaboration, G. Aad et al., “Search for heavy resonances decaying into a pair of Z bosons in the $\ell^{+}\ell^{-}\ell^{\prime+}\ell^{\prime-}$ and $\ell^{+}\ell^{-}\nu{{\overline{\nu}}}$ final states using 139 $\mathrm{fb}^{-1}$ of proton–proton collisions at $\sqrt{s}=13\,$TeV with the ATLAS detector,” Eur. Phys. J. C 81 no. 4, (2021) 332, arXiv:2009.14791 [hep-ex].
  * [60] ATLAS Collaboration, G. Aad et al., “Search for heavy Higgs bosons decaying into two tau leptons with the ATLAS detector using $pp$ collisions at $\sqrt{s}=13$ TeV,” Phys. Rev. Lett. 125 no. 5, (2020) 051801, arXiv:2002.12223 [hep-ex].
  * [61] CMS Collaboration, A. M. Sirunyan et al., “Search for a massive resonance decaying to a pair of Higgs bosons in the four b quark final state in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Phys. Lett. B 781 (2018) 244–269, arXiv:1710.04960 [hep-ex].
  * [62] CMS Collaboration, A. M. Sirunyan et al., “Search for Higgs boson pair production in events with two bottom quarks and two tau leptons in proton–proton collisions at $\sqrt{s}$ =13TeV,” Phys. Lett. B 778 (2018) 101–127, arXiv:1707.02909 [hep-ex].
  * [63] CMS Collaboration, A. M. Sirunyan et al., “Combination of searches for Higgs boson pair production in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Phys. Rev. Lett. 122 no. 12, (2019) 121803, arXiv:1811.09689 [hep-ex].
  * [64] CMS Collaboration, A. M. Sirunyan et al., “Search for a heavy Higgs boson decaying to a pair of W bosons in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” JHEP 03 (2020) 034, arXiv:1912.01594 [hep-ex].
  * [65] CMS Collaboration, A. M. Sirunyan et al., “Search for a heavy pseudoscalar Higgs boson decaying into a 125 GeV Higgs boson and a Z boson in final states with two tau and two light leptons at $\sqrt{s}=$ 13 TeV,” JHEP 03 (2020) 065, arXiv:1910.11634 [hep-ex].
  * [66] CMS Collaboration, A. M. Sirunyan et al., “Search for new neutral Higgs bosons through the H$\to$ ZA $\to\ell^{+}\ell^{-}\mathrm{b\overline{b}}$ process in pp collisions at $\sqrt{s}=$ 13 TeV,” JHEP 03 (2020) 055, arXiv:1911.03781 [hep-ex].
  * [67] CMS Collaboration, A. M. Sirunyan et al., “Search for heavy Higgs bosons decaying to a top quark pair in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” JHEP 04 (2020) 171, arXiv:1908.01115 [hep-ex]. [Erratum: JHEP 03, 187 (2022)].
  * [68] CMS Collaboration, A. M. Sirunyan et al., “Search for a heavy pseudoscalar boson decaying to a Z and a Higgs boson at $\sqrt{s}=$ 13 TeV,” Eur. Phys. J. C 79 no. 7, (2019) 564, arXiv:1903.00941 [hep-ex].
  * [69] CMS Collaboration, A. M. Sirunyan et al., “Search for resonant pair production of Higgs bosons in the $bbZZ$ channel in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Phys. Rev. D 102 no. 3, (2020) 032003, arXiv:2006.06391 [hep-ex].
  * [70] ATLAS Collaboration, G. Aad et al., “Search for dijet resonances in events with an isolated charged lepton using $\sqrt{s}=13$ TeV proton-proton collision data collected by the ATLAS detector,” JHEP 06 (2020) 151, arXiv:2002.11325 [hep-ex].
  * [71] ATLAS Collaboration, G. Aad et al., “Search for charged Higgs bosons decaying into a top quark and a bottom quark at $\sqrt{\mathrm{s}}$ = 13 TeV with the ATLAS detector,” JHEP 06 (2021) 145, arXiv:2102.10076 [hep-ex].
  * [72] CMS Collaboration, A. M. Sirunyan et al., “Search for charged Higgs bosons in the H± $\to$ $\tau^{\pm}\nu_{\tau}$ decay channel in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” JHEP 07 (2019) 142, arXiv:1903.04560 [hep-ex].
  * [73] CMS Collaboration, A. M. Sirunyan et al., “Search for physics beyond the standard model in multilepton final states in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” JHEP 03 (2020) 051, arXiv:1911.04968 [hep-ex].
  * [74] L. Wang, J. M. Yang, and Y. Zhang, “Two-Higgs-doublet models in light of current experiments: a brief review,” Commun. Theor. Phys. 74 no. 9, (2022) 097202, arXiv:2203.07244 [hep-ph].
  * [75] O. Atkinson, M. Black, A. Lenz, A. Rusov, and J. Wynne, “Cornering the Two Higgs Doublet Model Type II,” JHEP 04 (2022) 172, arXiv:2107.05650 [hep-ph].
  * [76] O. Atkinson, M. Black, C. Englert, A. Lenz, A. Rusov, and J. Wynne, “The flavourful present and future of 2HDMs at the collider energy frontier,” JHEP 11 (2022) 139, arXiv:2202.08807 [hep-ph].
  * [77] K. S. Babu and S. Jana, “Enhanced Di-Higgs Production in the Two Higgs Doublet Model,” JHEP 02 (2019) 193, arXiv:1812.11943 [hep-ph].
  * [78] T. P. Cheng and M. Sher, “Mass Matrix Ansatz and Flavor Nonconservation in Models with Multiple Higgs Doublets,” Phys. Rev. D 35 (1987) 3484.
  * [79] M. Ciuchini et al., “Delta M(K) and epsilon(K) in SUSY at the next-to-leading order,” JHEP 10 (1998) 008, arXiv:hep-ph/9808328.
  * [80] D. Becirevic, M. Ciuchini, E. Franco, V. Gimenez, G. Martinelli, A. Masiero, M. Papinutto, J. Reyes, and L. Silvestrini, “$B_{d}-\overline{B}_{d}$ mixing and the $B_{d}\to J/\psi K_{s}$ asymmetry in general SUSY models,” Nucl. Phys. B 634 (2002) 105–119, arXiv:hep-ph/0112303.
  * [81] UTfit Collaboration, M. Bona et al., “Model-independent constraints on $\Delta F=2$ operators and the scale of new physics,” JHEP 03 (2008) 049, arXiv:0707.0636 [hep-ph].
  * [82] L. Di Luzio, M. Kirk, A. Lenz, and T. Rauh, “$\Delta M_{s}$ theory precision confronts flavour anomalies,” JHEP 12 (2019) 009, arXiv:1909.11087 [hep-ph].
  * [83] Flavour Lattice Averaging Group (FLAG) Collaboration, Y. Aoki et al., “FLAG Review 2021,” Eur. Phys. J. C 82 no. 10, (2022) 869, arXiv:2111.09849 [hep-lat].
  * [84] J. Brod and M. Gorbahn, “Next-to-Next-to-Leading-Order Charm-Quark Contribution to the $CP$ Violation Parameter $\epsilon_{K}$ and $\Delta M_{K}$,” Phys. Rev. Lett. 108 (2012) 121801, arXiv:1108.2036 [hep-ph].
  * [85] G. C. Branco, W. Grimus, and L. Lavoura, “Relating the scalar flavor changing neutral couplings to the CKM matrix,” Phys. Lett. B 380 (1996) 119–126, arXiv:hep-ph/9601383.
  * [86] G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, “Theory and phenomenology of two-Higgs-doublet models,” Phys. Rept. 516 (2012) 1–102, arXiv:1106.0034 [hep-ph].
  * [87] A. Crivellin, A. Kokulu, and C. Greub, “Flavor-phenomenology of two-Higgs-doublet models with generic Yukawa structure,” Phys. Rev. D 87 no. 9, (2013) 094031, arXiv:1303.5877 [hep-ph].
  * [88] A. Crivellin, S. Najjari, and J. Rosiek, “Lepton Flavor Violation in the Standard Model with general Dimension-Six Operators,” JHEP 04 (2014) 167, arXiv:1312.0634 [hep-ph].
  * [89] BaBar Collaboration, B. Aubert et al., “Searches for Lepton Flavor Violation in the Decays tau+- —$>$ e+- gamma and tau+- —$>$ mu+- gamma,” Phys. Rev. Lett. 104 (2010) 021802, arXiv:0908.2381 [hep-ex].
  * [90] MEG Collaboration, J. Adam et al., “New constraint on the existence of the $\mu^{+}\to e^{+}\gamma$ decay,” Phys. Rev. Lett. 110 (2013) 201801, arXiv:1303.0754 [hep-ex].
  * [91] MEG Collaboration, A. M. Baldini et al., “Search for the lepton flavour violating decay $\mu^{+}\rightarrow\mathrm{e}^{+}\gamma$ with the full dataset of the MEG experiment,” Eur. Phys. J. C 76 no. 8, (2016) 434, arXiv:1605.05081 [hep-ex].
  * [92] Belle Collaboration, A. Abdesselam et al., “Search for lepton-flavor-violating tau-lepton decays to $\ell\gamma$ at Belle,” JHEP 10 (2021) 19, arXiv:2103.12994 [hep-ex].
  * [93] SINDRUM Collaboration, U. Bellgardt et al., “Search for the Decay mu+ —$>$ e+ e+ e-,” Nucl. Phys. B 299 (1988) 1–6.
  * [94] K. Hayasaka et al., “Search for Lepton Flavor Violating Tau Decays into Three Leptons with 719 Million Produced Tau+Tau- Pairs,” Phys. Lett. B 687 (2010) 139–143, arXiv:1001.3221 [hep-ex].
  * [95] CMS Collaboration, A. M. Sirunyan et al., “Search for lepton-flavor violating decays of the Higgs boson in the $\mu\tau$ and e$\tau$ final states in proton-proton collisions at $\sqrt{s}$ = 13 TeV,” Phys. Rev. D 104 no. 3, (2021) 032013, arXiv:2105.03007 [hep-ex].
  * [96] ATLAS Collaboration, “Searches for lepton-flavour-violating decays of the Higgs boson into $e\tau$ and $\mu\tau$ in $\sqrt{s}=13$ TeV $pp$ collisions with the ATLAS detector,” arXiv:2302.05225 [hep-ex].
  * [97] R. Harnik, J. Kopp, and J. Zupan, “Flavor Violating Higgs Decays,” JHEP 03 (2013) 026, arXiv:1209.1397 [hep-ph].
  * [98] Belle-II Collaboration, I. Adachi et al., “Search for lepton-flavor-violating $\tau$ decays to a lepton and an invisible boson at Belle II,” arXiv:2212.03634 [hep-ex].
  * [99] S. Bertolini, L. Di Luzio, H. Kolešová, and M. Malinský, “Massive neutrinos and invisible axion minimally connected,” Phys. Rev. D 91 no. 5, (2015) 055014, arXiv:1412.7105 [hep-ph].
  * [100] D. Espriu, F. Mescia, and A. Renau, “Axion-Higgs interplay in the two Higgs-doublet model,” Phys. Rev. D 92 no. 9, (2015) 095013, arXiv:1503.02953 [hep-ph].
  * [101] R. R. Volkas, A. J. Davies, and G. C. Joshi, “Naturalness of the Invisible Axion Model,” Phys. Lett. B 215 (1988) 133–138.
  * [102] R. Foot, A. Kobakhidze, K. L. McDonald, and R. R. Volkas, “Poincaré protection for a natural electroweak scale,” Phys. Rev. D 89 no. 11, (2014) 115018, arXiv:1310.0223 [hep-ph].
  * [103] S. Bertolini, L. Di Luzio, H. Kolešová, M. Malinský, and J. C. Vasquez, “Neutrino-axion-dilaton interconnection,” Phys. Rev. D 93 no. 1, (2016) 015009, arXiv:1510.03668 [hep-ph].
  * [104] F. Herren and M. Steinhauser, “Version 3 of RunDec and CRunDec,” Comput. Phys. Commun. 224 (2018) 333–345, arXiv:1703.03751 [hep-ph].
  * [105] Particle Data Group Collaboration, P. A. Zyla et al., “Review of Particle Physics,” PTEP 2020 no. 8, (2020) 083C01.
  * [106] R. J. Dowdall, C. T. H. Davies, R. R. Horgan, G. P. Lepage, C. J. Monahan, J. Shigemitsu, and M. Wingate, “Neutral B-meson mixing from full lattice QCD at the physical point,” Phys. Rev. D 100 no. 9, (2019) 094508, arXiv:1907.01025 [hep-lat].
  * [107] N. Carrasco et al., “$D^{0}-\overline{D}^{0}$ mixing in the standard model and beyond from $N_{f}$ =2 twisted mass QCD,” Phys. Rev. D 90 no. 1, (2014) 014502, arXiv:1403.7302 [hep-lat].