Charming ALPs

Axion-like particles (ALPs) are ubiquitous in models of new physics explaining some of the most pressing puzzles of the Standard Model. However, until relatively recently, little attention has been paid to its interplay with flavour. In this work, we study in detail the phenomenology of ALPs that exclusively interact with up-type quarks at the tree-level, which arise in some well-motivated ultra-violet completions such as QCD-like dark sectors or Froggatt-Nielsen type models of flavour. Our study is performed in the low-energy effective theory to highlight the key features of these scenarios in a model independent way. We derive all the existing constraints on these models and demonstrate how upcoming experiments at fixed-target facilities and the LHC can probe a vast region of the parameter space, which is currently not excluded by cosmological and astrophysical bounds. We also emphasize how a future measurement of the currently unavailable meson decay $D \to \pi + \rm{invisible}$ could complement these upcoming searches and help to probe a large unexplored region of their parameter space.


I. INTRODUCTION
One of the outstanding open questions in particle physics is the nature of dark matter (DM) and whether it is part of a larger dark sector that we yet have to discover. Most realistic models require some form of nongravitational interaction between us and the dark sector in order to satisfy cosmological constraints. These "portals" then offer the possibility to probe the dark sector through laboratory experiments or astrophysical observations.
An intriguing possibility is that the portal to the dark sector is flavour sensitive or even connected to an ultraviolet (UV) theory of flavour. Simple flavoured dark matter scenarios have received a lot of attention recently [1][2][3][4], since they allow probes of dark sector physics using low energy flavour observables. For strongly coupled dark sectors, a flavoured portal imprints the Standard Model (SM) flavour structure on the dark sector, leading to a variety of new phenomena [5,6].
So far only couplings to down-type quarks were considered in this context, therefore it is natural to ask what new aspects arise if the portal couples to up-type quarks instead. 1 The key feature in this case is the emergence of a light pseudo Nambu-Goldstone boson (pNGB) which dominantly couples to up-type quarks, and which we therefore dub a charming ALPs.
Our main goal in this work is to develop the effective theory of the charming ALP and its phenomenologi-cal profile, independently of its embedding in different UV scenarios. Besides QCD-like dark sectors [5,[10][11][12][13], these particles can arise e.g. in specific Froggatt-Nielsen (FN) models [14] where only right-handed (RH) up-quarks have non-zero charges. The different UV completions provide some guidance for the structure of the effective couplings, which we use to define benchmark scenarios for the charming ALP.
We will study the phenomenology of these different benchmark models through their low-energy effective field theory (EFT), taking into account the effects of QCD confinement for low enough ALP masses. A lot of work has been done in this arena, see e.g [15][16][17][18][19][20][21] for the study of ALP collider signatures, [22][23][24] and [25][26][27] for the study of flavour-changing neutral currents (FC-NCs) in the quark and lepton sector, respectively, as well as [28][29][30][31] for the calculation of the one-loop running. It is worth to mention also [32,33] regarding CP-violating probes of ALPs.
The work is organized as follows: We introduce the effective Lagrangian describing the charming ALPs and motivate briefly the four particular benchmarks models studied in this work in section II. In section III we examine the different flavour constraints relevant for charming ALPs. More specifically, we consider D −D mixing, exotic decays of D, B and K mesons as well as the decay J/ψ → aγ. We discuss the bounds arising from astrophysical observables and cosmology in section IV, describing also the different ALP decay channels. In section V we study the different collider probes on the models at hand, including upcoming fix-target experiments as well as those at LHC forward detectors. We combine all these different bounds and discuss the resulting constraints on the parameter space of the models in section VI. We conclude in section VII. Finally, we present in some detail the particular UV completions considered arXiv:2101.07803v2 [hep-ph] 29 Apr 2021 in this work in appendices A and B.

II. CHARMING ALP EFT
We consider a general ALP, which we will denote a, with flavour-violating couplings to RH up-quarks. The most general EFT describing such a system is given by the following Lagrangian [34,35] where g 1 , g 2 and g 3 are the gauge couplings of U (1) Y , SU (2) L and SU (3), respectively, whereas B µν , W I µν , I = 1, 2, 3, and G a µν , a = 1, . . . , 8, are their corresponding field-strength tensors.
Furthermore,B µν = 1 2 ε µναβ B αβ , . . . , denote their corresponding duals, while H stands for the SM Higgs doublet. The Wilson coefficients (WCs) c g , c W , c B and c H ∈ R, whereas c u R is a hermitian matrix. In order to write down the above Lagrangian we have assumed that a is the pNGB of the spontaneous breaking of some global U (1) symmetry, which is softly broken and may be anomalous. We have also assumed that the couplings to leptons, SM quark doublets and RH down-type quarks vanish. Here, contrary to the QCD axion case, we will treat m a and f a as independent parameters. Using field redefinitions, we can trade the operator O H = (∂ µ a/f a )H † i ← → D µ H by the flavour-blind and chirality conserving one (see e.g. [16,34]) Together with this operator induces flavour-changing neutral currents (FCNCs) at one-loop, which have been studied in [24] in the framework of B and K-meson decays. Here we will just assume that both WCs are small enough so that the leading flavour-violating effects are parametrized by c u R . Furthermore, after integrating by parts and using equations of motion, one can express this chirality-conserving operator as a function of plus some extra contributions to the anomalous terms in (1), whereH = iσ 2 H * and Y u is the up Yukawa matrix, Note that, without any loss of generality, we can always choose a basis where with λ u,d diagonal matrices with real and positive entries andṼ a unitary matrix. In the case where there is no extra contribution to the fermion masses,Ṽ is just the CKM mixing matrix V and λ u = and v = 246 GeV the Higgs vacuum expectation value. In this basis, the RH up-quarks do not need to be rotated to diagonalize the mass matrices generated after electroweak symmetry breaking. Indeed, one could just take Henceforth, we will assume that the above EFT Lagrangian is defined in such a basis. Then, if we denote the WC of the operators in (4) by C, one has that For small ALP masses, m a 1 GeV, a will mostly decay to hadrons. These decays will proceed through the following Lagrangian [34][35][36][37] where B 0 is a constant, f π ≈ 93 MeV is the pion decay constant and m q is the quark mass matrix m q = diag(m u , m d , m s ). In the above equation, where is the Goldstone matrix describing the spontaneous symmetry breaking SU with e = g 2 g 1 / g 2 1 + g 2 2 the electric charge and Q q = 1/3 diag(2, −1, −1). One should note that, in order to get the above Lagrangian, we had to get rid of the gluon coupling by the following chiral transformation Note that this Lagrangian gives an irreducible contribution to the ALP mass where Kinetic mixing arising from the last term in eq. (9) induces a mass mixing between the different neutral pions. In particular, we obtain where and As already mentioned, we are interested in scenarios where the ALP only interacts with RH up-quarks, and it is likely to mediate flavour-changing processes. In order to explore how the different experimental constraints are intertwined and the best way to probe these models, we will consider four different benchmarks. In practice, each of these scenarios corresponds to a particular choice of the matrix c u R defined above.
The first two benchmarks are motivated by theories of 'dark QCD', where the SM is extended with a confining dark sector, composed of n d dark flavours transforming under SU (N d ). These models constitute a particular UV completion of the scenarios we have in mind, where light pNGB bosons do only interact at leading order with the SM RH up-quarks. Indeed, this is a natural outcome when both sectors are mediated by a scalar, bifundamental of both confining groups, with hypercharge −2/3. We refer the reader to appendix A and references therein for more details. At the end of the day, after integrating out the heavy scalar mediator, and assuming confinement in the dark QCD group, we obtain almost degenerate, parametrically light scalars with a Lagrangian along the lines of (1). When studying phenomena where the interplay of the different scalars is not relevant, one can examine the phenomenological impact of the different degrees of freedom separately. In particular, focusing on the 'diagonal' dark pions π D3 and π D8 , we obtain 5c 12 s 12 0 5c 12 s 12 4c 2 12 + 9s 2 12 as well as c H = 0 and c g = c W = c B = 0 at tree level. In the equation above κ 0 ∈ R + and c 12 = cos θ 12 , s 12 = sin θ 12 , with θ 12 ∈ [0, π]. For the sake of concreteness, following [5], we fix θ 12 = 0.022. There is another class of models that can naturally UV complete these scenarios. It is the case of FN models where only RH up-quarks have non-zero charges. Such models are attractive because they naturally result in enhanced Yukawa couplings, while still being in agreement with existing flavour bounds. They where already considered in a slightly different context [38][39][40] but without paying attention to the phenomenology of the light scalar degree of freedom, the flavon. We refer the reader to appendix B and references therein for more details. One can reproduce the hierarchical up-quark masses with RH up-quark charges n u = (2, 1, 0) under the U (1) flavour group, assuming that the vev of the scalar breaking such symmetry is ∼ m c /m t times its UV cutoff. At the end of the day, this setup leads to an ALP Lagrangian along the lines of equation (1) with whereas c H = 0. The concrete values of the anomalous couplings c g , c W , c B depend on the specific UV completion of the FN model and may all be zero, which is the case we will consider in the following. We define the FNmotivated benchmark with c u R given as above, whereas the other WCs are zero. Finally, we also consider an infra-red (IR) motivated scenario, where (c u R ) = 1, ∀i, j at the scale f a . This is representative of the anarchic limit, where no hierarchies are present in c u R and all the entries are of the same order. Similarly to the previous cases, we assume that the anomalous gauge couplings and c H are negligible.

III. FLAVOUR CONSTRAINTS
The presence of flavour-changing ALP couplings to RH up-quarks will induce several flavour violating processes, constraining significantly the parameter space. In particular, we will consider As shown in fig. 2, in the models at hand, exotic D meson decays are tree level processes while B and K decays can only happen at one loop. In addition, ALPs also contribute to radiative J/ψ decays (c.f. fig. 3).c The effective Hamiltonian relevant for D meson mixing reads [41] where andÕ 1,2,3 are obtained from O i after exchanging both chiralities, i.e., L ↔ R. Henceforth, we will just be focusing on the new physics contribution to these WCs, i.e., Depending on the specific mass of the ALP, such contributions will involve either short or long-distance physics. The first case occurs when integrating out a heavy enough ALP, m a m c , whereas the second one is the consequence of applying naively the operator product expansion (OPE) in powers of ∼ 1/m c to the D −D system, when m a m c . In the first case, one obtains and zero elsewhere, while in the second onẽ with all other WCs vanishing. In general, where D 0 |Õ i |D 0 = D 0 |O i |D 0 , for i = 1, 2, 3, due to parity conservation of QCD and m D = 1.865 GeV [42]. For small ALP masses, we get [43] |M where O 2 = −0.1561 GeV 4 at µ = 3 GeV, see [43], and we have neglected O(m u /m c ) corrections. However, when studying short-distance physics, it is also necessary to run the different WCs from the scale of integra- Due to the running, to evaluate M 12 we also need O 3 at µ = 3 GeV, which reads 0.0464 GeV 4 [43]. As an example, for Λ = m a = 2 TeV and µ = 3 GeV, we obtain We demand that the new physics contribution to x 12 = 2|M 12 |/Γ does not exceed its upper bound at 95% confidence level (CL) [45], i.e., where Γ D = 1.60497 · 10 −12 GeV [42].

B. Exotic D, B and K decays
We study ∆F = 1 decays of the form M → N a with M = D ±,0 , B ±,0 , K ±,0 and N = π ±,0 , K ±,0 . The corresponding parton level Feynman diagram are shown in fig. 2. The associated matrix element can be decomposed as with k µ = (p − p ) µ the momentum transfer and q i and q j the relevant quarks for the decay at the parton level. The scalar form factor is then defined as and the resulting decay width is given by where κ M N is defined by In the model at hand and neglecting small isospinbreaking effects, the exotic D meson decay D ±,0 → π ±,0 a is induced by the dimension-5 operator The width for such decay channel can be read from equation (39) by simply replacing κ M N with (c u R ) 12 , m M = m D , m N = m π and using f Dπ 0 (m 2 a ) from [46]. On the other hand, the one-loop running of c u R from the UV scale f a to the IR scale µ will generate a term 2 at low energies. Indeed, one obtains [28][29][30][31] After EWSB, this operator leads to where This WC is responsible for the exotic decays B → Ka, B → πa and K → πa, where B = B ±,0 , K = K ±,0 and π = π ±,0 . 3 The corresponding expressions can be read from equation (39) after taking µ ∼ m t and There are no constraints on the branching ratio Br(D → π + invisible) to date. However, there are measurements of the three-body meson decay D + → (τ + → π + ν)ν [50,51]. Since these analyses show the event distribution as a function of the missing mass squared M 2 miss (which would correspond to m 2 a in the ALP case), one could recast them to constrain the branching ratio D + → π + a. 4 This was done e.g. in [29] for the massless axion by concentrating on the bins with M 2 miss ≤ 0.05 GeV 2 . Since, as we will see, the total width of the ALP in all the benchmarks under consideration is very small, one can safely produce a similar bound for different values of m a by comparing the observed number of events with the predicted background for every bin having M 2 miss ≥ 0. More precisely, we derive 90% CL on BR(D → πa) by using the TLimit class of ROOT [53], which implements the CL s method [54] and allows to include systematic errors in the background and signal. The bounds arising from [50] turn out to be stronger than those resulting from the use of the more recent experimental analysis in [51].
Regarding exotic meson decays involving down quarks, there are several analysis focused on a massless axion X 0 , like K + → π + X 0 [55] and B ± → π ± X 0 , B ± → K ± X 0 [56]. There are no searches for a massive ALP a, with the exception of the recent analysis in [57], where bounds on K + → π + a as a function of m a were presented. Similarly to the D ± → π ± a case, we fill this gap by recasting existing searches on three-body decays, where the relevant kinematic information is provided. In particular, we derive constraints on B → Ka and B → πa by recasting the searches performed in [58] and [59], respectively. More specifically, we set 90% CL on B → Ka by combining the observed number of events and the predicted background for B + → K + νν and B 0 → K 0 νν for every figure 5 of [58] with the CL s method. Finally, for the case of B → πa, we derive 90% CL limits on B → πa with the CL s method by comparing the observed number of events with the predicted background for every p π ≡ p 2 π bin in the right panel of figure 4 in [59] (which is in one-to-one correspondence with the ALP mass via . On the other hand, it is expected that Belle II will be sensitive to the SM Br(B ± → K ± νν) = (4.0±0.5)×10 −6 [60] at 10% accuracy with 50 ab −1 of data [61], whereas NA62 will measure the branching ratio Br(K ± → π ± νν) to within 10% of its SM value Br(K ± → π ± νν) = (8.4 ± 4.1) × 10 −11 [62,63]. Regardless of whether these collaborations publish limits directly on a two-body decay or a recast of the three-body decay analysis is needed, such numbers represent a great improvement with respect to current bounds. Diagonal ALP couplings to charm quarks are strongly constrained by charmonium decays like J/ψ → aγ, as first proposed by [65] and later studied by many others, see e.g. [66][67][68][69][70][71]. The parton-level diagram prompting such decay is shown in figure 3. In order to absorb some of the QCD uncertainties of the calculation, it is convenient to normalize the corresponding branching ratio by the one of J/ψ → µ + µ − , which is accurately measured Br(J/ψ → µ − µ + ) = 5.973 % [72]. One then obtains where F ∼ O(1/2) is a correction factor accounting for QCD effects [73,74], contributions related to bound-state formation [75,76] as well as relativistic corrections [77]. For the sake of concreteness, we will assume that F = 1/2 henceforth. This leads to This decay has been searched for by the CLEO collaboration [78], which we will use to constrain the benchmark models at hand.

IV. ASTROPHYSICAL AND COSMOLOGICAL BOUNDS
A. Bounds from supernova SN1987a The observed neutrino burst due to the core-collapse supernova SN1987a can impose constraints on the ALP parameter space. Since neutrino emission constitutes the main cooling mechanism for the proto-neutron star resulting from the collapse, a too large ALP emission could compete with this cooling mechanism and eventually conflict the observed amount of neutrinos. Following [79], we will impose that the ALP luminosity in the protoneutron star L a does not exceed the neutrino one L ν , i.e., L a ≤ L ν = 3 · 10 52 erg/s. 6 The ALP luminosity in the proto-neutron star is given by [81,82] where ω is the ALP energy and R ν ∼ O(40 km) is the radius of the neutrinosphere, beyond which neutrinos free stream until arriving to the Earth, and we have taken into account the probability e −τ for an ALP produced within the neutrinosphere to reach R far ∼ O(100 − 1000 km), after which neutrinos are not produced efficiently. If this is not the case, ALPs being produced within the neutrinosphere get 'trapped' due to their large couplings and their energy is eventually converted back into neutrinos. Such probability is computed with the help of the optical depth τ = τ (m a , ω, r, R far ), for which we will take R far = 100 km [81]. In the above expression, dP a /dV dω is the ALP differential power. For the charming ALPs considered here, the main channel will be the bremsstrahlung process N + N → N + N + a, since the sole tree-level couplings are those to up-type quarks and therefore to nucleons. Such differential power is given by [79,83] dP a dV dω where T is the temperature as a function of the radius, β is a phase space factor β = 1 − m 2 a /ω 2 and Γ a is the ALP absorption width. The latter is given by with [82] In the above equation, c aN N is the ALP-nucleon coupling, which reads (see appendix C for more details) while Y N ( ) is the mass fraction of the nucleon N ( ) , n B is the baryon density, n B = ρ/m N , and σ npπ is given by with α π ≈ 15. For concreteness we take Y p = 0.3 and Y n = 1 − Y p = 0.7. Moreover [84], while we use γ p = s(n B , Y N , ω/T, m π /T ) with s given by eq. (49) of [85]. 7 On the other hand, following [64], we assume that [86] γ h = −0.0726502 ln(ρ) + 10 10 /ρ 0.9395710 + 2.5558616, where the density ρ is expressed in g cm −3 . Similarly to [64], we assume for ρ(r) and T (r) and the "fiducial" profiles of [81] ρ with k ρ = 0.2, k T = −0.5, ν = 5, R c = 10 km, T c = 30 MeV and ρ c = 3 · 10 14 g/cm 3 . We define R ν as the distance at which the temperature is 3 MeV, obtaining R ν = 39.81 km. Finally, for the optical depth we take [64] Rν r dr Γ a (r). (59)

B. Bounds from red giant burst
The one-loop running of the dimension-5 effective Lagrangian will generate ALP couplings to electrons at low energy. Such couplings face very strong astrophysical bounds for small values of m a . This effect can be particularly relevant when the ALP couples to the top quark, since it will contribute significantly to the running [28][29][30][31]: which leads after field redefinition to After EWSB, these operators lead to

C. ALP lifetime, branching ratios and cosmological bounds
Cosmological bounds are very sensitive to the total decay width and the different branching ratios of the ALP, that we discuss in the following. We do this both for small values of the ALP mass, where QCD is confined and one can use chiral perturbation theory, as well as for larger values where the dominant ALP decays can be computed using quark-hadron duality [90,91]. Following [5], we determine the energy scale separating both pictures by demanding that the total decay width to hadrons or SM quarks is of the same order in both regimes, which leads to ∼ 1 GeV for the benchmark models at hand.
For small masses, m a 1 GeV, the ALP decays to two photons via the mixing (18) and the subsequent decay π → γγ, plus one-loop contributions coming from the integration of heavy quarks. This leads to [18,21] where This decay channel will be the only one present, whenever the one-loop lepton decays a → + − are kinematically closed. The leptonic decay widths can be written as The diphoton final state will dominate over the e + e − and µ + µ − decays in models where (c u R ) 33 is absent or negligible. Otherwise, once a → e + e − and a → µ + µ − open up, they will become the main ALP decay channel, at least for large values of f a leading to log-enhanced g a couplings. At any rate, the m 3 a dependence of the diphoton decay will make Γ(a → γγ) increase faster than the dilepton decay width with increasing values of m a . In some cases, this can turn such decay channel into the leading one, once more, for larger ALP masses before a → 3π opens kinematically. Such decay channel will always dominate the a → γγ final state in our models, with its decay width reading [18] Γ(a → π a π b π 0 ) = π 12 with and λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc). However, when (c u R ) 33 is sizeable and/or f a high enough to have significant g a couplings, a → µ + µ − can be the dominant decay channel in this region of ALP masses.
For masses m a 1 GeV, equation (9) becomes invalid and the dominant ALP decays into hadrons can be computed using quark-hadron duality. The ALP will then decay into gluons and quarks. While the first decay is loop induced, the tree-level decay into quarks will be suppressed by the light Yukawa couplings. Close to threshold, a →q ( ) q will be the leading decay channel, whereas for larger values of m a a → gg will take over. In this regime, the a → γγ decay will always be sub-leading and read We show in figure 4 the different branching ratios for the models at hand for f a = 10 4 TeV, assuming that κ 0 = 1 in the dark-QCD motivated benchmarks. We can see that in the cases where the ALP coupling to top quarks is absent or suppressed (top-left and bottom-right plots), a → γγ is the leading decay mode for most of the region m a 1 GeV until a → 3π is kinematically allowed. In the cases where this coupling is present (topright and bottom-left plots), and for this choice of f a , a → e + e − becomes the leading channel until a → µ + µ − opens up. For ALP masses 1 GeV, decays into hadrons are by far the dominant channels, with a → gg leading at large masses. Now that we have computed the different branching ratios and lifetimes, we can evaluate the impact of the cosmological constraints. It is important to stress that most of them are derived assuming that a only interacts with photons. However, as discussed in [92], one can apply these cosmological bounds to the more general case where other couplings are present. At the end of the day, we will be able to recast the limits from [92][93][94] by using the ALP lifetime 1/Γ, with Γ the ALP total decay width. These bounds include the possible impact on N eff , potential distortions of the cosmic microwave-background spectrum as well as modifications of the predicted bigbang nucleosynthesis, see [92][93][94].
Indeed, the bounds can directly be applied when the decay to lepton pairs is dominant, i.e. when there is a sizeable coupling of the ALP with the top quark, provided one interprets 1/Γ as the total lifetime. 8 In the cases where a → 3π dominates, bounds from 4 He overproduction -the dominant constraint in this regionwill still hold regardless of the changes in the branching ratios, since only a minimal amount of charged pions is enough for this bound to apply. For even larger masses, ALP decays into hadrons will eventually make its lifetime shorter than a second, making nucleosynthesis constraints harmless. Therefore, even for these masses, we can apply the corresponding bounds if we interpret τ as the total lifetime.

A. Fixed target experiments
The main production mode for charming ALPs at fixed-target experiments is the decay of D mesons. We consider NA62 [95] and the proposed SHiP experiment [96] as possible detection experiments. We also consider the bounds imposed by the CHARM experiment in [97]. The geometrical outlines of the experiments are listed in table I. NA62 operating in beam dump mode, meaning the target is lifted so that the 400 GeV proton beam hits the Cu collimator located 20 m downstream, can be used to search for hidden sector particles [98]. A short run in beam dump mode in November 2016 provided useful information about the backgrounds. It was found that an upstream veto in front of the decay volume could reduce the background to nearly zero [99]. The layout of the SHiP detector is proposed with the aim of reducing the beam-induced backgrounds to 0.1 events [96,100], so that 3 decay events correspond to the expected exclusion region at over 95% CL.
The total number of dark pions decaying inside the respective decay volume is with ε geom the geometric acceptance, defined as the fraction of ALPs with lab frame momentum at the acceptance angle of the respective detector and F decay the fraction of ALPs that decay inside the decay volume of the respective detector. F decay and ε geom are calculated following [5]. The D meson momentum distribution for SHiP is taken from [101]. The same distribution is used for the NA62 case as the proton beam is the same. The CHARM experiment searched for ALPs decaying into pairs of photons, electrons and muons. In [97] no events where found, so that we set a bound at 90% CL at N obs = 2.3 events [102]. We assume that for our model the main production channel for ALPs is the decay D ± → π ± a. CHARM also has a 400 GeV proton beam, so we again use the same momentum distribution for the D mesons. The number of protons on target is 2.4 × 10 18 [97] and such a proton beam has a probability of 1.7 × 10 −3 to produce a pair of c quarks [103], leading to 4.08 × 10 15 produced D mesons. Multiplying eq. (71) by i=γ,e,µ Br(a → ii) to take into account that in [97] only the γγ, ee and µµ final states were searched for we use the same procedure as for NA62 and SHiP to impose the 90% CL bound from CHARM.

B. LHC forward detectors
FASER [104] and the proposed MATHUSLA detector [105,106] are designed to detect long-lived particles produced in proton-proton collisions at LHC. Hadron collisions have the advantage of additional production modes, FIG. 4. Branching ratios of the ALP as a function of its mass, ma, for the different benchmark models. The top panels correspond to the dark-QCD inspired case with a = πD 3 (top left) and a = πD 8 (top right), whereas the bottom panels show the anarchic scenario (bottom left) and the FN motivated benchmark (bottom right), respectively. In all cases, we have assumed fa = 10 4 TeV. Moreover, for the dark-QCD motivated scenarios illustrated in the top panels we have taken κ0 = 1. We illustrate with a gray narrow band around 1 GeV, the matching between the calculations performed using chiral perturbation theory and with quark-hadron duality.
such as production via gluon fusion. The detector parameters are given in table I. We focus again on the production via D meson decays. The meson momentum distribution was simulated using FONLL with CTEQ6.6 [107]. As for SHiP and NA62 the number of dark pions decaying inside the decay volume can be calculated using eq. (71). For LHC run-3 N D = 1.1 × 10 15 , which increases by a factor 20 at the High Luminosity (HL)-LHC. FASER will operate at LHC run-3, while FASER2 and MATHUSLA are under consideration for the HL-LHC. Following the procedure as described for SHiP and NA62 the number of dark pions decaying inside FASERs, FASER2s and MATHUSLAs decay volume is calculated. At least three events must decay inside the respective decay volume for a discovery.

VI. RESULTS
The resulting constraints on the parameter space of the four benchmarks models under consideration are displayed in figures 5 and 6. More specifically, we show in figure 5 the different bounds for the dark-QCD motivated benchmarks, with left panels corresponding to the case a = π D3 and the right ones to a = π D8 , respectively. For these two benchmarks we further assume κ 0 = 1 in order to evaluate the logarithms coming from the one-loop running, which only depend on f a . Still, different choices  of κ 0 ∼ O(1) will not significantly affect the resulting bounds shown in the figure. On the other hand, we show in figure 6 the corresponding bounds for the anarchic (left panels) and FN (right panels) scenarios.
In both figures, we represent in dark yellow the region of the parameter space excluded by D meson mixing, whereas the bounds from exotic meson decays B → Ka, B → πa and K → πa are displayed in green, blue and red, respectively. Moreover, for the sake of illustration, one can roughly interpret the expected bounds on B → Kνν from Belle II as sensitivity to B → Ka, which we show in light green. Similarly, one could do the same thing in the case of K → πνν from NA62, which is represented in light red. A proper recast of these two upcoming experimental results will most likely result in slightly weaker bounds. On the other hand, we show the limits arising from the recast of D + → (τ + → π + ν)ν as a red line. Finally, we show in purple the bounds resulting from J/ψ → aγ searches by the CLEO collaboration. The cosmological bounds discussed previously are shown in gray in both figures, with the constraints arising from red giant bursts and from SN1987a exhibited in lilac and light blue, respectively.
The region where the ALP mass lies above the D and B meson masses is only weakly constrained by flavour observables, and is open for probes at the energy frontier, i.e. the LHC and future colliders. Below the D meson mass, some viable regions remain which can be probed by the upcoming or proposed collider experiments discussed in section V. The projected bounds from SHiP and NA62 are represented by dark red and orange dashed contours, respectively. Above these lines more than three events are expected. Furthermore, the discovery lines for FASER with LHC run-3 and FASER2 at HL-LHC are displayed in dashed light and dark blue, respectively, whereas the detection line for MATHUSLA is pictured as a turquoise dashed contour line.
In order to appreciate better the region which can be probed by these upcoming experiments, we show in the lower panels of figures 5 and 6 a zoomed-in view of the area which may be reached. While FASER will mainly validate the constraints from the CHARM experiment shown in light yellow, FASER2 at the HL-LHC as well as NA62 will probe new regions of parameter space, with SHiP and MATHUSLA covering the remaining unexplored areas below the charm mass threshold.
A possible measurement of Br(D → π+invisible) could provide a complementary test of these parts of the parameter space, and might be crucial to probe the region close to the charm mass at relatively large coupling. In the FN inspired model as well as the π 3 D scenario this region of parameter space is not otherwise accessible, while in the π 8 D scenario it will be probed by Belle II, and in the anarchic case it is already excluded. Similarly, the low mass region of the π 3 D scenario is only accessible via Br(D → π + invisible) and future NA62 measurements. To display the discovery potential of such a measurement, we show black contour lines corresponding to values of Br(D ± → π ± invisible) ∈ {10 −8 , 10 −7 , 10 −6 , 10 −5 , 10 −4 , 10 −3 , 10 −2 , 10 −1 }. This demonstrates that providing an experimental measurement of the exotic meson decay D → πa is paramount to leave no stone unturned in the quest for well motivated extensions of the SM that feature charming ALPs.

VII. CONCLUSIONS
In this work we have presented several examples of models featuring charming ALPs, i.e., light pNGBs having off-diagonal couplings with SM up-quarks, and studied in detail the phenomenology associated with their low-energy EFTs. More specifically, we have studied the constraints arising from flavour experiments, astrophysics and cosmology as well as planned fixed-target and collider experiments in four benchmark models. We have shown that such scenarios have still a large unexplored parameter space. We have also demonstrated how future collider and fixed-target experiments can probe these models and that they could be perfectly complemented by the measurement of the exotic decay D → π + invisible, which is currently unavailable. We thus encourage our experimental colleagues to proceed with such measurement. In the absence of dedicated searches, we have also derived bounds on the parameter space of the models by recasting three-body meson decays like D + → (τ + → π + ν)ν or B → K/π νν.
The scenarios considered here can be the low-energy EFT of several compelling UV completions. We have presented two of them: the case of a QCD-like dark sector interacting with the SM via a heavy scalar mediator with hypercharge −2/3, and a FN model of flavour where only RH up-quarks and a heavy scalar have nonzero charges with respect to an spontaneously broken global U (1) symmetry. Some of the phenomenology studied here may change when considering the whole picture in the dark-QCD case, since there might be a non-trivial interplay between the complete set of pNGBs in some regions of the parameter space. In the present work, we have focused on the phenomenological aspects expected to hold when singling out one of such light states. The complete dark-QCD model will be studied in an upcoming work, where a detailed study of the full charming dark sector including the possible connection with dark matter will be presented. A particularly interesting aspect of such scenarios to be studied is the collider phe- nomenology involving rare top decays t → ca as well as the phenomenology of 'charming' emerging jets.
A final region of parameter space that was left unexplored in our study is that of very small ALP masses and couplings, i.e. the lower left regions of our figures. There the charming ALP would be a stable dark matter (DM) candidate. The freeze-out of such a DM candidate is excluded for m a 100 eV due to DM overproduction [92] and in the whole parameter region if structure formation bounds are also taken into account [108]. On the other hand, if the charming ALP DM is produced via a "freeze-in" mechanism, a region of parameter space remains viable [108] in principle. A more detailed exploration of this region of the charming ALP is left for future work.

Appendix A: A dark QCD UV completion
One particular class of theories which can provide a UV completion to the charming ALP EFT is what can be collectively denoted by 'dark QCD', see e.g. [5,11,12]. In these scenarios, one assumes that the SM is extended with a new dark QCD-like gauge group SU (N d ) d with n d Dirac fermions, Q α , α = 1, . . . , n d , singlets of the SM gauge group and transforming in the fundamental representation of SU (N d ) d . For concreteness one can assume N d = 3 = n d , which allows in particular the QCD-like sector to confine at a scale Λ dQCD . Both sectors talk to each other through the coupling to a heavy scalar mediator, X , transforming as a (3,3) under SU (3) ⊗ SU (3) d , which is also charged under SU (2) L ⊗ U (1) Y as 1 −2/3 . Such scalar mediator is naturally heavy, also in agreement with LHC bounds. A similar setup was presented first in [11] but with a different assignment of hypercharge, Y = 1/3, such that only couplings to RH downlike quarks were allowed. It was shown in [5,12] that this scenario lead to emerging jets and its flavour phenomenology was studied in [5]. Here, we focus on a different hypercharge assignment, which leads to a distinct phenomenology. The structure of the model is sketched in figure 7.
The Lagrangian of the dark sector reads where we use greek (latin) indices to denote the flavour indices in the dark (visible) QCD sector and we do not need to explicitly write the corresponding covariant derivatives. Similarly, we did not make explicit color or dark color indices with the exception of G a µν , a = 1, . . . , 8, the field-strength tensor for the dark QCD. In general, the coupling matrix κ αi can be expressed as where both U and V are 3×3 unitary matrices and D is a 3×3 diagonal matrix. In the case where m Qαβ = m Q δ αβ , there is a U (3) d flavour symmetry in the dark sector, which can be used to rotate away V . We will assume that this is the case henceforth. The matrix U can be further expressed as the following product where U ij are unitary rotations between the flavours i and j. For example, U 12 reads where c 12 = cos θ 12 , s 12 = sin θ 12 . Moreover, following [2,5] one can express D as where κ 0 ≥ |κ 1 + κ 2 |. For the sake of concreteness, we will consider the benchmark model defined by We also assume that all this happens in a basis where the SM up Yukawa matrix is diagonal. If the mass of scalar mediator is much heavier than Λ dQCD and Λ QCD , one gets in addition to the SM Lagrangian after integrating out X and using Fierz identities. Similarly to QCD, in the limit m Q → 0 and m X → ∞, the dark sector features a global dark chiral symmetry SU (3) dL ⊗ SU (3) dR , which we assume is spontaneously broken to its diagonal group SU (3) dV by the dark QCD condensate Q α Q β ∝ δ αβ Λ 3 dQCD . This spontaneous symmetry breaking delivers 8 Nambu-Goldstone bosons, π D1 , . . . , π D8 , which will become pNGBs once we switch on m Q . We will consider the case where m Q Λ dQCD so that the pNGBs are parametrically lighter than the rest of particles in the spectrum. With the exception of the lightest baryonic bound states carrying a conserved dark baryon number, the rest of the particles of the spectrum will undergo fast decays to dark pions. Therefore, it is a reasonable approximation to just consider such dark pions. In the case at hand where n d = 3 = N d , one can write down an effective theory for the resulting eight pNGBs along the lines of the QCD case of pions and kaons. Using the basis of Gell-Mann matrices, λ a , a = 1, . . . , 8, one can write The corresponding Goldstone matrix transforming nonlinearly under SU (3) dL ⊗ SU (3) dR and linearly under SU (3) dV can be written as where f d is the dark pion decay constant, in principle a free parameter. In analogy to chiral perturbation theory (ChPT) for QCD, the Lagrangian describing the dark pions in the absence of interaction with the SM is given by where B d is a constant related to the dark pion mass. More precisely As one can see, they are all degenerate in mass, but small splittings will be induced by their interactions with the SM. Such radiative corrections will define new mass eigenstates with π D3 and π D8 unchanged. If the dark pions are light enough, m π D 4πf π , their decays are better described by ChPT for the SM quarks. The part containing only SM fields is described by whereas the interaction with the dark QCD sector is described by where the projection matrices c αβ and c ij are defined as c mn αβ = δ m α δ β n , α, β = 1, 2, 3, c mn ij = δ m i δ n j , i, j = 1, being zero otherwise. For larger dark pion masses, one can use quark-hadron duality [90,91], with the relevant Lagrangian being For the benchmark model at hand, only π D3 , π D8 and π (1,2) D decay at tree-level. The rest of dark pions, which decay through loop-induced processes, are therefore longlived. For this reason, we focus on the first dark pions, and in particular, on the two real ones π D3 and π D8 . 9 In particular, the mixing terms in eqs. (A16) and (A18) involving π D3 and π D8 read and L mix ⊃ − f d 2m 2 X a=3,8 αβ κ αi κ * βj (λ a ) αβ ∂ µ π Da (ū Ri γ µ u Rj ), respectively. Comparing these equations with the mixing terms present in eqs. (1) and (9) for a = π D3 and π D8 , respectively, we obtain f a = m 2 X /f d and Another motivation for the ALPs considered here corresponds to what is generically known by the name of flavons or familions (see [87,[109][110][111][112][113] and e.g. [114][115][116][117][118][119][120][121][122][123] for more recent implementations), pNGBs of some spontaneously broken flavour symmetry, which may be anomalous, and that generically feature flavour-violating couplings to quarks or leptons. One particular setup leading to the scenario we have in mind, i.e., the effective Lagrangian (1) in addition to c H = 0, is given by FN models when only RH up-quarks have non-zero charges. Specifically, one considers a global U (1) flavour symmetry, spontaneously broken by the vacuum expectation value of some extra scalar S = f a , where and has charge −1 under this new U (1). If only u Ri are charged under such global symmetry, with charges n u i , Yukawa couplings for up-quarks will be higherdimensional where one typically assumes that f a < Λ. At the end of the day, such term in the Lagrangian will generate interactions like (4) where Y u is the effective up Yukawa matrix, If we assume that (y u ) ij = O(1) and take f a /Λ = ∼ m c /m t , we can get the correct up quark masses by choosing n u = (2, 1, 0) since In this case, we can diagonalize Y u by making with This leads, after going to the basis where Y u is diagonal to (B9) N = p, n, which read L int ⊃ 1 12 ∂ µ a f a (c u R ) 11 S − 4D (nγ µ γ 5 n) + 2D + 6F + S (pγ µ γ 5 p) + 6pγ µ p (C4) The last term in the equation above is a total derivative which can be neglected. In this case c ann = (c u R ) 11 1 6 S − 2 3 D = (c u R ) 11 (−0.51 ± 0.03) . (C7)