Exploring new possibilities to discover a light pseudo-scalar at LHCb

We study the possibility of observing a light pseudo-scalar a at LHCb. We target the mass region 2.5GeV≲ma≲60GeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.5\,\mathrm{GeV}\lesssim m_a\lesssim 60\,\mathrm{GeV}$$\end{document} and various decay channels, some of which have never been considered before: muon pairs, tau pairs, D meson pairs, and di-photon. We interpret the results in the context of models of 4D Composite Higgs and Partial Compositeness in particular.


Introduction
The search for resonantly produced particles Beyond the Standard Model (BSM) is a high-priority goal at the LHC. Many searches focus on resonant production and decays into di-bosons (see, e.g., [1,2]) or di-fermions (see, e.g., [3][4][5]) in the high-mass region. This has led the LHC experiments to push the exclusion to mass scales, in most cases, well above the TeV. Yet, new states with small masses may still be allowed and lie in unconstrained oases of the BSM parameter space. One example is provided by electrically neutral, colorless, scalars with a mass below the Z pole and above the heavy meson mass scales [6]. Rare meson decays, in fact, provide an additional class of strong bounds [7][8][9][10]. At the LHC, low masses are mainly constrained by di-muon [11][12][13][14][15] and di-photon searches [16][17][18][19][20][21]. A light spin-0 state could emerge in many BSM scenarios like supersymmetry, Higgs sector extensions, and models based on composite dynamics. Specifically, a pseudo-scalar a with a mass much lighter than the BSM scale is naturally realized as pseudo-Nambu-Goldstone boson (pNGB) associated with the spontaneous breaking of an approximate global symmetry. A time-honored example is provided by axions emerging from the Peccei-Quinn solution to the strong CP problem of QCD [22][23][24][25]. Other models of different nature featuring a light pNGB fall under the generic class of models of axion-like particles (ALPs) [26,27]. ALPs can be found in supersymmetry [28], models of composite electroweak symmetry breaking [29][30][31], and models with extended scalar sectors, like multiple Higgs doublet models [32,33] (including 2HDMs), type-II see-saw models for neutrino masses [34], and models with custodial triplets [35,36]. Among them, we focus specifically on composite Higgs models with an underlying fermionic UV description [37]. In particular, a light ALP is ubiquitous [6,38,39] in models with two species of confining fermions, needed to implement top partial compositeness [40,41].
The physics of ALPs can be encoded in a generic effective Lagrangian at low energy. The pNGB nature of the pseudo-scalar bears additional information on its coupling to Standard Model (SM) fermions (via derivative interactions that yield couplings proportional to the fermion masses) and gauge bosons (via Wess-Zumino-Witten terms). Moreover, in composite models, the coefficients are related to each other [38], as they emerge from the same underlying dynamics, and thus the branching ratios of a into SM particles can be correlated and predicted. Hence, con-trary to a generic ALP scenario, detection prospects in different decay channels can be compared. A pseudo-scalar pNGB with a mass below the Z pole is expected to dominantly decay into the heaviest accessible fermion pairs (bb, cc, τ + τ − ) or gluons gg (i.e. light hadrons), for which no low-mass LHC searches are available, so far (see Ref. [39] for a proposal of a low-mass di-tau search). The branching ratios into the experimentally tested μ + μ − and γ γ channels are typically small. In this mass range, the production mode at LHC is completely dominated by gluon fusion and this is the only production mode considered in this work. For studies of ALP production at lepton colliders see e.g. Refs. [42][43][44][45].
In this article, we present the first study of LHCb prospects to observe a composite ALP a in the cc and τ + τ − channels. For comparison with existing bounds, we also re-interpret searches in the μ + μ − and projections for the γ γ channel. As benchmarks, we consider the 12 composite Higgs models (M1-M12) defined in Refs. [6,38,41]. However, results are also presented model-independently and can be applied to any other light pseudo-scalar model. At the LHC, the LHCb detector [46] is a forward spectrometer, whose special features make it appropriate for the types of signatures described in this work [47]. This includes the capability of triggering on soft objects, excellent vertex reconstruction that is useful to distinguish shorter lifetime objects, such as τ leptons, and very good invariant mass resolution that provides advantages for the discrimination against large continuous backgrounds. This last feature is crucial in the case of cc. LHCb is currently undergoing an upgrade [48][49][50], after which it is expected to collect 15 fb −1 over the next three years. Overall, the experiment will collect 300 fb −1 in its whole lifetime [51,52].
The paper is organized as follows. In Sect. 2 we briefly introduce the effective Lagrangian and the benchmark models used in this article. In Sect. 3 we present the recast and projections for the existing di-muon searches. In Sect. 4 we compare to the reach obtainable in the di-photon final state. In the following two sections, 5 and 6, we describe in detail new search proposals for final states containing taus and D mesons. Finally, we offer our conclusion and summary plots in Sect. 7.

Model and simulation
The phenomenology of a light ALP a can be generically described by the following effective Lagrangian Table 1 The numerical values of the coefficients in the lagrangian (1) for the 12 models considered. Only the sum K γ = K W + K B is relevant for the mass range considered in this work. The coupling to the top quark C t can take a discrete set of values depending on the spurionic charge assignments. In this work we chose the value that leads to the largest constructive interference to the gluon coupling and thus the largest production cross section via gluon fusion where ψ are the SM fermions, F μν the SM gauge field strengths (F = G a , W i , B),F μν = 1 2 μνρσ F ρσ , and f the ALP decay constant. In models addressing the hierarchy problem, like composite Higgs ones, the scale f is typically assumed to be in the TeV range. For recent systematic studies of the dynamics of the Lagrangian above see Refs. [53][54][55][56][57]. Note that all couplings are considered of order 1 for generic ALP scenarios.
The presence of an explicit mass for the ALP evades the usual constraints that exclude the Peccei-Quinn-Weinberg-Wilczek (PQWW) axion [22][23][24][25] for a TeV scale f but also precludes its use to solve the strong CP problem. Conversely, this scenario arises naturally in models of composite Higgs with top partial compositeness, emerging from an underlying gauge theory with fermions [41]. An ALP state, potentially light, is an unavoidable byproduct of the global symmetries broken by the condensates [6,38], which generate both a composite Higgs and composite top partners.
A main feature of this class of composite ALPs is that the coefficients in the effective Lagrangian can be computed in terms of the underlying theory, hence rendering the theory highly predictive. In this article, we follow the set of benchmark models M1-M12 that are defined in Ref. [6]. They yield predictions for the couplings C ψ , K g , K W , K B summarized in Table 1. We emphasize that the couplings are computed from the underlying theory and not arbitrarily chosen. The couplings to gauge bosons K g,W,B are generated by anomalies, and do not include the effects of loops of SM fermions, which are computed separately. In particular, the coupling to the photon K γ = K W + K B can take values between −13.7 and 20.3 yielding very different branching ratios in this channel. For a detailed discussion of the models and the derivation of the coupling constants, we refer to Ref. [6]. In the following, we will only summarize the aspects relevant for the LHCb study presented here.
Contrary to the generic scenario with all order 1 couplings, the composite ALP scenarios M1-M12 offer cases where specific couplings could be substantially enhanced or suppressed. This opens up the possibility of discovery in what would be generically considered sub-leading channels or vice versa. As a concrete example, one finds models like M3 and M8, where the photon channel is suppressed. This feature arises from the specific electroweak couplings of the confining fermions. In fact, the main difference of Eq. (1) from similar Lagrangians arising in the context of 2HDMs is that, for the models at hand, the K g,W,B constants denote the anomalous contributions from the (confined) hyperquarks of the UV theory and are thus non-zero even before integrating out the heavy SM quarks (mainly the top and bottom, for the mass ranges we consider in this study). Note that such terms could also be present in 2HDMs and other extensions of the scalar sector of the SM if heavy non-SM fermions are included.
In the remainder of the paper, we focus on the four most promising signatures for the search of composite ALPs at LHCb, with special focus on the role it can play in the coming runs. We consider, in turn, the decay modes: The most studied decay channel, and the one likely to be dominant under generic assumptions on the couplings, is the di-muon channel. Here, we offer a straightforward recasting of the existing searches targeting 2HDMs [11][12][13][14][15] and convert the bounds on the mixing angle among Higgses to those of the tuning parameter v/ f , where v = 246 GeV is the Higgs vacuum expectation value.
In the di-photon channel, we use the projections in Ref. [16], (obtained from inclusive diphoton cross-section measurements imposing that the signal events are less than the total measured events plus twice their uncertainty), which can be treated similarly to the di-muon case. In our models, this channel is expected to give relatively weak exclusion bounds.
The di-tau and charm channels have not been considered before in the LHCb context. The di-tau channel can potentially cover a significant mass range (from 14 to 40 GeV). It benefits from a branching ratio enhanced by a factor (m τ /m μ ) 2 ≈ 283 compared to the di-muon (for C τ ≈ C μ ) but suffers from the presence of neutrinos and hadrons in the  Table 1 for f = v. The cross sections scale as (v/ f ) 2 ) and are computed at NNLO in QCD using the HIGLU program [62,63]. The mass range starts at 2.5 GeV to stay away from the non-perturbative region (see Ref. [64]) final states. The a → cc decay mode is relevant in the mass range 3.8 GeV m a 6 GeV. (For m a < m J/ψ the nonobservation of the J/ψ → a γ process puts strong bounds on f [58].) Besides the recast of existing bounds, we will offer projected LHCb reaches for all fours channels for integrated luminosities of 15 fb −1 and 300 fb −1 .
For the signal simulations, we compute the ALP total production cross sections σ ( pp → a), dominated by the gluonfusion channel, with the HIGLU program [59] at NNLO in QCD, using the NNPDF 3.1 containing LHCb data [60] for 14 TeV in the pp center-of-mass. The numerical results are shown in Fig. 1 for each of the 12 models for f = v. The cross sections always scale as (v/ f ) 2 . In the plot we use renormalization and factorization scales μ F = μ R = m a . It should be noted that at low masses (few GeVs) a large scale dependence is present and the predictions start to become unreliable. In Fig. 1 we only show the central value obtained for each model and refer to Ref. [61] for a discussion of the expected error estimated by varying μ F and μ R . We use these values only for limit setting. The branching ratios of a into the four decay channels (2) are shown in Fig. 2. We use analytical expressions for the partial widths into μ + μ − , τ + τ − , cc, bb, γγ [6]. For the a → gg decay channel we use instead the HIGLU program [62,63].

Recast and projections for a → μ + μ −
If the ALP couples to the muon with a typical strength C μ ≈ 1, we expect the muon channel to give strong bounds for a wide range of masses. Searches in this channel have already been performed by various collaborations [11][12][13][14][15] and we start by presenting a recast of these results. We use the summary plot in fig. 10 of Ref. [11], presenting upper limits at  (2) for the 12 models. Note that in these models the ALP is always promptly decaying with a narrow width, ranging, in this mass region, from few keV to few MeV for f = 1 TeV. The HIGLU program has been used for the computation of the a → gg decay channel [62,63] 90% confidence level on the mixing angle sin θ H between the pseudo-scalar of a 2HDM and the imaginary component of a complex singlet. The plot is obtained from the total di-muon cross-section of the previous searches [11][12][13][14][15] and tests the type IV 2HDM model of [61,65] at tan β = 0.5. The bounds are set on the mixing angle sin θ H as a function of the ALP mass, ranging between 1-60 GeV.
The recast is easily implemented thanks to the following two observations. The first one is that, for all models, the narrow width approximation holds very well and thus σ ( Note however that in these models the ALP is always promptly decaying since the width ranges from few keV to few MeV for f = 1 TeV. The second observation is that all σ ( p p → a) are proportional to sin 2 θ H in the 2HDM and proportional to v 2 / f 2 in the composite ALP models of Eq. (1), while all the branching ratios are independent on these quantities. Thus, denoting the cross-section of the 2HDM model with sin θ H = 1 byσ 2HDM ( p p → a) and, similarly, the cross-section of the models of Table 1 with f = v byσ Mi ( p p → a), the bounds on v/ f are obtained from the bounds on sin θ H in Ref. [11] by the simple rescaling Using this procedure one could also recast the LHCb search [11] to obtain bounds on any of the four types of 2HDMs (I, II, III, and IV) for any value of tan β. If one were content with working at leading order, employing the detailed balance σ ( p p → a) ∝ Γ (a → gg) one could replace the cross-sections in Eq. (3) with the partial width and obtain an analytic formula valid to within few % in the mass region m a > 15 GeV. The relevant widths are listed in Ref. [65] for the 2HDM and in Ref. [6] for the ALP models. For the total width we sum over the channels μ + μ − , τ + τ − , cc, bb, gg, γγ. For the a → gg decay, giving the full hadronic decay at lower masses, we used HIGLU and compared the results with analytic ones including one-loop renormalization of the quark masses and the gauge couplings as well as the finite part of the QCD corrections [62,63]. We find good agreement between the two. More precisely, the mean deviation between the two estimates, averaging over the 12 models, is Bounds on v/ f as a function of m a for the di-muon channel at 90% C.L. This is a recasting of Fig. 10 in Ref. [11] using the envelope of the exclusion curves from the searches [11][12][13][14][15]. In order to smooth out some of the variability of the original exclusion curve we performed a moving average over the 10 nearby mass point for each point on the axis 20%, 4.1%, 1.5%, and 0.56% for ALP masses of 2.5, 5, 30, and 60 GeV respectively. The non-perturbative aspects of ALP decay into light hadrons become crucial for lighter ALPs and are discussed in Ref. [64]. The cross-sections are also computed numerically at NNLO with HIGLU as described in Sect. 2.
We reiterate that the important difference between 2HDMs and Eq. (1) is the presence in the latter of contact terms to gluons and photons coming from the anomaly of the hyperfermions, which are absent in 2HDMs and lead to an enhancement of the gg coupling strength. Just like for the value of tan β in the 2HDM, we need to fix the coupling strengths C ψ to the fermions in our models. These are given in Table 1.
The coefficients for all fermions other than the top quark are fixed in the underlying theory as explained in Refs. [6,38]. The variability comes from different discrete choices for the coupling of the top quark due to different spurion charge assignments. Throughout the paper we present the bounds arising from the choice of C t giving the largest effective coupling between the ALP and the gluon. This leads to the strongest bounds and sensitivities for three of the four channels in Eq. (2) due to the constructive interference between the hyper-fermion anomaly and the top quark coupling. The exception is the di-photon channel, discussed in more detail in Sect. 4. Figure 3 shows the bounds on v/ f for the 12 benchmark models that arise from the rescaling of the di-muon bounds of Ref. [11] as described in Eq. (3). For each ALP mass hypothesis we use the strongest bound of the four experimental analyses: BaBar [12], CMS Run 1 [13], CMS Run 2 [14], LHCb Run 1 [15], and LHCb Run 2 [11]. These bounds are the strongest constraints on v/ f for these models in the considered m a mass range to-date, showing that this channel is the most sensitive one under the assumption that the ALP couples to the muon with standard strength. However, for a "muon-phobic" ALP the other channels become relevant and should be considered in order to broaden the reach.
We conclude this section by presenting in Fig. 4 the projections for the di-muon channel obtained by rescaling the LHCb results [11] (Run 2 only, top panel) from 5.1 fb −1 to 15 fb −1 (middle panel), and 300 fb −1 (bottom panel). We expect the results to be dominated by statistics. Hence, noting that S/ √ B ∝ √ L and that S ∝ (v/ f ) 2 , the exclusion in v/ f scales like L −1/4 . Note that the projections do not account for the removal of the first trigger level at LHCb [48], based on hardware, which is expected to happen from Run 3 onward. While doing this would yield more stringent exclusions, we choose to extrapolate from the existing LHCb results, based on data, which provides more realistic estimates of the background.

Projections for a → γ γ
We continue by examining the di-gamma decay channel. We restrict ourselves to the low mass region 2.5 GeV < m a < 20 GeV where the strongest bounds derive from the analysis in Ref. [16] of the data from ATLAS, CMS, LHCb, and BaBar [17][18][19][20][21]. The exclusions [16] translate into fairly weak current bounds for our models and we thus simply present the projected LHCb bounds obtained from Ref. [66,67] in the same spirit as for the di-muon channel.
Reference [66] investigated the γ γ channel for a fully fermiophobic model (C ψ = 0) whose anomaly coefficients can be written in our notation (1) as K g = 10 and K W + K W = 80/3. The same reasoning that yields Eq. (3) now gives ("fph" stands for fermiophobic) the barred cross-sections denoting the values with f = v.
In Fig. 5 we show the projected sensitivity for the 12 models for the LHCb run with integrated luminosities of 15 fb −1 (top) and 300 fb −1 (bottom) respectively. For the calculation of the partial width of a → γ γ , we take into account the one-loop corrections from t, b, c, and τ loops, which yield sizable corrections, in particular for low m a . The remaining calculation is performed exactly as explained in the di-muon section.
In Fig. 5 we still choose the value of C t leading to the largest effective coupling of the ALP to the gluon and thus to the largest production cross section (within each benchmark model). However, contrary to all other channels, the C t with largest ALP production cross section does not necessarily lead to the strongest exclusion bound in the di-photon channel. This is because, contrary to the fermionic decay channels, the decay width Γ (a → γ γ ) depends strongly on the value of C t , as explained in Ref. [38].
Specifically, the top loop correction to the effective coupling of the gluon and the photon are, respectively  Fig. 1 of Ref. [66] using (4) with values for the couplings in Table 1.
Note that for L = 15 fb −1 many models are basically unconstrained where K g , K W , K B are the anomaly coefficients in Table 1 from the hyper-fermions, p(τ ) = τ arctan 2 ( √ 1/(τ − 1)) and the dots represent contributions from the lighter SM fermions. For all models one always has K g < 0, and C t ranges in an interval containing positive and negative values. The largest |K g,eff. | 2 is thus attained by picking the spurion charges giving the largest (positive) C t for each model. On the contrary, K W + K B can have positive or negative values and there can be destructive interference if one picks the largest C t .
For a consistent comparison with (projected) bounds in the di-muon, di-tau, and di-charm channel, we use the same values of C t for the analysis of the di-photon search, but it should be noted that other choices of C t can lead to altered sensitivity in the di-photon channel, up to an order of magnitude.

New search in a → τ + τ −
In this section we show the potential of the LHCb experiment to search for a → τ + τ − decays. Despite the presence of neutrinos and the fact of being a non-hermetic spectrometer, LHCb has already shown its potential to search for decay modes including τ leptons in the final state. Notable examples include τ leptons produced either from a displaced low-mass vertex (semileptonic decays of B mesons [68] or pairs from the leptonic decay of a B 0 s meson [69], where τ leptons are reconstructed using three charged pions in the final state), or from a prompt high-mass vertex (decay of a Z boson into τ + τ − [70], where several combinations of the decay modes of both τ leptons -hadronic, semileptonic, and fully leptonic -are explored).
The signal topology described in this paper is challenging for LHCb, due to the fact that both τ leptons are produced from the decay of a prompt object, which has a relatively low invariant-mass. However, the excellent capabilities of the detector to reconstruct soft objects in the final state help suppressing most of the dominant background components that would pollute our signal. A previous proposal for a a → τ + τ − search at CMS and ATLAS can be found in Ref. [39].

Simulation and analysis strategy
The signal process is defined by the production of a prompt light pseudo-scalar in proton-proton collisions at a center-ofmass energy of 14 TeV, decaying into a pair of tau leptons with opposite charges. It is simulated with the Pythia 8.305 program [71] with fully spin-correlated tau decays [72].
A signal fiducial region is defined in order to account for differences in the production mechanism between this simplified Pythia model and the NLO model described in more details in Sect. 6. This region is defined by selecting pseudoscalar Pythia objects with a pseudo-rapidity between 2 and 4.5 and a transverse momentum p T between 15 and 150 GeV. We also impose requirements on the two τ leptons: a pseudorapidity between 1.5 and 5, at least one τ with p T > 7.5 GeV and the second tau with p T > 5 GeV. These requirements are also part of the selection imposed to the reconstructed objects, as described in the following paragraphs. We checked that any leak in selected data from events not in the fiducial region is negligible.
Charge-conjugate final states are left understood.
Three main SM background processes are expected to contaminate the signal selection: -QCD multijet production, pp → j j with j standing for light quarks, gluons, and charm quarks. -QCD heavy-flavor pp → bb production. Most of the background leptons are expected to originate from bhadron decays. -Drell-Yan pp → τ + τ − production.
Other background components (such as other Drell-Yan productions and Z /W boson pair production with leptonic and semileptonic decays) were found to be negligible. All background processes have been simulated with Pythia8.
The QCD multijet background is challenging to estimate and is expected to be dominant for the fully hadronic and semileptonic channels. On the other hand, its contribution to the fully leptonic mode is estimated to be at most 10% of the dominant bb background. We also expect the fully leptonic channel to be the most sensitive one, see Appendix A for more details, hence we will neglect the other channels in computing the limits. We therefore only consider the two dominant backgrounds sources in the analysis of the fully leptonic mode, namely, heavy-flavor bb and Drell-Yan productions.
We limit our study to the mass region above 14 GeV and below 40 GeV. Below 14 GeV, the QCD background becomes unacceptably large and the Υ resonances are present, severely limiting the sensitivity at LHCb. Conversely, above 40 GeV the signal efficiency becomes compromised due to acceptance limitations.

Analysis of the eμ mode
All charged stable tracks under consideration (muons, electrons, pions, kaons, and protons) are required to be inside the LHCb pseudo-rapidity acceptance, 2 < η < 5. A requirement of having the particles produced within the Vertex Locator (VELO) region is also imposed, that is, a minimum p T of 0.5 GeV, with V r < 30 mm and V z < 200 mm, where (V z , V r ) denotes the spatial position of their production vertices, expressed in cylindrical coordinates. Kaons, pions, and protons are later used only to define the quantities used to define the electron and muon track isolation.
Electrons and muons are required to have a p T greater than 3.5 GeV, a minimum energy of 10 GeV, a minimum impact parameter (IP, defined below) of 0.01 mm for muons (and 0.03 mm for electrons), and to be well isolated from other tracks in the event. For this purpose, a quantity I to measure the isolation of a track is defined as the fraction of its p T over the sum of the p T of all stable charged tracks (muons, electrons, pions, kaons, and protons) inside a cone of a certain ΔR 2 = Δφ 2 + Δη 2 value, built around the track of interest. We require at least one of the muons or electrons to be well isolated, imposing a tight cut to one of them, this is, max(I μ , I e ) > 0.99, for ΔR 2 = 0.05. Note that the discrimination provided by isolation tends to be overestimated in simulation with respect to data. Since this is an effect that is hard to estimate, we use the discrimination provided by our simulation, but we acknowledge this is a limitation of our study.
ALP a candidates are reconstructed by summing the 4momenta of the selected (e, μ) pair. The reconstructed a is required to be prompt with a maximum distance of flight of 1 mm, to have an IP smaller than 0.2 mm, a pseudo-rapidity between 2 and 4.5, and a transverse momentum between 15 and 150 GeV (otherwise, the signal would be polluted by low p T QCD background contributions). The muon-electron pairs from the decay of a are required to have a maximum distance of closest approach (DOCA) of 0.4 mm. Each lepton must have a minimum p T > 5 GeV and at least one of them must have p T > 7.5 GeV. The DOCA is the minimum distance between two different trajectories, defined as |V × ρ|/|ρ| where V = V 1 − V 2 and ρ = p 1 × p 2 . Here, V j and p j are a vertex 3D position and the tri-momentum associated to a track j ( j = 1, 2), respectively. With the same vertex and tri-momentum definitions per track, the IP of a track j is defined as Moreover, for each signal mass hypothesis, the invariant mass of the eμ system is required to be in a range that enhances the signal over background ratio. The ranges are shown in Table 2 for each ALP mass hypothesis.

Computation of efficiencies and bounds
Signal, bb, and DY background efficiencies, AL P , bb , and DY respectively, are obtained for several a mass hypotheses, and shown in Tables 2 and 3. The signal efficiency AL P is the product of two contributions: the first one, computed with the NLO model, is the ratio of the cross section in the signal fiducial region and the inclusive cross section pp → a → τ + τ − ; the second one, computed with the Pythia8 LO model, is the efficiency of the fully leptonic analysis in the fiducial region.
For each mass hypothesis, the number of signal and background events are   Table 3 Signal efficiencies for the τ + τ − channel in the eμ reconstruction mode, considering the NLO production model in the full acceptance. Mass window requirements are imposed on top of the selection, as described in Table 2  where σ is the production cross-sections, B the branching fractions of a and τ , and L the integrated luminosity. The bb production cross-section, σ ( pp → bb) = (562 ± 82) × 10 9 fb, is taken from Ref. [73]. The Drell-Yan cross-section for the different mass windows are taken from Ref. [74] to be σ ( pp → τ + τ − ) = (4.494 ± 0.237) × 10 6 fb, and the values of the branching fractions of the different τ decay modes are taken from Ref. [75].
Apart from discriminating against the background, our tight selection is defined to maximize the reconstruction and trigger efficiencies at LHCb. The removal of the first trigger level at LHCb [48], mentioned above, makes this assumption more realistic. Therefore, we assume this efficiency to be 100%. Moreover, we neglect experimental resolution effects that would affect the computation of the invariant masses, since the inaccuracy on this is dominated by the presence of neutrinos. The signal and background di-lepton invariant mass m(eμ) distributions are shown in Fig. 6, for L = 300 fb −1 , and model M1 with v/ f = 0.1 for the signal.
We also neglect systematic uncertainties and provide expected limits at 90% C.L based purely on statistical uncer-   [76,77]. An actual experimental search would be needed to control systematic effects.
In Fig. 7 we show the projected bounds on the signal cross section σ ( pp → a) × B(a → τ + τ − ) as a function of the ALP mass, assuming an integrated luminosity of L = 15 fb −1 (light blue) and L = 300 fb −1 (dark blue). Signal In Fig. 8 we show the projected bounds on v/ f for the 12 benchmark models, assuming an integrated luminosity of L = 15 fb −1 (top) and L = 300 fb −1 (bottom).

New search with D mesons targeting a → cc
The cc channel is only relevant for a small range of low ALP masses, 3.8 GeV m a 6 GeV. It is especially motivated in scenarios where the muon obtains its mass from a different mechanism and a couples only to quarks of the up type, similarly to the scenarios in Ref. [58], although we will not consider flavor violating couplings (which have also been studied recently in Refs. [78,79]) nor the mass range relevant for J/ψ decay.
To estimate the LHCb sensitivity in this channel we perform a novel dedicated analysis. We generate signal events using MG5_aMC@NLO [80] with the Higgs Characterization model [81] and pass them through Pythia8 [71] for showering, hadronization, and decays. For the purpose of computing efficiencies, the signal event generation is performed at NLO in QCD.
The background is expected to be fully dominated by the QCD production of cc. We use the total background cross section σ B (cc) = 7.1 ± 3.4 mb [82] and simulate cc events with Pythia8 at LO. The idea is that slightly above threshold (m a ≥ 3.8 GeV) the ALP decay leads to a fully reconstructable D + D − pair whose invariant mass fits into a very narrow bin (≈ 40 MeV) allowing to overcome the huge background.
The LHCb capabilities in terms of particle identification allow identifying D ± with nearly 100% efficiency if they decay into K − 2π + (K + 2π − ) [83], so we select events with at least one D + and one D − each decaying into this mode. The rate of events with at least one D + D − pair compared to the total cc production is denoted by f cc→D + D − . We select events in which all the six decay products K ± and 2 × π ± are within LHCb coverage 2 < η < 5 and p T > 0.25 GeV. The corresponding acceptance is denoted by f Acc . The values obtained for the background QCD cc productions are The corresponding values for the signal for different values of m a are shown in Table 4. Lastly, we require the events to fulfill m(D + D − ) = m a ±20 MeV, given the LHCb high invariant mass resolution reconstruction of the D + D − system. This last cut, denoted by f mass is almost fully efficient for the signal in the low mass region and allows to dramatically reduce the background rates. This resolution corresponds to approximately ±2σ , where σ ∼ 9 MeV is the D + D − invariant mass resolution, taken from the LHCb measurement of the B 0 → D + D − decay [84]. The invariant masses of the D + D − system with a Gaussian smearing according to this resolution are shown, for model M1 with v/ f = 1 and with L = 300 fb −1 , in Fig. 9. The cumulative efficiencies of this final selection, are shown in Table 4. The signal efficiency drops rapidly for m a 5.5 GeV due to the opening of other decay channels. Figure 10 shows the resulting bound, obtained by the CLs method at 90% C.L., on the signal cross section σ ( pp → a) × B(a → cc) as a function of the ALP mass, assuming an integrated luminosity of L = 15 fb −1 (light blue) and L = 300 fb −1 (dark blue). Following Ref. [85], systematic uncertainties are expected to be below 0.01%, since the background yields can be also extrapolated here from the invariant mass side bands, and therefore are neglected. For reference, in Fig. 10 (top) we show the central prediction for the signal cross section σ ( pp → a) × B(a → cc) of the benchmark models M1, M4, M9, and M11 with v/ f = 1 Table 4 Cumulative efficiencies (in %) for both signal and background. The background fragmentation fraction f B cc→D + D − and acceptance f B Acc are independent on the mass points and are given in (9).
mass . The total number of signal events are given by S = σ (pp → a) × B(a → cc) × B(D → K ππ) 2 × S × L and similarly for the background. B(D → K ππ) = 9.38% is the branching ratio D + → K − 2π + (and D − → K + 2π − ), and L the integrated total luminosity   Figure 11 shows the projected bounds on v/ f in the 12 benchmark models, assuming an integrated luminosity of L = 15 fb −1 (top) and L = 300 fb −1 (bottom), which result from simple scaling of the projected bounds in Fig. 10.
Having discussed the decay into charm pairs as a possible discovery channel, one may wonder why not considering bottom pairs as well. Unfortunately applying the same strategy to the a → bb channel is not viable, in spite of the larger ALP branching ratio B(a → bb) ≈ 21 × B(a → cc), due

Conclusion
Light pseudo-scalar particles are present in many Standard Model extensions. In particular pseudo-Nambu-Goldstone bosons resulting from a spontaneously broken underlying global symmetry may successfully evade current searches even if they are as light as a few GeV. Their couplings to gauge bosons are suppressed as they arise at loop level or through anomaly contributions, while their couplings to fermions are proportional to the fermion masses.
In this article, we provided a first study for the prospects of detecting a light pseudo-scalar at LHCb in τ + τ − and in cc (D + D − ) channels, with √ s = 14 TeV and luminosities of 15 fb −1 or 300 fb −1 . We also compared these new channels to the projected reach of existing searches in μ + μ − and prospects in di-photons (the main production mode in all these channels is gluon fusion). As benchmarks, we use 12 models of composite Higgs with top partial compositeness, which lead to calculable couplings of the ALP to the SM fermions and gauge bosons. The results can, however, be For the a → τ + τ − channel discussed in Sect. 5 we designed an analysis strategy targeting prompt a → τ + τ − with subsequent di-tau decays in four categories: fully hadronic, semileptonic with a muon, semileptonic with an electron, and fully leptonic (eμ). However, for the computation of the limits we only considered the fully leptonic mode, which is found to be highly dominant over the other three decay categories. We find the exclusion reach on σ ( pp → a) × B a → τ + τ − given in Fig. 7 for a luminosity of 15 fb −1 and 300 fb −1 . A dedicated study on the signal efficiencies of the hadronic and semileptonic modes is discussed in Appendix A.
For the a → cc channel discussed in Sect. 6, we focused on the exclusive final state D + D − , which is fully reconstructable at LHCb. This final state is only relevant for masses right above the threshold of 3.8 GeV. The limits quickly deteriorate above 5.5 GeV, yielding limits on σ ( pp → a) × B (a → cc) in Fig. 10.
Comparing prospects to find a pseudo-scalar a in the newly proposed a → τ + τ − or a → cc channels to the already studied a → μ + μ − or a → γ γ channels is inher-  ently model-dependent, as the comparison depends on the branching ratios of a. For the 12 models under consideration we focused on the top couplings that maximise the production cross section in gluon fusion and the overall sensitivity in the fermionic channels at the price of reducing the sensitivity to di-photon, for some models.
As summary, in Figs. 12 and 13 we show the projected limits for a luminosity of 300 fb −1 for the 12 models, where models with the same global symmetries at low energy are shown in the same panel (with the exception of M5 and M12 in the last one. See Ref. [38] for details). The plots reveal that the muon channel remains the dominant one, but the new ones, particularly the di-tau, provide comparable and complementary information albeit in a more limited mass range. The photon channel is always minor, even though it provides unique limits in the mass ranges not covered by the muons due to the presence of large background from J/ψ and Υ resonances.
We conclude that LHCb has excellent prospects to investigate light di-tau resonances. For models with an identical coupling to μ and τ (C μ = C τ in Eq. (1)), the first feasibility study presented here promises an exclusion range which is comparable to the well-established and highly optimized dimuon resonance searches. The di-charm resonance search is applicable only in a small mass range near the D + D − thresh-old, and yields weaker bounds than the di-muon search (under the assumption C c = C μ ), but offers bounds which can partially cover the gap left in the di-muon search near the J/ψ resonance.
Finally, we emphasize that searches in the di-tau and dicharm channel are of course interesting in their own right.
The assumption of uniform pseudo-scalar fermion coupling is theoretically well motivated in these models but not guaranteed in general (see e.g. Refs. [58,78] for counter examples).

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This sensitivity study is solely based on Pythia and MadGraph simulations, which are easily reproducible with the information in the paper.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Hadronic and semileptonic τ modes
As mentioned in Sect. 5, we have also conducted a study of the hadronic and semileptonic modes of τ decays.
The QCD component of the background for these modes, dominated by cc pairs and jets produced from gluons and light quarks, is overwhelming. Given the fact that our simulation and reconstruction framework becomes substantially slower when the τ 3-prong decay modes are involved (due to the additional selections and due to the vertex reconstruction of the three pions), it is unfeasible to simulate enough QCD background events that pass our full selection. Instead, a proper study of the h 3 h 3 , h 3 μ, and h 3 e categories should be done by using a minimum-bias LHCb dataset of protonproton collisions, to which we have no access for this study.
Furthermore, our limits are fully dominated by the eμ decay mode. We have tested this by computing the bound in σ ( pp → a) × B(a → τ + τ − ) with the CL s method, using the four categories h 3 h 3 , h 3 μ, h 3 e, and eμ. As a conservative check, we have scaled the bb contribution by a large factor to account for the missing sources of background in the QCD component as previously mentioned, leading to almost negligible changes in the combination.
For all the above reasons, the bounds in Sect. 5 are obtained using only the eμ channel. Nevertheless, we believe it is worth presenting in this appendix the selection and reconstruction procedure, as well as the signal efficiencies, of the hadronic and semileptonic categories to provide a useful input for potential future studies for these modes using LHCb data.
In discussing the semileptonic channel, pions are required to have a minimum p T of 1 GeV (just as for the electrons and muons), while for the fully hadronic channel this requirement Table 5 Signal efficiencies for the hadronic and semileptonic modes of τ + τ − , considering the NLO production model in the full acceptance. Mass window requirements are imposed on top of the selection, as described in the caption of Table 6 Mass (GeV) h3h3  is loosened to 0.5 GeV. In all cases a minimum IP of 0.01 mm, and a minimum momentum of 2 GeV are required. The pions are then combined for all these channels considering all possible three-body combinations of appropriate charge. The reconstruction procedure goes as follows: the threepion combination is required to have an invariant mass below 1.7 GeV, a minimum p T of 10 GeV (2.5 GeV for the h 3 h 3 category), an IP smaller than 0.2 mm (0.1 mm for the h 3 h 3 category), and a corrected mass between 1.2 GeV and 2.5 GeV. In order to determine these quantities, a τ decay vertex is defined as the point in space that minimizes the sum of the distances to the three daughter pions. In addition, a maximum DOCA of 0.05 mm for all two-body combinations of these sets of tracks is required. Finally, for the h 3 e and h 3 μ modes we also require I > 0.99 for ΔR 2 = 0.05.
The corrected mass [86] is defined as m 2 (πππ) + p 2 T (πππ)+ p T (πππ), where the p T is computed with respect to the τ direction of flight. This quantity, which necessarily requires to know the τ decay vertex and its direction of flight, serves as a good proxy to the real invariant mass of the tau, accounting for the presence of an invisible massless particle, and has been widely used in LHCb analyses involving a neutrino in the final state.
We then combine pairs of daughters to reconstruct candidates, computing (pseudo-)decay vertices of the a candidate for the h 3 h 3 , h 3 μ, and h 3 e categories, in a similar way as for the eμ case. The candidates are reconstructed using two τ leptons for the h 3 h 3 category, or one τ lepton and a selected e(μ) for the h 3 e(h 3 μ) category. The cuts are very similar to those of the eμ candidates. The only exception is the h 3 h 3 category, with tighter cuts, being the maximum distance of flight is 0.25 mm and the IP smaller than 0.1 mm. As for the a daughters, the cuts are again the same except for the h 3 h 3 . For these, the DOCA should not exceed 0.1 mm and both τ leptons are required to be produced promptly, that is, 0.1mm < V r < 5 mm, being V r the radial position of their production vertices.
Signal efficiencies and mass windows are reported in Tables 5 and 6.