We present a novel data-driven method for determining the hadronic interaction strengths of axion-like particles (ALPs) with QCD-scale masses. Using our method, it is possible to calculate the hadronic production and decay rates of ALPs, along with many of the largest ALP decay rates to exclusive final states. To illustrate the impact on QCD-scale ALP phenomenology, we consider the scenario where the ALP-gluon coupling is dominant over the ALP coupling to photons, electroweak bosons, and all fermions for mπ ma 3 GeV. We emphasize, however, that our method can easily be generalized to any set of ALP couplings to SM particles. Finally, using the approach developed here, we provide calculations for the branching fractions of ηc → V V decays, i.e. ηc decays into two vector mesons, which are consistent with the known experimental values.Axion-like particles (ALPs) are hypothetical pseudoscalars whose couplings to the gauge bosons of the Standard Model (SM)-the gluons, photons, and electroweak bosons-are highly suppressed at low energies by a large cut-off scale Λ. ALPs are found in many proposed extensions to the SM (see ), since they naturally address such puzzles as the Strong CP [5-8] and Hierarchy problems [9]. Moreover, ALPs may explain the long-standing anomaly with the magnetic moment of the muon [10], and could provide a portal connecting SM particles to dark matter [11][12][13][14].ALPs are pseudo-Nambu-Goldstone bosons, and therefore, their masses, m a , are expected to be m a Λ. Recently, MeV-to-GeV scale, henceforth QCD-scale, ALPs have received considerable interest [15][16][17][18][19][20][21][22][23][24][25]; however, the phenomenological impact of ALP-gluon interactions is not well understood for QCD-scale ALPs. The effective Lagrangian describing such interactions iswhere c g is the dimensionless agg vertex coupling constant andG µν ≡ 1 2 µναβ G αβ . In this Letter, we present a novel data-driven method for determining the hadronic interaction strengths of QCD-scale ALPs. Using our method, it is possible to calculate the hadronic production and decay rates of ALPs, along with many of the largest ALP decay branching fractions to exclusive final states. To illustrate the impact on QCD-scale ALP phenomenology of c g = 0, we considerfor m π m a 3 GeV; i.e. the scenario where the ALP-gluon coupling is dominant over the ALP coupling to photons (c γ ), electroweak bosons (c EW ), and all * Electronic address: daniel.aloni@weizmann.ac.il † Electronic address: yotam.soreq@cern.ch ‡ Electronic address: mwill@mit.edu fermions (c f ). We emphasize, however, that our method can easily be generalized to any ALP couplings to SM particles. The impact of ALP couplings to photons, electroweak bosons, leptons, and heavy quarks is known [26], while additional direct couplings to light quarks are easily handled within our framework (see the Supplemental Material to this Letter). We begin by noting that ALP-lepton couplings arise at the 3-loop order in this scenario, and therefore, are neglected throughout. ALP couplings t...