Leading Logs in QCD Axion Effective Field Theory

The axion is much lighter than all other degrees of freedom introduced by the Peccei-Quinn mechanism to solve the strong CP problem. It is therefore natural to use an effective field theory (EFT) to describe its interactions. Loop processes calculated in the EFT may, however, explicitly depend on the ultraviolet cutoff. In general the UV cutoff is not uniquely defined, but the dimensionful couplings suggest to identify it with the Peccei-Quinn symmetry-breaking scale. An example are $K \rightarrow \pi + a$ decays that will soon be tested to improved precision in NA62 and KOTO and whose amplitude is dominated by the term logarithmically dependent on the cutoff. In this paper, we critically examine the adequacy of using such a naive EFT approach to study loop processes by comparing EFT calculations with ones performed in complete QCD axion models. In DFSZ models, for example, the cutoff is found to be set by additional Higgs degrees of freedom and to therefore be much closer to the electroweak scale than to the Peccei-Quinn scale. In fact, there are non-trivial requirements on axion models where the cutoff scale of loop processes is close to the Peccei-Quinn scale, such that the naive EFT result is reproduced. This suggests that the existence of a suitable UV embedding may impose restrictions on axion EFTs. We provide an explicit construction of a model with suitable fermion couplings and find promising prospects for NA62 and IAXO.


Introduction
Effective Field Theories (EFTs) are extremely useful tools to test new physics scenarios in a model-independent way. Their main ingredients are an often unknown cutoff energy scale and a set of low-energy fields and symmetries from which the operator expansion is constructed. Some examples of EFTs include the Weak Effective Theory [1][2][3], the Standard Model Effective Field Theory [4,5], and the Higgs Effective Field Theory [6][7][8].
Nearly all phenomenological studies of axions and axion-like particles 1 are performed within an EFT where one considers the Standard Model (SM) extended by a pseudoscalar particle -the axion -that interacts with the SM via dimension-5 operators involving either two gauge bosons or two SM fermions.
The rationale behind the employed EFTs is that the (QCD) axion is a pseudo-Goldstone boson arising from the Peccei-Quinn (PQ) solution to the Strong CP Problem [90][91][92][93]. It is expected that the underlying U (1) PQ symmetry is broken spontaneously at a very high energy scale while the axion remains very light. The pseudo-Goldstone boson nature of the axion then sets the guiding principles for the construction of its EFT. Firstly, the shift symmetry a → a + const, resulting from the underlying U (1) PQ , restricts the interactions of the axion to be of purely derivative form. Secondly, the U (1) PQ symmetry must be anomalous under QCD in order for the strong CP problem to be solved. Non-perturbative effects then generate an axion potential and in particular a non-vanishing mass. Moreover, the QCD anomaly and similar anomalies under the other gauge groups result in characteristic couplings of the axion to two gauge bosons. Finally, it seems natural to identify the cutoff scale of the EFT with the PQ symmetry-breaking scale, which is closely related to the axion decay constant, commonly denoted by f a . This choice is suggested also by this scale appearing as the mass scale in the couplings of the dimension-5 derivative and two gauge boson interactions.
The adequacy of the EFT language to study axion phenomenology is based on the assumption that the presence of the axion is the only manifestation of the PQ construction at energies below the scale f a . At tree level, it is usually sufficient that any other particles are sufficiently heavy so that they cannot be produced at the energy scale relevant to the process under consideration. However, higher energy scales can become relevant at loop level where the contribution from off-shell particles in the loop needs to be included. In particular, when the EFT calculation yields divergent contributions, it may be more appropriate to identify the EFT cutoff scale with the mass of the lightest particle heavier than the axion instead of with f a .
In order to construct a suitable axion EFT it is therefore essential to check whether the underlying model features additional relevant degrees of freedom below f a . The goal of this paper is to do this in a number of full-fledged QCD axion models, thereby assessing the validity of studies performed using an EFT framework. As a particularly pertinent example, we focus on loop-induced flavour-violating decays involving the axion. Indeed, some of the one-loop diagrams involved in these processes are logarithmically divergent and therefore sensitive to the physics around the cutoff of the EFT. We show that the leading log prescription, where one focuses on the logarithmically divergent term and identifies the cutoff scale with f a , commonly employed in EFT calculations, produces results which can differ qualitatively and quantitatively from the ones obtained using full models. The reason for this discrepancy, as will become clear, is that most popular (DFSZ-and KSVZ-type) QCD axion models [94][95][96][97] do not satisfy the assumptions that implicitly enter the EFTbased loop calculation. Our results are in line with the earlier work in a similar direction presented in [87]. There, a leading-order RG evolution was employed to show that simple electroweak-complete axion models (e.g. DSFZ-type models explicitly including the two Higgs doublets) feature an absence of a large logarithm in the flavour-changing effects. Our loop-based calculations are in agreement with these results.
In this situation, the question arises whether there exists an explicit field theoretic 3 UVcomplete QCD axion model that allows for large logarithms in flavour-changing observables. A positive answer amounts to finding ways to generate tree-level couplings between the axion and SM fermions without introducing new states below the scale of PQ symmetry breaking. We address this in two steps. First, we consider an effective model valid up to a scale f a , which coincides with a specific charge assignment of the Lagrangian considered in ref. [87], where it was however not further explored. Our effective model also features some similarities with the flavoured axion models of [98,99]. Crucially however, the charge assignments in our case are not flavour-dependent and thus tree-level flavour-violating interactions are not present in our setup. In the next step and drawing inspiration from Froggatt-Nielsen models [100], we construct a UV completion where the logarithms are explicitly calculable. The result is a QCD axion model that features flavour-conserving SM fermion couplings without the need for an extended Higgs sector, which therefore provides an affirmative answer to the question posed above.
The logarithmic enhancement of the K + → π + + a decay rate in this model makes it a particularly good candidate to be tested at experiments like NA62 [101] and KOTO [102] as well as future experiments such as KLEVER [103]. 4 Furthermore, the potential existence of an axion-electron coupling increases the production of axions in the solar interior. This enhancement results in a larger expected axion flux in helioscopes like IAXO [11], which are projected to have sensitivity to the model in the O(10 meV) mass range.
In essence, the goal of the present study is to asses whether a given combination of axion EFT and cutoff is consistent with an embedding into a more fundamental theory. In the context of quantum gravity and string theory, the terms landscape and swampland have become popular to respectively denote the set of low energy effective theories that can and cannot consistently be incorporated into a theory of gravity [104]. Although our considerations are purely field-theoretic (see [105] for a swampland conjecture that is independent of the Planck scale), the landscape/swampland analogy serves as a motivation to raise the question of whether there are additional constraints on embeddable axion EFTs that have not yet been fully appreciated. In the case at hand, we find a promising avenue to construct the desired UV completion. That said, our model hints at some potentially non-trivial requisites. For example, achieving the absence of tree-level flavour-changing interactions in the EFT requires a non-trivial choice of the model parameters and potential tuning. Furthermore, realizing a sufficiently large top Yukawa coupling may constrain the separation of scales between the UV cutoff and f a , as was noted in [87]. This suggests that mild additional assumptions (e.g. the absence of fine-tuning) may give strong constraints on the possible embeddings of axion EFTs into full UV models. This paper is structured as follows. In section 2, we set up the QCD axion EFT under examination and calculate the expected rate for the K + → π + + a decay within this framework. This result is compared with detailed calculations in DFSZ-and KSVZtype QCD axion models in section 3. After establishing the discrepancy between the EFT prescription and existing models, in section 4 we present a new QCD axion model that satisfies the assumptions under which the EFT is constructed. The phenomenology of this model is explored in sections 5 and 6, after which we conclude in section 7. A discussion of CP violation in axion couplings and most of the technical details can be found in the appendices A-E.

Loop-induced rare decays in QCD axion EFT
The most general effective Lagrangian describing all possible interactions of an axion (or axion-like particle) with the Standard Model fields before electroweak (EW) symmetry breaking and involving dimension-5 operators can be written in a compact form as [54,106] where the first sum runs over SM gauge bosons with couplings c F F . In the second sum, the SM chiral fermion multiplets are summarized in χ = (Q L , L L , u R , d R , e R ) and each of the C χ is a matrix in generation space that allows for flavour changing effects (see appendix A for a discussion of the CP properties of C χ ). It is also customary to define the gauge boson couplings g aF F = c F F α F /(2πf a ), but for our purposes the dimensionless couplings c F F are more convenient. The dimensionful quantity f a , usually known as the axion decay constant, serves as the large scale for the expansion of the operators in the EFT. In Eq. (2.1), there is an ambiguity in the definition of f a and the EFT coefficients. As long as the coupling to gluons c gg is nonzero, as is necessarily the case for any QCD axion, we can solve this ambiguity by normalizing f a in such a way that c gg ≡ 1, and we do so in the rest of this work.
The other relevant feature of the Lagrangian is that a potential for the axion is generated. For a nonzero coupling to gluons, non-perturbative QCD dynamics provide for a mass term [92,95,107,108] At energies below the electroweak scale, it is convenient to rewrite the axion-fermion couplings in Eq. (2.1) using the mass eigenbasis for the SM fermions, where the sum runs over the up, down, charged lepton, and neutrino flavour-triplets u, d, , and ν. Each f thus summarizes the three generations. The coupling matrices c f,L and c f,R for left-and right-handed fields are related to the initial C χ matrices in Eq. (2.1) by unitary matrices [54]. In the left-handed quark sector, we have the additional relation  Figure 1. One-loop diagrams inducing the flavour-violating s → d + a transition. We include the self-energy contribution on the external strange-quark leg (analogously for the down-quark line) arising from the renormalization of quark fields, see appendix D. All Feynman diagrams throughout this work are drawn using the TikZ-Feynman [111] package.
coupling structures. While some axion models [98,99] feature flavour-dependent and even flavour non-diagonal interactions (see [71] for a summary of the associated phenomenology), they are not the main focus of this work. That said, as such couplings may also arise in our explicit model construction, we discuss them in some detail later. For now, let us consider flavour-universal couplings in the quark and lepton sectors. This is already sufficient to demonstrate our main point, namely that the leading logarithm in FCNCs obtained in the naive EFT does, in many cases, not agree with the result of more complete models. In this case, we can remove the vectorial part of the quark and lepton couplings by performing universal vectorial phase rotations of the quark and lepton fields. These leave the mass terms and the charged current interactions with W ± bosons invariant and therefore no scalar Yukawa couplings to a appear. The only possible effects are twogauge-boson couplings arising for the chirally coupled electroweak fields, as noted in [71], corresponding to the anomaly for the baryon and lepton number U(1) symmetries. The effect can therefore be absorbed in a redefinition of the corresponding two-gauge-boson coupling coefficient, most notably c W W . As we discuss below, the contribution from c W W to the loop processes studied in this work is finite and therefore a shift in this coefficient does not play a role when evaluating the leading logarithm. We thus only need to consider axial-vector couplings of the axion to quarks and leptons, 5 where the neutrino fields are understood to be purely left handed, ν = ν L . With these assumptions, no inter-generation fermion couplings and therefore no flavour-changing processes are present at tree level. But even if these are absent, flavour-changing neutral currents (FCNCs) appear at one loop through the diagrams depicted in figure 1. Rare decays induced by such transitions have been exploited in the literature to test the axion paradigm [59,60,63,64,[69][70][71], and lepton-flavour violation has been studied in similar detail [109,110].
In this work, we perform a detailed study of the FCNC transition in figure 1, with a focus on the impact of the possible UV completions of the axion effective Lagrangian in Eq. (2.4). For concreteness and also because it is of significant experimental relevance, we concentrate on the kaon decay K + → π + + a. That said, our results are also applicable to other FCNC-induced processes.
The s → d + a transition can be described using the effective Hamiltonian The hadronic matrix element for the K + → π + + a decay 6 can then be parameterized as [112] where P = p + p and q = p − p so that the momentum transfer is q 2 = m 2 a . For a sufficiently light axion, we only need the first form factor at q 2 0. Recent Lattice QCD evaluations [113] give f + (0) = 0.9706 (27). In this limit, the decay width can be expressed as where we have introduced the notation λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + xz + yz).
Within the framework of Eq. (2.4), the coefficients h S,P ds can be expressed in terms of the parameters of the axion EFT. Given that, as we have seen, the QCD matrix element for h P ds vanishes due to parity, we only need [60,114] where we have substituted the divergence arising from the calculation of the diagrams in figure 1 by a leading logarithm depending on an a priori undetermined UV scale Λ. In an EFT sense, this identification can also be understood by considering a leading-order RG evolution between the high and the low scale [87,89]. We have also introduced the loop function [60] f and expanded the amplitude at leading order in the external momenta (corresponding to dropping terms suppressed by m d /m W , m s /m W , and m a /m W ) in order to determine the leading finite contribution. This expansion, however, has barely any effects compared to the "leading log". Finally, we have neglected the contribution from the gluon anomalous coupling, which can be computed based on mixing effects [47,115] but is subleading to the leading logarithmic term. The result above has been computed in unitary gauge using the basis of Eq. (2.4), where the QCD axion couples only derivatively. An axion-dependent rotation of the fermion fields allows to trade the derivative terms for pseudoscalar ones ∂ µ af γ µ γ 5 f → 2m f af iγ 5 f , importantly along with anomaly terms. One can check that the pseudoscalar interactions give the same leading log result as the axial-vector ones. In Feynman gauge, a factor of 4 difference between the axial-vector and pseudoscalar interaction has been reported (see footnote 3 in ref. [67]), but we have confirmed that the claimed discrepancy disappears when the EW Goldstone bosons are taken into account.
Combining the above equations, one can evaluate the decay rate as a function of the EFT parameters and, importantly, of the UV cutoff Λ. This can in principle be used to place bounds on f a using the experimental constraint Br(K + → π + a) < 7.3·10 −11 from the E787 and E949 experiments at BNL [116], which NA62 expects to improve by an order of magnitude by 2025 [117,118]. In practice, the reach of the constraint depends significantly on the choice of the value of Λ. It is common to identify Λ with the only other scale that is available in the EFT description, which is the decay constant f a [63,67,114]. However, this is an ad hoc choice which, as we will see, is not reproduced in the most popular UV completions of the axion EFT. In order to clarify this issue, we now turn to study the s → d + a transition in the benchmark KSVZ and DFSZ QCD axion models.
3 Loop-induced rare decays in QCD axion models

DFSZ-type models
Let us compute the flavour-violating decay of the kaon through the diagrams of figure 1 in DFSZ [96,97] models. Crucially, we use a complete DFSZ model explicitly including the two Higgs doublets in order to regularize the logarithmic divergence and identify the cutoff scale Λ. The results are obtained from [59], where the B → K + a process was studied, by applying the obvious substitutions. Before going into the details, let us note that these calculations are performed in the basis where the axion interactions are of Yukawa and not of derivative type. As we briefly discuss below Eq. (3.8), this avoids complications arising when separate chiral rotations of up-and down-type quarks are performed.
Our starting point is the simplest DFSZ model, which is a two-Higgs-doublet model (2HDM) of Type II extended by a singlet complex scalar Φ that transmits its PQ charge via the operator The PQ symmetry is broken when Φ acquires a vacuum expectation value (VEV) f a , which we assume to be much larger than the EW VEV v = (v 2 u + v 2 d ) 1/2 . As usual, we denote tan β = v u /v d . The axion inherits couplings to the SM quarks due to its mixing with the pseudoscalar Higgs, resulting in axion Yukawa interactions with the up ∼ (m u,i /f a ) cos 2 β/3 and down ∼ (m d,i /f a ) sin 2 β/3 quarks. With this, the value of h S ds in   Figure 2. Expected NA62 reach on f a from K + → π + + a in flavour-universal DFSZ models, in terms of the 2HDM parameters m H ± and tan β. We use the projection of [117] for Br(K + → π + +a). the effective Hamiltonian Eq. (2.5) describing the s → d + a transition is found to be The sum runs over all the up-type quark flavours and the (finite) loop functions are given by [59] in terms of the charged Higgs boson mass m H ± . The different sign of the leading logarithmic term in the limit of large m H ± compared to eq. (2.8) is due to a negative c q in the DFSZ model. We can then use Eq. (2.7) to compute the expected K + → π + + a branching ratio in DFSZ models. The expected NA62 constraints [117] as a function of m H ± and tan β are shown in figure 2. Note that for each value of tan β, there is a value of m ± H for which h S ds in Eq. (3.2) changes sign, which causes the funnel in figure 2 along which the bound disappears.
With this computation, it becomes explicit that the value of the cutoff scale of the leading logarithm for the DFSZ model is the charged Higgs boson mass, Λ DFSZ = m H ± , which is orders of magnitude smaller than f a in the interesting region of parameter space. Moreover, cancellations in the decay rate can occur for specific parameter points, an effect that is not captured within the EFT approach. Finally, FCNCs can also be suppressed by large values of tan β [87]. Overall, this means that the K + → π + + a branching ratio in DFSZ models is much smaller than expected based on the EFT computation.
As a consequence, the actual bounds and projections from experiments such as NA62 on the DFSZ axion parameter become correspondingly weaker. This comparison can be seen in figure 5 in section 6, where the dark region represents the branching ratio for a benchmark value of m H ± = 800 GeV and 1 ≤ tan β ≤ 5, roughly representative of the experimentally allowed 2HDM parameter space (see e.g. [119,120]). We also display a lighter region which corresponds to only requiring perturbativity of the Yukawa couplings, 0.25 ≤ tan β ≤ 170 [121].

KSVZ-type models
The computation above shows that the EFT calculation with the naive leading-log approximation is not appropriate for DFSZ axion models. The reason for this is that DFSZ axion models introduce new degrees of freedom (the second Higgs doublet) at a scale below the naive EFT cutoff f a . In KSVZ [94,95] models, one expects that the EFT leading log term should not appear either. The reason is that in KSVZ models no tree-level fermion couplings are present, and even if they can be generated by an axion-dependent chiral rotation of fermion fields, the logarithmically divergent contribution from the 1-loop process involving the derivative coupling ∂ µ a/(2f a )f γ µ γ 5 f exactly cancels with the one coming from the pseudoscalar coupling ia/f a m ff γ 5 f , so that only a finite piece remains. From a different perspective, this cancellation reflects the equivalence of the linear and polar representation used for the complex scalar field containing the axion [122].
To check this statement, we start with and perform a chiral rotation [54] on all SM quarks, i.e. q = (u, d, s, c, b, t) T , with where κ q is a 6 × 6 diagonal matrix. Taking, e.g., Tr[κ q ] = 1 removes the gluon coupling completely [54]. As the one-loop process for s → d + a is O(g 2 ), we have to collect all terms that contribute at the same order. The relevant terms before the rotation are yielding after the chiral rotation (expanded at LO in the ALP field) where M q is the quark mass matrix and κ u and κ d are the 3 × 3 submatrices of κ q including only the elements for the up-type and down-type quarks, respectively. After tracing the influence of these terms on s → d + a at one loop, one notices that the UV-divergent terms cancel exactly. Note that in this context the third term in Eq. (3.8) provides crucial contributions by axion emissions off the W − q − q vertices, which do not vanish unless κ q is fully flavour-universal. 7 We can therefore robustly conclude that a leading log does not appear in KSVZ models at the one-loop level.
4 An EFT-inspired QCD axion model

Low-energy effective model
Because the naive EFT result is not recovered in common axion models like KSVZ or DFSZ, one wonders if there is any UV completion of the axion EFT where the leading-log result is actually applicable, meaning that the full result of the kaon decay width is also enhanced by as large a logarithm (or an even larger one) as in (2.8). From our discussion above, we can identify two key requirements for such a model. First, it should not contain any additional new states below the new-physics scale Λ apart from the axion (unlike the DFSZ model). And second, it should feature physical couplings to quarks that cannot be removed by a chiral rotation (unlike the KSVZ model).
A general class of effective QCD axion models that satisfy both of these requirements is given by the Lagrangian 8 where Q L , u R , d R , and H denote the usual SM chiral fermion fields and Higgs doublet before EW symmetry breaking, and the indices i, j label the three generations. The complex scalar Φ = φ/ √ 2 e ia/ φ contains the axion as its angular degree of freedom, and we have therefore normalized its PQ charge to χ Φ = 1. The Higgs doublet has no PQ charge 9 and the charges of the quarks therefore satisfy In order to avoid having to impose any restrictions on the coupling matrices Λ u,d , we assume that the PQ charges are flavour independent, χ L = χ Q Li and χ R = χ d Ri = χ u Ri for all generations. (4.2) Note that this is the essential difference to the models considered in [98,99], where the charges were assumed to be flavour-dependent in order to explain the observed flavour structure of the Standard Model Yukawa couplings. Such a choice generically leads to tree-level flavour-violating interactions for the axion. As our focus is on the loop-induced flavour effects, we choose flavour-independent charges to avoid these complications.
The QCD anomaly induced by the SM quarks is as in the usual DFSZ models. The electromagnetic anomaly coefficient, can either be 8 or 5 depending on whether the leptons are PQ-charged or not (leptons can also couple to Φ * , resulting in a negative contribution to E). Once the PQ scalar acquires . Because the axion is automatically aligned with the Yukawas, there are no tree-level flavour-violating axion couplings in this model.
After EW symmetry breaking and CKM diagonalization, the usual pseudoscalar couplings of the axion to quarks, are recovered, meaning that the couplings are strictly proportional to the fermion masses.
The presence of axion-quark couplings in this model implies the existence of a divergent contribution to flavour-changing amplitudes at one loop from diagrams as the ones shown in figure 1. When replacing this divergence by a leading logarithm, the only possibilities are log( φ 2 /m 2 q ) or log(Λ ij /m 2 q ). Since all of these scales are at least as large as the PQ breaking scale, we therefore necessarily obtain the large logarithmic enhancement of flavour-changing processes anticipated from the axion EFT.
Since we started from a non-renormalizable Lagrangian in Eq. (4.1), it is not possible to obtain an exact result for the leading logarithm and as such one also lacks its physical interpretation. It is therefore instructive to understand how the effective interactions discussed above can be embedded in a UV-complete model, which allows to calculate the flavour-changing processes in detail.

UV completion of the model
In order to embed the effective model introduced above into a renormalizable UV completion, we add three generations of heavy coloured up-and down-type fermions, F u i and F d i , which resemble the messenger fields in the Froggatt-Nielsen mechanism [100]. For simplicity, we first study the case of a single generation, and later generalize it to three generations.
For one generation, the relevant renormalizable Lagrangian involves the terms where α u,d and β u,d denote the couplings in the up-and down-type sectors. 10 More generally, one can allow the possibility to use either Φ or Φ * . This leads to slightly different models, which are described in appendix B. We focus on the example of the model as specified in Eq. (4.7), but the other options can be treated in a similar way. For all terms to be gauge invariant, F u and F d must transform under the SM gauge groups like u R and d R , respectively, and left-and right-handed F fields (we omit the u/d superscript when referring to any of the two fields) must transform identically. The terms in Eq. (4.7) always enforce the conditions χ Q L − χ F R = 0 and χ F L − χ q R = 1. It is precisely the sum of these equations that appears in the anomaly coefficients N and E in Eqs. (4.3) and (4.4), as up-and down-type F and SM quarks couple similarly to gluons and photons. Hence, the additional SU(3) quarks do not change the anomaly coefficients N = 3 and E = 8 or 5. Furthermore, a bare mass term for the F quarks, can be added. For easier extension to the multi-generation case, we have parameterized the mass in terms of the VEV of a real spurion field Σ for which χ Σ = 0. This imposes the additional condition χ F L − χ F R = 0. With this, all axial charges of all combinations of SM and F quarks have been fixed and the only remaining freedom is a shift of all charges by an arbitrary constant. It is clear that after PQ and EW symmetry breaking, the terms in Eq. (4.7) induce mass mixing between the different fields involved. In unitary gauge, we can write the mixing in matrix form as where M depends on the VEVs of all scalar fields in our theory (H, Φ, Σ), the corresponding coupling matrices (α, β, λ), and the axion field a.

Light-and heavy-quark mass matrix diagonalization
In order to recover the quark masses and couplings, we proceed exactly as in the Standard Model and diagonalize the mass matrices M u and M d by unitary transformations of leftand right-handed fields. Only the main steps and results are mentioned here: we refer to appendix C for detailed calculations. Throughout this and the next subsection, we drop the labels u/d whenever the procedure is identical for both types of quarks. We consider at every point all three generations of SM and F quarks even though the generation indices are also omitted.
The full mass matrix M is given by where we have introduced two parameters with the Higgs VEV v = 246 GeV, which highlight the hierarchies of scales in our model. We can safely assume that the parameter is much smaller than one. In contrast, may but does not need to be smaller than one for a general axion model, where extra degrees of freedom other than Φ may appear close to the PQ scale. Importantly, such additional degrees of freedom close to the PQ scale can lead to sizable corrections to our desired effective model in Eq. (4.5), introducing e.g. tree-level flavour-violating axion couplings.
In order to avoid this complication, we assume that the additional degrees of freedom at hand (the F quarks) are heavy enough so that these effects are sufficiently suppressed. In practice, this means that we require where min m F i is the mass of the lightest fermion that is not part of the SM. More details, such as quantifying how small has to be in the most general case as well as how tree-level flavour violation can be avoided even for sizeable , can be found in appendix C.
Continuing with the diagonalization, we define U to be a unitary matrix which diagonalizes the hermitian product M M † , (4.13) Here, M q and M F are the diagonal mass matrices of Q and F quarks. Transforming the left-handed quark fields with U and the right-handed ones with exactly gives us a basis of fields in which the mass matrix is diagonal and independent of a. Note that to lighten the notation we are using identical symbols for fields in the original Lagrangian, Eq. (4.7), and for the mass eigenstates after the diagonalization.
By expanding U in and , we can express the quark masses in terms of the UV parameters. At leading order, the masses are where U δ and U ξ are unitary matrices which diagonalize the hermitian matrices in parentheses. It is important to keep in mind that all matrices in (4.15) and (4.16) exist for upand down-type particles, giving a total of 12 different quark masses. Any physical realization of the coupling matrices needs to reproduce the masses of SM quarks in Eq. (4.15).
Since the mass of the top quark is close to H = v/ √ 2, it may be difficult or impossible to recover such a high value when all F quarks are much heavier than the PQ scale and when simultaneously requiring perturbativity of all Yukawa couplings. A similar observation was made in [87], where it was concluded that the cutoff of a Lagrangian as the one in (4.1) can at most be of order f a . However, in our explicit UV completion the masses of the six F quarks can have a large intrinsic hierarchy, which can in principle lead to logarithmic enhancements of flavour-violating effects by factors larger than log(f 2 a /m 2 q ) and possibly also tree-level flavour violation, and simultaneously reproduce the observed top Yukawa. These effects are part of the non-trivial flavour properties of our model, whose detailed investigation goes beyond the scope of this work. In the main part of this text, we work in the limit 1 while remaining agnostic about the exact flavour structure that may realize this hierarchy. As a proof of principle, in appendix C we give an explicit example with = 0.2 for which all of our calculations apply and the SM is reproduced.

Axion couplings to gauge bosons
The matrices U and S introduced in the previous section generally contain axion-dependent phases. As a consequence, the quark-field transformations in Eq. (4.14), which are required to diagonalize the mass matrix, source axion interaction terms. To understand these better, it is helpful to split the diagonalization procedure into two subsequent field redefinitions, of which only the first one depends on the axion field.
Looking at the full UV Lagrangian in (4.7), it is clear that the a dependence of the mass-mixing terms can be absorbed into the quark fields by the redefinitions This transformation removes the axion field a from M in Eq. (4.10), while the quark kinetic terms generate derivative couplings of the axion to right-handed quarks. Hence, both diagonalization matrices U and S can subsequently be chosen to be independent of the axion field. The resulting derivative interactions are explicitly written out below in terms of the mass-diagonal quark fields (see Eq. (4.22)).
In addition to the kinetic terms, the path integral measure is not invariant under the transformation in (4.17) and anomalous interaction terms between the axion and gauge bosons arise. The U (1) A × G × G anomaly, with G being either the strong or the EW gauge group, sources the interaction terms [90,91,[123][124][125] where s W and c W denote the sine and cosine of the Weinberg angle θ W . Note that there are no anomalous W -couplings because we are only transforming right-handed fields, which are singlets under SU (2). From here, it is customary to define 19) where N DW ≡ 2N = 6 is the domain-wall number counting the number of inequivalent vacua in the QCD-induced axion potential. The consequences of N DW = 1 in our model are briefly discussed in section 7. The definition of f a results in the conventional normalization of the axion couplings,

Axion couplings to light and heavy quarks
After absorbing the axionic phase of the mass matrix in the q R fields, M becomes independent of a. Hence, U and S can also be chosen to be independent of a when proceeding with the subsequent steps in the diagonalization as described in Eq. (4.14). After the diagonalization is completed, the derivative axion-quark couplings induced by the transformation Eq. (4.17) can be expressed in terms of the mass-diagonal fields as where A and B are coupling matrices which, to leading order in , are given by Note that the reason why the axion only couples to right-handed quarks is that we chose to only rotate these chiral components in Eq. (4.17). This is clearly an arbitrary choice: we can perform an axion-dependent vector (non-axial) rotation in the quark field to include derivative couplings to left-handed quarks.
The field basis in which axions only couple derivatively is convenient for the calculation of loop processes as the one in section 5.1. But if we want to recover the effective terms in Eq. (4.6), we have to rewrite the first term in Eq. (4.22) by performing a rotation of the right-handed fields (neglecting interactions with the Higgs and CP-conserving anomalous axion terms), Hence, we exactly arrive at purely pseudoscalar couplings to SM quarks which are proportional to the quark masses, as anticipated in Eq. (4.6). Note that this implies that the right-handed flavour-diagonal coupling structure conserves CP, which can equally be shown by applying the CP transformation explicitly, see appendix A.
In contrast, applying the same steps to the second term in Eq.
This demonstrates that the coupling of heavy F -quarks is not proportional to their masses, which is due to the fact that their mass is (mostly) generated by the spurion field Σ, whose VEV does not break the PQ symmetry. Furthermore, the interaction with F quarks is not of purely pseudoscalar form, which, as is discussed in section 5.2, has important implications for CP violation in our model. Finally, for the last term in (4.22) we make use of the fact that q R and F R only differ in their respective masses but transform identically under the SM gauge group, and therefore where the indices i and j do not imply a sum but only a specific quark. This allows us to rewrite the last term in (4.22) as Once more, the interaction contains both pseudoscalar and scalar parts which can in general be CP violating.

Low-energy effective couplings
To conclude this section, let us summarize how the model introduced above maps onto the generic axion EFT. In this context, it is important to note that the low-energy effective model does not contain any free parameter other than the axion decay constant and the mass of the heavy F -fermions, the latter only being relevant as a cut-off scale for loop processes (see the next section). All the couplings of the axion to SM particles are therefore fully determined once f a is fixed. There is only one discrete choice, namely whether the SM leptons are charged under the PQ symmetry or not, leading to the two variations of the model. Table 1 summarizes the values that the EFT coefficients introduced in Eqs. (2.1) and (2.3) take for both variations.  .3). The general dependence on the high-energy parameters is shown, together with the explicit value for each of the two variations of the model, the one with and the one without tree-level couplings to leptons. In the three cases c γγ , c γZ , and c ZZ , the coupling term should be written as in (4.20) with the electromagnetic fine structure constant α EM . Further, the fermion coupling matrices in the last two rows are proportional to the unit matrix and 1/6 − 0 indicates that only right-handed fields couple at tree-level.

Flavour and CP effects in the EFT-inspired QCD axion model
In the previous section, we have introduced a QCD axion model with tree-level couplings to SM fermions and without any additional new states below the scale f a other than the axion itself. As a consequence, the model fully reproduces the phenomenological expectations based on the generic EFT description introduced in section 2, even when one-loop processes are involved. To confirm this, in this section, we delineate the main phenomenological features of this QCD axion, including a calculation of the K + → π + + a decay rate in the full model.
In the previous section, we have computed the effective interactions between the axion and the SM quarks, which are summarized in  amplitude at one loop and we would again be confronted with the same problem of UV divergences as before. Crucially, however, there are indeed further relevant interactions arising in our UV model that are not contained in table 1 and which render the loop computation finite. These operators are, as expected, related to the heavy F quarks and are hence not captured by the EFT approach.
To obtain these additional vertices, it is important to discuss further implications of the unitary transformation of the quark fields that we performed in Eq. (4.14). As in the SM, this rotation generates new flavour-dependent interactions of neutral and charged hadronic currents with the weak gauge bosons. The detailed calculations are included in appendix C. The leading-order W interactions in the mass-diagonal basis can be written as where we have reintroduced the labels for up-and down-type quarks as well as coupling matrices, and V denotes the CKM matrix as before. At leading order, V is given by V = U u † δ U d δ and turns out to be unitary up to terms of the order ( ) 2 . Similarly, the leading-order Z interactions become 2) Here, the upper (lower) sign refers to up-(down-) type quarks, P L denotes the projector onto left-handed fields and Q is the electromagnetic charge of each field. We find that only the left-handed coupling structure is affected by the transformation to the mass-diagonal basis.
With these additional interactions at hand, we can move forward to perform the computation of the loop processes shown in figure 3. The first diagram in figure 3 represents the sum of four individual contributions because each internal fermion propagator can either be a SM or an F quark. To calculate these diagrams, we make use of the axion-quark interactions as given in Eq. (4.21) and the W couplings to hadronic currents in Eq. (5.1). Expressing the interactions in this way means that the W boson and the axion only interact with left-handed and right-handed fields, respectively. Hence, we need mass insertions by each of the internal propagators. The corresponding factors of fermion masses are however cancelled by the inverse Λ matrices in the axion interaction. Putting everything together, it is easy to see that all four diagrams have the same flavour structure and are of the same order in and . We can write their combined contribution to the s → d + a amplitude in the compact form where m Q refers to the up-type quarks. We evaluate the loop integral using Package-X [126] and under the simplifying assumption that all F quarks have equal masses. Only considering the leading-order terms in the down and strange fermion masses, we arrive at where we expanded the finite contributions in terms of 1/m F . This gives a contribution to h S ds = h ds as defined in the effective interaction Hamiltonian in (2.5). The logarithmically enhanced term is identical to the EFT result for c q = 1 6 and Λ = m F . Note that the finite term is different from the EFT result as the loops with internal F -quarks also yield relevant contributions.
In principle, we also need to include diagrams exchanging the W boson in the first diagram in figure 3 with a Z boson. However, even if the couplings of the Z boson do have non-diagonal entries, they do not induce s → d + a at one-loop at the order in that we are considering. This is due to the matrix structure ABA † appearing in the relevant amplitudes, which is flavour-diagonal as it is identical to the matrix that determines the quark masses in Eq. (4.15).
Finally, there are also counterterm contributions [127] from the renormalization of quark fields as depicted on the right of figure 3. Performing the calculations, however, one notices that the relevant contributions to Eq. (5.4) cancel each other at linear order in the down-and strange-quark masses (see appendix D for more details). We can therefore also discard these counterterm diagrams and work with Eq. (5.5) as our effective coefficient.

CP violation in axion interactions
From Eq. (4.22), we observe that the axion has a right-handed coupling structure to all quarks in our theory. As discussed in appendix A, CP violation in flavour-violating couplings occurs when the corresponding coupling is complex valued. The axion interactions in Eq. (4.22) can therefore induce CP violation since A and/or B can have imaginary entries. The fact that the interactions of the QCD axion can be CP violating may seem worrisome at first, since the QCD axion is introduced precisely to eliminate the CP violation in the strong sector. As shown in this section, however, the effects of this kind of CP violation on the electric dipole moment (EDM) of the neutron are comparable or even smaller to those induced by the phase of the CKM matrix, and are therefore not a threat to the axion solution of the strong CP problem. We point the reader to refs. [128,129] for recent discussions on CP-violating axion interactions.
To see this, firstly note that this CP violation is of explicit nature, like the CKM phase in the SM, and not due to any spontaneous breaking. This can be understood considering the full UV theory in Eq. (4.7). If the couplings α, β, and λ are real-valued, no CP-violating terms are generated before or after EW and PQ symmetry breaking. This is because in that case, U δ and U ξ appearing in Eqs. (4.15) and (4.16) can be chosen to be orthogonal and real because the matrices that they diagonalize are real and symmetric to begin with. This then results in A = U T δ αλ −1 U ξ and B = U T ξ ββ T U ξ being real-valued as well. Hence, Eq. (4.22) is CP conserving in this case. However, starting from complex-valued α, β, or λ couplings generically makes A and/or B complex valued as well. Therefore, the CP violation present in the axion interactions of Eq. (4.22) is a reflection of potential explicit CP-breaking in the UV model. As a matter of fact, given that the CKM matrix, which is constructed from α, β, and λ, is experimentally confirmed to be complex valued, there is no reason to expect the coefficients A and B to be real without invoking an ad hoc cancellation.
The relevant question is therefore whether this CP violation is problematic, that is, whether it is in conflict with any observations. To answer this, we consider a very sensitive observable to CP violation: the neutron EDM d n .
As shown in appendix E, one-loop processes do not generate quark EDMs in our UV model. Hence, a neutron EDM based on free quark EDMs may occur at the earliest at the two-loop level. Without delving into involved two-loop calculations, we instead derive a very conservative upper bound based on the expected scaling of the contributions. Assuming that at least two new physics vertices appear (such as internal emission and absorption of the axion), each of which carries a coupling suppression of f a , we get the rough estimate (5.6) Here, m n denotes the neutron mass, which we use as the characteristic scale, and we have inserted f a = 4 · 10 6 GeV as a lower-end value for the range of axion decay constants of interest. The above estimate is already very conservative in the sense that we have not included any suppression due to new physics contributions to W ± and/or Z gauge boson interactions and/or the heavy F -quark mass scale Σ , since at least one of the three has to participate to generate a quark EDM. Moreover, we have not invoked any additional electric/weak coupling insertion or Yukawa/light quark mass suppression, the latter being expected to appear given the discussion above regarding the origin of the CP-violating interactions.
Interactions between the quark constituents within the neutron have also been shown to contribute to the neutron EDM (see ref. [130] for a review). That said, these effects are already below present sensitivity for purely SM-related interactions. This does not change in our UV model, as the new physics couplings/particles appearing would only lead to a stronger suppression. We can therefore safely neglect such contributions.
All in all and despite being very conservative, the estimate in Eq. (5.6) is still much smaller than the current bounds [131]  6 Phenomenology of the model and discovery opportunities Because our QCD axion model was constructed to reproduce the general features of the EFT setup, one may be sceptical that the model could posses any particular phenomenological feature that may serve as a handle to probe it. But interestingly, what at first seems like a lack of features does in fact have interesting consequences for experimental searches, which we now describe.

Astrophysical limits and helioscope searches
As any other QCD axion model, the one presented in this work is subject to constraints coming from its couplings to photons, electrons and nucleons. Figure 4 summarizes the bounds as a function of f a and in the usual m a -g aγγ parameter space.
The axion-photon coupling is tested by a combination of laboratory experiments and astrophysical observations. The negative results from the solar axion search performed by CAST [10], together with the absence of exotic cooling in Horizontal Branch (HB)  Figure 4. Constraints on the two variations of the model introduced in section 4, which we denote EFT-and EFT depending on whether tree-level axion-lepton couplings are present or not. On the left, we present the bounds as a function of the QCD axion decay constant. For the K + → π + +a constraints, the lighter shading of orange represents the range of possible values of the cutoff scale for the leading log, Λ ∈ (30f a , M Pl ) as discussed in the main text. The vertical hatching indicates uncertainties in the supernova limit. On the right, we show the models predictions in the usual mass vs. photon coupling parameter space. The band of QCD axion models corresponds to E/N ∈ (44/3, 5/3) as defined in [132]. Note that IAXO becomes particularly sensitive to the EFTmodel due to solar axion production through axion-electron interactions. Hence, the predicted IAXO sensitivity in dark green is based on the interactions with both electrons and photons. stars in globular clusters [133], constrain the photon coupling to be g aγγ 7 × 10 −11 . This translates into a bound on f a which is different for the leptonic and the non-leptonic models due to their different values of E/N (note that in the leptonic case, this bound only considers the axion flux generated by interactions with photons and not electrons; it is therefore overly conservative and only illustrative). Furthermore, the proposed successor of CAST, the IAXO solar telescope [11], is expected to improve this limit by over an order of magnitude in its upgraded version IAXO+. Unfortunately, the loss of sensitivity for axion masses above ∼ 10 −2 eV, combined with an accidental cancellation 11 in the photon coupling, renders the non-leptonic model out of reach for IAXO+, as can be seen in the right panel of figure 4. The model with tree-level couplings to SM leptons, however, is within reach of the sensitivity forecast, as can be appreciated in the left panel of figure 4. The reason for this is the enhanced production of axions in the solar interior caused by the existence of a large electron coupling. Finally, haloscopes searches have an exquisite sensitivity to axions in the 1 − 100 µeV range, provided these particles were to make up the dark matter of the Universe. The currently excluded region, as compiled in [30,69] and including new ADMX data [28,29], is shaded in yellow in the right panel of Fig. 4.
Axions coupling to electrons induce additional cooling mechanisms in stars. Currently, the strongest limits come from observations of the brightness of the tip of the red-giant branch (RGB) in globular clusters, which exclude values of f a /c e ≥ 3.9 × 10 9 GeV, as derived in [134]. 12 The axion-electron coupling is present at tree level in the leptonic 11 The low-energy QCD axion-photon coupling is a combination of the model-dependent anomaly coefficient E/N and a model-independent piece coming from hadronic contributions [107]. 12 Slightly weaker limits come from measurements of the R parameter in globular clusters [135,136] and variant of the model, and only loop induced (and therefore suppressed) in the non-leptonic one. As a consequence, the previously mentioned bound is only competitive for the former variation of the model, and suppressed by a loop factor [121,139] ∼ α 2 EM /(π 2 ) log ∼ 10 −4 in the latter. Indeed, in the presence of tree-level leptonic couplings, RGB observations currently place the strongest constraints on the model. As mentioned above, these are only expected to be improved once IAXO has reached its full sensitivity.
Finally, we discuss the constraints arising from the effective coupling of axions to nucleons, which arises at low energies from the interactions with gluons and quarks. It is customary to define it in a way analogous to the other fermionic couplings as The expression for the coefficients c N in terms of the quark couplings can be found in [107], which shows good agreement with a recent reevaluation [140]. For the model at hand, they are c p = −0.39875 and c n = 0.05125. These couplings are best tested by studying their impact on the extreme dynamics of the proto-neutron star that forms in the course of a core-collapse supernova. Building on the seminal reference [12], the most recent limits are given in [83][84][85]. While there is overall agreement given the uncertainties, the bound of [84] on f a is a factor of ∼ 3.5 stronger than the most conservative ones in [83,85]. In figure 4, we choose to use vertical hatching to showcase the spread in the different evaluations of the limit.

Constraints and opportunities from K
Perhaps the most important feature of the model at hand is that quark-flavour violating transitions q → q + a, which are induced at the one-loop level, are enhanced by a large logarithm as discussed in the previous sections. This means that processes involving such transitions are particularly suitable for testing this QCD axion model. Among all such processes, the s → d + a decay is the one that offers the best experimental perspectives, thanks to the very precise measurements that can be performed at kaon facilities. The predicted decay rate for K + → π + + a is given in Eq. (5.4), where we can see that the leading-log piece matches the EFT result with the identification Λ = m F . The corresponding branching ratio thus depends on the value of m F . In principle, the only requirement for the mass of the heavy F -fermions is for them to be above the axion decay constant, m F ≥ f a . However, given that SM Yukawa couplings are proportional to φ /m F , it is possible to obtain some rough limits on the largest possible values of m F that reproduce the SM quark masses. Assuming perturbativity of all dimensionless couplings, reproducing the top quark Yukawa requires at least one of the F fermions to have a mass m F j 10 3 f a for some j. However, the other F fermions could be much heavier, and without a complete construction of the 3-dimensional coupling matrices α, β, and γ, it is not possible to point from observations of white dwarfs, which exclude values of fa/ce ≥ 1.9 × 10 9 GeV, as derived in [137], although the data seems to prefer some amount of extra cooling compatible with a non-vanishing electron (and possibly photon) coupling [135,138]. down the contribution of each of the F fermions to s → d + a and therefore to the loop cutoff Λ. Therefore, and although it is reasonable to expect that Λ 10 3 f a , one cannot rule out the possibility of larger values.
Taking into account the discussion above, we can derive the bounds and prospects for this model in past and present kaon facilities. The K + → π + + a branching ratio as a function of f a for different values of Λ corresponds to the orange lines in figure 5. A robust exclusion can be placed by assuming a conservative value for the cutoff, Λ = 30 f a . Note that this choice is equivalent to 0.2 and hence corresponds to a slight scale separation between φ and Σ . This choice may appear in tension with the assumption 1 employed in our discussion so far. However, as shown in appendix C for an explicit choice of coupling matrices, our expression for the relevant effective s → d + a coefficient derived in eq. (5.5) still holds to a very good degree in this case. With this, the E787 and E949 experiments rule out f a < 1.6 × 10 6 GeV, and NA62 could push this up to f a < 5.6 × 10 6 GeV (left panel of figure 4, dark orange) assuming an order of magnitude better sensitivity in the branching ratio. That said, given that Λ could be much larger, the discovery potential for these experiments extends somewhat further. Only in this way, by identifying m F as the scale appearing in the logarithm and treating it as an independent parameter from the PQ breaking scale φ , can the NA62 sensitivity go beyond the robust limits set by the cooling of HB stars for the EFT scenario. For the most extreme value of Λ = M pl , the potential reach of NA62 extends to f a 1.6 × 10 7 GeV, which as can be seen in the left panel of figure 4 (light orange) comes closer to the SN1987a limit.

Conclusions and Discussion
Effective Field Theory tools are a cornerstone of studies of the QCD axion and general axion-like particles. However, the validity of an EFT is limited by the existence of a cutoff. This is evident in log-divergent loop processes, where the cutoff explicitly appears in the final result. It is tempting to infer the cutoff scale from the scale set by the dimensionful coupling constants of the higher dimensional operators in the axion EFT. Yet, one may wonder if this choice is validated by embeddings into simple axion models with renormalizable couplings. In this work, we have compared calculations made using an EFT setup with those obtained using complete QCD axion models. Our results show that matching EFT predictions to full models is far from trivial when loop processes are involved.
Owing to its experimental relevance, we have chosen to focus on the rare decay K + → π + +a for our study. We have shown that the K + → π + +a decay rate in the most popular QCD axion models can differ by orders of magnitude from naive EFT predictions, similar to what was found in [87] based on effective models that respect electroweak symmetry. The origin of this discrepancy lies in the UV sensitivity that is induced by the logarithmic divergence of this loop-induced decay. For DFSZ models, it is the existence of new degrees of freedom much below the PQ symmetry-breaking scale which violates the assumptions on which the EFT loop calculation is based. On the other hand, for KSVZ (or hadronic) QCD axion models, which do not have extra degrees of freedom below the PQ scale, suitable tree-level couplings to SM fermions are absent.
The above considerations appear to question the existence of a complete axion model which can successfully reproduce the EFT result, i.e. the large logarithmic enhancement found in the EFT calculation for loop-induced rare decays involving an axion. We have addressed this issue by explicitly constructing such a model, whose low-energy dynamics are dictated by the Lagrangian in Eq. (4.1). Crucially, this new model allows for tree-level axion-SM fermion couplings without introducing any extra degree of freedom below the PQ symmetry-breaking scale other than the axion itself.
While very minimal at energies below f a , the EFT-inspired QCD axion model constructed in this work presents rich dynamics above the PQ scale. In addition to the PQ complex scalar, three generations of heavy coloured fermions have to be introduced in order to obtain a renormalizable model. Through mixing with SM quarks, these heavy fermions induce modifications of the low-energy axion couplings. Ultimately, it is the presence of these additional states which tames the logarithmic divergence of the K + → π + + a rate, therefore providing a physically-motivated cutoff as opposed to the artificial one that is employed in the EFT.
The implementation that we have followed in this work is only one of the possibilities to UV complete the low-energy Lagrangian Eq. (4.1). Some other scenarios are delineated in appendix B. These may feature qualitatively different dynamics for the interactions of the QCD axion with the SM fields. One example that is briefly described in appendix C is the alternative that the mass of the heavy quarks lies close (or even below) the PQbreaking scale. In this case, the QCD axion generically enjoys tree-level flavour-violating interactions with the SM quarks, which significantly change the flavour phenomenology (see, e.g. [71]). Other interesting possibilities, like potential CP-violating axion interactions or connections between the axion dynamics and the flavour puzzle are interesting avenues of further research. Finally, the large logarithms appearing together with the top-quark Yukawa suggests upgrading our loop calculations with RG improvements along the lines of [87][88][89] to achieve better accuracy in the phenomenological limits.
The EFT-inspired QCD axion model has a rich phenomenology, the main features of which are summarized in figure 4. As is also the case in other axion models, astrophysical observables currently place the strongest constraints on axions with sub-eV masses. In particular, the duration of the neutrino burst in the supernova 1987a excludes decay constants f a 10 8 GeV. However, owing to the logarithmic enhancement of the K + → π + + a rate, K + → π + + inv searches at the NA62 experiment have the potential to supersede all astrophysical bounds except for the SN1987a one. Given the uncertainties in the modelling of axion emission in supernovae, 13 a confirmation of the exclusion from a terrestrial experiment is highly desirable. Additionally, the presence of additional colour-charged fermions with PQ charge could act to decrease the value of the QCD anomaly coefficient N , in which case the SN1987a limit would be weakened compared to the NA62 one.
The possibility to couple the axion to the SM leptons in a similar way as we have done for quarks offers even more exciting possibilities. Due to the particular balance between the photon and the electron couplings in our model, the IAXO helioscope (in its upgraded version IAXO+) is expected to become the most sensitive probe of sub-eV mass QCD axions, with a reach close to 10 9 GeV in f a (corresponding to a mass as small as 0.01 eV). This would supersede all astrophysical limits, including the previously mentioned SN1987a one and the one derived from observations of the brightness of the tip of the red giant branch in globular clusters.
The dark matter phenomenology of the investigated models is similar to those in the DFSZ and KSVZ benchmark models, with two notable differences. One is that the photon coupling can be quite small due to an electromagnetic anomaly coefficient of 5/3. Second, in the models that we have investigated the domain-wall number is greater than one, favouring the scenario where the PQ symmetry is broken during inflation and not restored afterwards. Alternatively, one may have an additional explicit breaking of the discrete Z N symmetry. While this can be useful to allow for the correct dark matter abundance at smaller values of the PQ scale f a [144], potentially accessible with experiments such as NA62, it may also lead to a certain degree of tuning.
The findings of this work showcase the necessity of developing a detailed understanding of the strengths and limitations of EFT tools in axion studies. We have seen that it is far from trivial to construct a UV embedding of an axion EFT coupled to fermions that agrees with the leading log predictions for relevant experimental processes. More generally, this leads to the question of whether a generic axion EFT with a cutoff given by the (inverse of) the dimension-5 coupling constants can be obtained from a UV model. While our explicit construction answers this question in an affirmative manner for the specific case at hand, it also highlights potential issues such as the quite complicated structure of the model, the potential appearance of additional flavour-changing couplings as well as potentially non-perturbative Yukawas being required to reproduce the top-quark mass. 14 It therefore remains an intriguing question whether even mild extra assumptions may severely restrict the space of axions EFTs that are the low-energy manifestations of reasonable UV models, possibly giving rise to an ALPine swampland.

A CP violation in pseudoscalar interactions
In this appendix, we investigate under what conditions CP violation arises in the couplings studied in this work. In particular, we focus on axion-quark, W -quark, and Z-quark interactions, where quark generally refers to a quark of any type, i.e. SM or BSM. The flavour-diagonal structure of photon and gluon interactions is not altered in our UV model compared to the SM and hence does not require any further comments.

Axion-quark interactions
We begin with derivative axion-quark interactions, differentiating between flavour-violating and flavour-conserving operators. The latter can in full generality be written as Hermiticity of this operator restricts h S and h P to be real-valued. A CP transformation leaves (A.1) unchanged and hence this interaction does never violate CP. On the other hand, the flavour-violating version of this operator and its hermitian conjugate read Applying the CP operation on the first term, we obtain and hence CP violation is closely related to the imaginary parts of h S and h P . Importantly, CP violation does not arise even if the axion simultaneously enjoys scalar and pseudoscalar couplings in the derivative basis, as long as both coupling constants remain real. For 14 Problems in obtaining the top-quark mass while retaining perturbativity were already noted in [87].
example, in the EFT scenario in Eqs. (2.5) and (2.8), the CP violation is fully determined by the imaginary part of the CKM matrix elements.

Gauge boson-quark interactions
We start with the Z bosons-quark couplings. For flavour-conserving interactions, we havē where hermiticity requires both couplings V and A to be real-valued. As a consequence, invariance under a CP transformation is automatically guaranteed. In contrast, flavourviolating couplings involvē A CP transformation on the first operator gives and hence CP violation is again closely related to complex-valued couplings, i.e. phases of V and/or A. In our case the coupling structure for flavour-violating interactions of the Z bosons is always left-handed (see Eq. (C.47)) and hence we have V = −A. CP violation therefore also amounts here to a global phase attached to the operatorqγ µ Z µ P L q .
The situation of W bosons is simpler given that in general only left-handed fields couple to it. The coupling structure are always of the typē As is familiar from the CKM matrix in the SM, CP violation manifests itself when (V − A) is complex-valued.

B Variations of the model
Instead of the specific UV completion in Eq. (4.7), we could have allowed either or both of the up-and down-type fields to couple to Φ * instead of Φ. In total, there are four possible ways to choose between Φ and Φ * . The most general Lagrangian reads where the upper sign corresponds to Φ and the lower sign to Φ * in the Lagrangian (B.1). Because all q and F quarks are in the fundamental representation of SU(3), the sum of the four equations above appears in the anomaly coefficient, Hence, |N | = 3 as long as we choose either Φ or Φ * for both up-and down-type quarks. Otherwise, the two contributions cancel, there is no QCD anomaly and a is not a QCD axion. We therefore do not further pursue that last possibility. Therefore, PQ charges of up-and down-type fields are necessarily identical, and we drop the corresponding labels in what follows.
After this, the only other choice is the PQ charge assignment of the F quarks. This determines the possible origin of the F -quark masses, which can either originate from the PQ-scalar Φ or from the VEV of an additional spurion field Σ without PQ charge. The combination of the two possible charge assignments of q and F quarks results in four different models: 1. Use Φ in (B.1) and generate F masses from Σ . The additional mass term is F is vector-like with respect to the PQ symmetry and does not contribute to the anomaly. This is exactly the model which is analyzed in great detail in sections 4 and 5.
2. Use Φ * in (B.1) and generate F masses from Σ . Using the same mass term as in (B.6), we get This option is of course equivalent to the first one after a redefinition of Φ ↔ Φ * .
3. Use Φ in (B.1) and generate F masses from Φ . The mass term of F quarks in this case is given by And consequently, Only F carries an axial PQ charge, which means that no symmetry forbids a Yukawa term as in the SM, The existence of these couplings slightly modifies the diagonalization of the mass matrix. Even though the q fields are not charged under PQ in the UV model, they inherit axion couplings through the mixing with the heavy quarks during the diagonalization. These are however parametrically suppressed by 2 = v 2 φ 2 .
4. Use Φ * in (B.1) and generate F masses from Φ . Using the same mass term as in (B.9), we get Both q and F quarks have axial PQ charges and therefore tree-level axion couplings in this variation of the model.
In each of the cases above, all axial charges of all combinations of SM and F quarks are fixed. The only remaining freedom is a shift of all charges by an arbitrary constant, which can be used to fix one vector-like charge.
Note that it only makes sense to work with the effective Lagrangian Eq. (4.1) in the first (or the equivalent second model) when φ m F i . Otherwise, the F fields cannot be integrated out at any scale above PQ symmetry breaking as is done in the main text. In models 3 and 4, a full diagonalization including all q and F quarks along the lines of what is done in appendix C is always necessary.

C Mass diagonalization and axion interactions
This appendix contains the full diagonalization procedure of the quark mass matrix as well as the resulting couplings of the axion and light and heavy quarks. For convenience, this includes some repetitions of intermediate steps and results which are also included in section 4.3.

Absorbing the axion field into the quarks
To have a field basis in which all fields are mass eigenstates, we have to diagonalize the mass matrix given by with the two expansion parameters As in the main text, we split the diagonalization into two parts, of which only the first one depends on the axion field. We start by absorbing the axion dependence of the quark mass matrix into the right-handed quarks, This transformation removes the axion field a from M in Eq. (C.1) while the quark kinetic terms generate derivative couplings of the axion to right-handed quarks, where q stands for both up-and down-type fields. The path integral measure is not invariant under this transformation and anomalous interaction terms between the axion and gauge bosons arise. These are [90,91,[123][124][125]] With Eq. (C.4) we have chosen to apply axion-dependent rotations only to right-handed fields which do not couple to SU (2) gauge fields. Therefore, W -couplings are absent. N and E are the anomaly coefficients defined as T (R f ) is the Dynkin index of the SU (3) representation and Q f are the electric charges. In our model, N = 3 and E = 5 or 8 depending on whether the leptons are also charged under PQ. We normalize the axion gluon coupling by defining where N DW ≡ 2N = 6 is the domain wall number counting the number of inequivalent vacua in the QCD-induced axion potential. This means that we can write the gauge-boson coupling terms as Comparing to the EFT setup defined in Eq. (2.1), we can identify c gg = 1, c W W = 0, c γγ = E/N , c ZZ = tan(θ W ) 2 E/N and c γZ = −2 tan(θ W ) E/N .
We have eliminated the a dependence of the mass matrix by the transformation in Eq. (C.3). The subsequent steps in the diagonalization procedure also contain axial transformations of the quark fields, which are however not axion dependent and only result in the usual shift of the theta angle of QCD θ QCD θ QCD → θ QCD + arg(det( M | a=0 )), (C.10) therefore simply displacing the location of the minimum of the axion potential.

Mass diagonalization
We continue with fully diagonalizing the matrices M u and M d by unitary transformations of left-and right handed fields. As before, we drop the labels u/d. At every point, we consider all three generations of SM and F quarks even though the indices are also omitted.
Because M M † is hermitian, it can be diagonalized by a unitary matrix U , where Λ 2 is a diagonal matrix with only real positive eigenvalues. We then define the unitary matrix S = M † U Λ −1 and a unitary transformation of fields Here, q can be either u or d. This transformation diagonalizes the mass matrix We now need to perturbatively find U , which diagonalizes In the last step we have defined the matrices δ, µ and ξ for notational convenience, of which δ and ξ are hermitian. To quadratic order in , U is given by where the unitary matrices U δ and U ξ are defined by the property that they diagonalize hermitian matrices to give the q and F masses. Respectively, (C. 18) In the region of interest for NA62, this expansion in can be done safely, as 10 −4 .
In order to map onto the effective description and to avoid sizable corrections to the effective description in (4.1), we need to ensure that m F i φ such that the F quarks can be integrated out at some scale above the PQ one. In other words, we need the scale separation between the mass of the F quarks and the PQ scale to be sufficiently large. This condition can be written as Here, min eig denotes the smallest eigenvalue of a matrix. In the second and third lines, the size of the eigenvalues are constrained by perturbativity. We see that 1 is only a necessary condition for the F quarks to be much heavier than the PQ scale, while the last line is a sufficient condition. We define one more expansion parameter as (C.22) When 1, we can expand ξ −1 as where we have used a Neumann series from the second to the third line and assumed that eigenvalues of ββ † are at most of order 1. Inserting the leading-order result in into (C.17), we find for the light quark masses In the last step, we have defined the coupling matrices

Axion-quark interactions
Because the axion field was absorbed entirely into the right-handed fields, the derivative terms with left-handed quark fields are not affected by the unitary transformation U . However, the derivative axion coupling to right-handed fields as in Eq. (C.4) does not transform trivially, To leading order in , we obtain for S from which we can determine the relevant coupling matrix to leading order in as where in the last line we have only kept the leading-order term in as in Eq. (C.27), which is of course only justified if 1. So as one would expect, the coupling of SM quarks to axions is only strictly proportional to their masses if the heavy messenger fields are well separated from the PQ scale and can be integrated out. In the final step, we have also used that Λ −1 = diag(M −1 q , M −1 F ). With this, we can finally write the axion couplings as 36) in the form that has been used in the main body of this work.
Tree-level contributions to s → d + a The previous leading-order expansion only induces flavour-diagonal tree-level couplings between the axion and SM quarks. In this case the s → d + a process is only induced at loop level as we have discussed in the main text. However, in general we can not exclude the possibility that higher orders in the expansion of both and , which induce nondiagonal coupling structures, are relevant. These flavour-violating couplings would trigger the s → d+a decay already at tree-level. It therefore becomes important to analyze whether these tree-level couplings can in fact spoil our analysis or, at the very least, whether they restrict the range of values of our expansion parameters.
Before going into the details, let us briefly summarize how large the axial-vector coupling, irrespective of whether it is induced at tree-or loop-level, is allowed to be without being in conflict with the bound BR(K + → π + + a) < 7.3 · 10 −11 [116]. Generally, a coupling of the form results in the decay width of Eq. (2.7) which leads to |h S ds | 1.5 · 10 −12 1 GeV . (C.39) From this we can readily conclude that the expansion is safe: The parameter = v/ φ is O(10 −4 ) for f a ∼ 10 6 GeV, which is the region where NA62 is sensitive. The leading-order axial-vector coupling between quarks and the axion in Eq. (C.36) therefore corresponds to 1/(2 φ ) ∼ 1/(12 · f a ) ∼ 1/(10 7 GeV). The NLO contribution to the axialvector coupling between the axion and SM quarks would be suppressed by an additional factor of 2 with respect to this leading-order coupling if we were to insert higher orders terms in S in Eq. (C.31). This therefore means that the NLO coupling is parametrically suppressed by 10 −7 · 10 −8 GeV −1 = 10 −15 GeV −1 , which is sufficiently far away from the E787 bound quoted in Eq. (C.39) and is also out of reach for NA62.
Next, we want to quantify how small has to be in order for significant flavourviolating axion couplings to be avoided at tree level. By expanding ξ −1 in as in Eq. (C.23), we can write S schematically as and estimate the parametric size of the coupling structure in Eq. (C.29) as Hence, we can put an upper bound on the tree-level flavour violating axion couplings to the SM quarks, To avoid the constraint in (C.39), it is sufficient to require 10 −3 . By further restricting to 10 −4 , the tree-level effect becomes negligible compared to the loop-induced effect. This is exactly the setup which is considered in the main text.
Two more comments are in order. First, note that flavour-violating axion couplings at tree-level are not necessarily a problem but could also be considered an interesting feature of our model. However, since we initially set out to UV complete the model in (4.1) and not to build an axion-flavour model, we choose to restrict ourselves to the case of large scale separations, where the tree-level flavour violation is negligible and flavour violation is induced only by loop processes.
Second, tree-level flavour violation can also be suppressed by the coupling matrices without requiring the masses of the additional fields to be much larger than the PQ scale.
For instance, when = 0.2, we can set β = λ = 1 and α = Y 1+0.2 2 0.2 2 , with Y being the SM Yukawa couplings. This results in vanishing flavour off-diagonal axion couplings to SM quarks in Eq. (C.33) at all orders in , as the inversion of ξ in (C.23) is trivial in this case and thus ABA † becomes proportional to M 2 q , which is exactly diagonal, without expanding in : Moreover, we also see that corrections to the relation in Eq. (C.27), crucial to obtain the effective coefficient Eq. (5.5), are of the order 2 = 2 = 4% such that our discussion of the K + → π + + a decay is still valid. Note however that α in this explicit realization is already very close to the perturbativity limit of Yukawa couplings.

Electroweak interactions
Right-handed q-and F -quarks are in identical representations of the SM gauge group (considering up-and down-type separately). Hence, the interactions of this chiral components with gauge bosons are unchanged under our transformation. This is not the case for left-handed fields. It is useful to write the 6x6 matrix U as a block matrix where each block is a 3x3 matrix. Note, however, that unitarity of U does not imply unitarity of any of the blocks. Because the transformation above mixes different representations of SU(2) × U(1) Y /U(1) EM , we have to check how the interactions with W and Z bosons are modified. Let us start with the Z bosons, In this expression, the upper (lower) sign refers to up-(down-) type quarks and Q is the electromagnetic charge. We see that Z can in principle couple to all available neutral currents, including ones involving light SM quarks of different flavour because A does not have to be unitary. By identifying the blocks A and B in our perturbative result of U in (C.16), we find the Z-interactions at leading order in to be (C.48) Tree-level flavour-changing couplings to SM quarks only appear at order 2 .
In the next step, we look at the W interactions. Because these mix up-and down-type quarks, we reintroduce the corresponding labels to write The W boson couples to all available charged currents. We identify the CKM matrix V as which unlike in the SM does not have to be unitary, but non-unitarity only appears at order 2 . When we again insert the perturbative results for A and B, we arrive at Radial modes of H and Φ So far we have not included the radial modes of the Higgs field H and of the PQ field Φ in our discussion. To capture the impact of the diagonalization procedure on their couplings, it is convenient to write The unitary transformation then results in As the expressions are lengthy, we quote each 3x3 blocks separatel u, c, t + ... = 0 Figure 6. Relevant diagrams for determining the renormalization constants contributing to s → d.
The cross marks the insertion of a counterterm. The dots stand for higher order diagrams in our UV model, such as Z boson-induced flavour changes or diagrams with heavy F quarks.
where we have defined C = U † δ αλ † U ξ 2 A M 2 F / Σ 2 . As we can see, modifications of the Higgs coupling to SM quarks always appear together with additional factors of 2 and are thus negligible. A contribution of φ to s → d + a, on the other hand, must involve internal F -quarks and is therefore proportional to 2 . Consequently, any contributions from φ loops, in particular the ones related to the unsuppressed coupling in Eq. (C.58), become negligible in the 1 limit (as long as this suppression also compensates for possible enhancements originating from the hierarchical coupling matrices). But even if is not small, the φ-mediated flavor violation can be suppressed given an explicit form for the α and β coupling matrices. For instance, in the realisation proposed above Eq. (C.43), the φ-induced loop processes become either flavour-diagonal (because the amplitude has the same flavour structure as the mass matrix) or numerically negligible. This is also true for Higgs-induced loops.

D Counterterm contribution to the kaon decay width
As mentioned in the main text, counterterm contributions are relevant for the computation of the s → d + a transition as a result of the renormalization of quark fields (see the right diagram in figure 3). Most important in this context is the renormalization of SM quark fields, while contributions related to heavy F -fields are of higher order in . Note that in the latter case one also encounters divergent loop integrals, which are parametrically suppressed by O( 2 ). Cancelling all of these divergences (and, more generally, all oneloop divergences in our UV model) requires a more complete renormalization discussion. However, in the following, we focus on the leading order processes that are not suppressed by 2 , and therefore avoid this complication.
We start by parameterizing the renormalized down-type quark fields as [145] The renormalization constants of relevance for the s → d + a transition are determined by demanding that the one-loop W ± contribution to s → d and the counterterm cancel each other [127], as depicted in figure 6. The renormalization constants are therefore O(g 2 ) as a result of canceling the W ± loop. Inserting the renormalized quark fields from Eq. (D.1) into the right-handed derivative axion-fermion interactions, we obtain the terms where the quark fields now denote the renormalized ones. These new terms involving the renormalization constants therefore induce s → d + a at tree-level. Note that this process is of order g 2 / φ , just like the one-loop FCNCs in Eq. (5.5). Hence, there is no reason to neglect this counterterm contribution at this stage. As is known in literature, the explicit renormalization calculation does not need to be performed. For this, as will be proven below, the diagrammatic equation in figure 7 can be shown to hold without knowing the precise expressions for the renormalization constants. Moreover, from the renormalization condition shown in figure 6, we can also conclude the two relations showcased in figure 8 .
Combining all the relations above, we end up with the key relation depicted in figure 9, which allows us to exchange the counterterm calculation with the computation of additional self-energy diagrams where the ALP is emitted from the external down-type quark legs [146], i.e. the diagrams on the right-hand side. From this point it is most straightforward to simply compute the additional loop diagrams instead of taking a detour to compute renormalization constants. Upon performing the calculations, one notices that the diagrams on the right-hand side of figure 9 add up to zero at linear order in the downand strange-quark masses. Therefore, the counterterm contribution is subdominant and can safely be neglected for all practical purposes.
Proof of the relation in figure 7 The derivative axion-quarks interaction can be generally written as T denote the bare down-type quark fields and G Q is a general coupling matrix of diagonal structure, i.e. G Q = diag(g d , g s , g b ). Inserting now the renormalized quark fields from Eq. (D.1), we obtain We also need to trace the influence on the kinetic and mass terms, where M = diag(m d , m s , m b ) and the counterterm Lagrangian L C is given by For the left diagram in figure 7, we use Eq. (D.4) to obtain where Z R * sd is a specific matrix element from δZ R † and * denotes complex conjugation. For the diagram in the middle, we can use Eqs. (D.4) and (D.9) to find where we have used the equation of motionū d (p 2 ) / p 2 =ū d (p 2 )m d . Using the equation of motion once more for / p 2 in the propagator yields which cancels the g s term in Eq. (D.10). Note that the terms involving δZ L cancel due to P R P L = 0.
We can proceed in a similar fashion for the right diagram in figure 7, which cancels the g d term in Eq. (D.10). We therefore conclude that the sum of all three diagrams vanishes, independently of the precise expression of the renormalization constants.

General considerations
In this appendix, we assess whether the one-loop contributions shown in Figs. 10 and 11 can generate quark EDMs and hence also a neutron EDM. Before computing the relevant diagrams explicitly, it is worthwhile to remember the discussion of appendix A. Similar to the case of weak interactions, CP-violating interactions can be parameterized by assigning complex phases to the corresponding operators. Therefore, once a specific operator is inserted twice in a given diagram, i.e. a specific vertex and its hermitian conjugate appear together, the CP phase necessarily drops out.
To be more specific, consider the diagram on the left in figure 10. Let us first write the corresponding a − q − F vertices generally as (∂ µ a)q i γ µ X ij P R F j + h.c. with an arbitrary 3×3 flavour matrix X, and as usual use q (F ) to denote the up-or down-type quark triplets in the SM (F-quark) sector. Note that the right-handed coupling structure is fixed by our diagonalization procedure, see appendix C. This operator violates CP if X has complexvalued entries. Note now that the external quarks are fixed in the Feynman diagram and hence the same a − q − F vertex appears twice in the one-loop diagram because the photon only couples flavour-diagonally to all orders in and . If we rewrite the relevant entry X ql as x ql e iθ ql , where l denotes the internal F -quark flavour and θ ql denotes the CP violating phase, we get e iθ ql · e −iθ ql = 1 in the Feynman amplitude and hence the CP violating phase drops out. As a result, the diagram on the left in figure 10 cannot produce any net CP-violating quark EDM. This is also consistent with a description based on Jarlskog invariants, discussed in ref. [129], upon taking the right-handed coupling structure into account.
An analogous argument can be applied to the remaining diagrams in figure 10, and also to the ones where the photon is radiated off an external quark line in those diagrams. None of these diagrams can therefore induce an EDM at any order in and . This is true even if further CP-violating coupling structures were to occur at higher orders, which could 'allow' for an internal SM quark q = q to appear in the loops (if q = q, no CP violating vertex is possible as discussed in appendix A).
Thus, only the Barr-Zee diagram on the left of figure 11 and its 'conjugate' with the internal gauge boson and axion line interchanged remain. 15 For these diagrams the  Figure 11. Barr-Zee diagrams relevant for the quark EDM computation in our UV model. reasoning is not as straightforward since no vertex appears twice in the diagram. For an internal photon line, however, one notices that no CP violating vertex appears in the diagram. As the photon only couples flavour-diagonally, the internal quark line has to be the same quark q as the external one. In this case, the axion-quark vertex cannot violate CP (see appendix A) and consequently this diagram is fully CP conserving. This reasoning does not apply to the Barr-Zee diagram with an internal Z-boson and hence a contribution to the quark EDM d q given by −(d q /2)F µνq σ µν iγ 5 q can be expected. For the explicit computation below, we follow the momentum flow as shown in the middle and right panels of figure 11.

Z-induced Barr-Zee diagram
We now calculate the EDM contribution based on the Barr-Zee diagram with an internal Z boson in figure 11. Considering the order in and to which we have expanded, to obtain a net CP violating contribution we have to focus on an internal F -quark for the fermionic line. We parameterize the flavour matrix structure of the a − q − F and Z − q − F vertices by X ql = v(M −1 q ABM −1 F ) ql /2 and Y ql = g 2c W A ql , respectively, as dictated by Eqs. (4.22) and (5.2). For both diagrams and for each individual F -quark with index l, the amplitude reads iM (±) d 4 k (2π) 4 ū q (p 2 ) iX ql i(/ k + / p 1 − / p 2 ) P R i(/ k + / p 1 + m F l ) Here, we have schematically used the Feynman rule due to (e/4π) 2 (a/f a )F µνZ µν resulting from the chiral rotation discussed in section 4.3. The global ± depends on whether q is an up-or down-type quark, see Eq. (5.2), and * α is the photon polarization vector. We can write the amplitude in a more compact form as

(E.2)
We then use Package-X [126] to compute the loop integrals and map out the relevant Lorentz structure for the quark EDM. Both of the Barr-Zee diagrams give the same loop functions f (m F l , m a , m Z , m q ) ≡ f (m F l , m q ) for the EDM operator, differing only in their sign. The relevant amplitude structure therefore reads iσ αβ γ 5 2m q * α (p 1 − p 2 ) β u q (p 1 ).

(E.4)
In the case considered in the main text, where all F quarks are assumed to be equally heavy, performing the sum over all internal F -quarks, (i.e. over the index l) leads to a vanishing amplitude. Higher orders in the expansion of the relevant Z − q − F and a − q − F vertices can in principle again lead to a non-zero result, but one that would be suppressed by higher powers of and would therefore be negligible.
As the Z-induced Barr-Zee diagram is formally of two-loop order and to motivate that the estimate for the upper limit of the neutron EDM in Eq. (5.6) is reasonable, we can, for the sake of an estimate, assume that all F quarks have different mass and hence no cancellation takes place. The loop function f (m F l , m q ) is UV divergent and hence exhibits a scale dependence µ. This dependence would of course disappear if further contributions were taken into account. Expanding in inverse powers of m F l , we find Without performing a full two-loop analysis, we conservatively assume the term in the brackets to be an O(1) factor. The remaining finite terms in f (m F l ), which are not expected to fully cancel, would actually be significantly smaller for m F ≥ 1 TeV. Moreover, if we parameterize the entries of (AB) ql and A ql by x ql e iθ ql and a ql e iφ ql , respectively, we obtain x ql e iθ ql a ql e −iφ ql − a ql e iφ ql x ql e −iθ ql iσ αβ γ 5 2m q * α (p 1 − p 2 ) β u q (p 1 ) . (E.6) The EDM operator to map onto is −(d q /2) F µνq σ µν iγ 5 q, where d q has to be real-valued for the operator to be hermitian. By comparing with Eq. (E.6) and including the contributions from all three F quarks, we have x ql a ql 2 sin(θ ql − φ ql ) (E.7) where we have identified x ql a ql ∼ 2 m 2 q /( 2 4 Σ 2 ) in the second line based on Eq. (C.27). For the up-and down-quarks respectively, fixing f a = 4 · 10 6 GeV and assuming sin(θ ql − φ ql ) ≈ 1, we get  [130]. Of course, it should be kept in mind that the result in Eq. (E.10) is only an estimate as well. However, given that the sensitivity of current neutron EDM experiments is away by several orders of magnitude, it is an estimate that is sufficient for our purposes.