# Machine learning assisted droplet trajectories extraction in dense emulsions
and their analysis

Mihir Durve  Center for Life Nano- & Neuro-Science, Fondazione Istituto
Italiano di Tecnologia (IIT), viale Regina Elena 295, 00161 Rome, Italy
Adriano Tiribocchi  Istituto per le Applicazioni del Calcolo del Consiglio
Nazionale delle Ricerche, via dei Taurini 19, 00185, Rome, Italy  Andrea
Montessori  Dipartimento di Ingegneria, UniversitÃ degli Studi Roma tre, via
Vito Volterra 62, Rome, 00146, Italy  Marco Lauricella  Istituto per le
Applicazioni del Calcolo del Consiglio Nazionale delle Ricerche, via dei
Taurini 19, 00185, Rome, Italy  Sauro Succi  Center for Life Nano- & Neuro-
Science, Fondazione Istituto Italiano di Tecnologia (IIT), viale Regina Elena
295, 00161 Rome, Italy  Istituto per le Applicazioni del Calcolo del Consiglio
Nazionale delle Ricerche, via dei Taurini 19, 00185, Rome, Italy  Department
of Physics, Harvard University, 17 Oxford St, Cambridge, MA 02138, United
States

###### Abstract

This work analyzes trajectories obtained by YOLO and DeepSORT algorithms of
dense emulsion systems simulated by Lattice Boltzmann methods. The results
indicate that the individual droplet’s moving direction is influenced more by
the droplets immediately behind it than the droplets in front of it. The
analysis also provide hints on constraints on writing down a dynamical model
of droplets for the dense emulsion in narrow channels.

##  1 Introduction

The last decades have witnessed an impressive rise of machine learning
methods, which have profoundly impacted many areas of science, ranging from
high energy physics and quantum computation to material design, and have
played a vital role in recent milestone discoveries [1, 2]. Many machine
learning-based applications are targeted to perform specific tasks such as
image recognition, natural language processing, sentiments analysis, and
handwriting recognition, to name a few [3, 4, 5, 6], with the aim of improving
accuracy, decreasing processing speed, and reducing human efforts to perform
tasks that are perceived as mundane or labor intensive. In microfluidics, in
particular, they have been used to study shapes and predict transport
properties of flowing droplets [7, 8, 9].

Deep neural networks have become one of the fundamental mathematical models
for implementing machine learning algorithms [10]. Computer vision is one such
application domain that uses deep neural networks to perform essentially two
tasks, object recognition and object tracking. For example, a computer vision
application connected to a camera observing a traffic junction could easily
count the number of cars passing through it and monitor their speed, thus
making it a widespread tool for smart traffic management [11].

Recently, we employed two state-of-the-art computer vision algorithms, namely
You Only Look Once (YOLO) and DeepSORT, to recognize and track droplets in
high internal phase dense emulsions and infer individual droplets
trajectories. These emulsions consist of highly ordered liquid droplets
arranged in crystal-like structures [12, 13, 14], and have shown promising
applications in electrochemical sensing and tissue engineering [15]. Such
materials represent a fundamental challenge to non-equilibrium thermodynamics
as they feature highly non-Newtonian mechanical and rheological properties.
Thus, studying the rich dynamics of these classical many-body systems is
crucial to optimize their disegn as well as the microfluidic devices employed
for their synthesis and applications.

In this work, we revisit the tools employed to infer the trajectories of the
droplets in dense flowing emulsions simulated using lattice Boltzmann methods,
presenting an analysis of such trajectories (obtained via computer vision
tools, see Fig.10 and video4.avi in Ref.[16]) in terms of quantities generally
adopted to characterize the behavior of active matter systems.

The paper is structured as follows. In the next Section we describe the system
under investigation; later in Section 3, we take a brief overview of the
computer vision algorithms deployed to get individual droplet trajectories and
in Section 4 we present the analysis of the inferred trajectories.

##  2 Physical system

The system under study is a soft granular material made of approximately
monodisperse fluid droplets (white region in Fig. 1) immersed within an inter-
droplet continuous phase (dark orange lines) and surrounded by an external
bulk phase (region outside the emulsion). Its structure closely resembles that
of a double emulsion with multi-core morphology [17, 15, 12], where a high
volume fraction of dispersed droplets, generally above the close packing limit
for spheres, arrange in a tightly packed configuration. This material is
produced within a microfluidic channel which, in our system, consists of an
inlet reservoir followed by a thinner channel connected to a further
downstream reservoir. Further details about the simulations, performed via
lattice Boltzmann methods [18, 19, 20], can be found in Ref. [21].

Figure 1:  Snapshot of a dense emulsion simulated by Lattice Boltzmann
methods.

In the next section, we briefly describe the algorithms employed to analyze
videos of the physical system.

##  3 Algorithms

Two algorithms were combined to achieve droplet tracking. The first one,
called You only look once (YOLO), is tasked with identifying the droplets in
an input image while the second one, called DeepSORT, for tracking the
droplets in sequential frames.

###  3.1 You Only Look Once (YOLO)

YOLO is a single-stage state-of-the-art object detection algorithm. The
current version of YOLO (YOLOv5) is the fastest and most accurate object
detector on two commonly used, general-purpose object detection datasets
called Pascal VOC (visual object classes) [22] and Microsoft COCO object
detection datasets [23]. The image analysis speed, i.e. inference speed of the
YOLOv5 networks, is at or above 60 frames per second (FPS) [24, 25] for
general object detection.

The YOLO algorithm is the fastest due to its smart operating procedure [26,
27]. The input is divided in a $S\times S$ grid, with each cell responsible
for detecting an object within the cell. Each grid cell then predicts $B$
bounding boxes with their confidence score for each detected object and $C$
conditional class probabilities for the given object belonging to a specific
class. This information is then combined to produce the final output as a
single bounding box around the detected object and the class of that object.
This final output is then passed to the object tracking DeepSORT algorithm.

The training dataset is used to train the YOLO networks for recognizing the
droplets. The training dataset contains several images of droplets and
associated label files, which include the location and dimensions of the
droplets in each image. The training data is used to get predictions and
update the network parameters based on true and predicted output in an
iterative process. The technical details of the training procedure to train a
YOLO network for droplet recognition is described here [16, 21].

###  3.2 DeepSORT

The DeepSORT algorithm constructs trajectories of all the detected objects by
analyzing sequential frames [28], employing a classical Simple Online Real-
Time tracking module [29] in the first stage. This module uses the Hungarian
algorithm [30] to distinguish detected objects in two consecutive frames and
assigns individual objects their unique identity. The module also uses a
Kalman filtering [31] for predicting the future position of the objects based
on their current positions. At the second level, the deep network learns
object descriptor features to minimize the identity switches as the object
moves in subsequent frames. The YOLO and the DeepSORT algorithms together
accomplish droplet recognition and tracking, thus allowing to obtain the
trajectories of individual droplets (see Fig.10 in Ref. [16]). These ones will
be analyzed in the next section.

##  4 Analysis of the droplet trajectories

Preliminary observations (see video4.avi in Ref. [16]) indicate that, in a
densely packed configuration, the droplets are moving along the same
direction, in a way akin to birds in highly polar flocks [32, 33]. To dig more
deeply into such analogy, we compute a polar order parameter $\psi$ which is
typically used to measure direction consensus in moving particles [34]. It is
computed as

$\psi=\frac{1}{N}\bigg{\lvert}\sum_{i=1}^{N}\frac{{\bf v}_{i}}{|{\bf
v}_{i}|}\bigg{\rvert},$ (1)

were N is the total number of droplets, ${\bf v}_{i}$ is the velocity vector
of the $i^{\text{th}}$ droplet and $|.|$ is the modulus of a vector. If
$\psi=1$ all the droplets move in the same direction while $\psi\approx 0$
means that droplets move randomly.

The polar order parameter is shown in Fig. 2(a). The high average value of
$\psi$ ($\psi>0.94$) indicates that the droplets are moving more or less in
the same direction except at a time interval when they are about to enter the
narrow channel. In Fig. 2(b), we compare an individual droplet’s heading
direction with the average moving direction of the emulsion. This comparison
is made by measuring the cosine of an angle $\theta$ between the velocity
vector of the emulsion (taken as the resultant velocity vector of all the
droplets) and the individual droplet’s velocity vector. Except before entering
the narrow channel, the droplets move in the direction of their bulk motion as
indicated by the cosine value close to unity. In a system containing self-
propelled particles (flock of birds, for example), such a high value of polar
order parameter and cosine of $\theta$ corresponds to a highly ordered state
in which the agents move in a common direction [32]. These observations show
apparent similarities between a flock of birds and the emulsion in terms of
dynamical properties. These similarities raise an interesting question, i.e.
whether the well-known models [34, 35, 36] to mimic the behavior of active
matter systems by taking into account only local interactions can account for
the dynamical properties of the emulsion observed here.

Figure 2: (a) Polar order parameter $\psi$ computed at each frame. It accounts
for moving direction consensus. The black line indicates average of $\psi$
over all the frames. (b) Cosine of angle between individual droplet’s moving
direction and average direction of all the droplets. Cosine of the angle
computed for individual droplet is shown with a unique color. Figure 3:
Scatter plot of individual droplet’s moving direction and average moving
direction of its neighbors. The left column shows neighbor selection in the
opposite direction of the motion of the droplet while the right column shows
neighbors selected in front of the moving droplet. The middle plot is computed
with a full circle with radius $R$ as the neighborhood of a droplet. The range
of $\phi$ is chosen as (0, 1), meaning the actual value of the angle is
$2\pi\phi$. The points on the red line correspond to the instance when
$\theta_{i}=\theta_{S}$, while points within the two black lines indicate
instances where $\theta_{i}$ and $\theta_{S}$ are within 10 degrees of each
other. Thus, points within the band bounded by two black lines show instances
when individual droplet’s moving direction is approximately same as it’s
neighbors.  Figure 4: The neighborhood of the droplet $i$, placed at the
center of a circle of radius $R$, is the sector of the circle enclosed by an
angle $2\pi\phi$. The “front” neighbors are the droplets whose center of mass
is within the sector shaded with dark gray color while the “back” neighbors
are the ones whose center of mass is within the sector shaded with light gray
color. The vector ${\bf v}_{i}$ is the speed of the droplet $i$.

Figure 3 shows a scatter-plot of an individual droplet’s moving direction with
the average direction of its neighbors. The neighbors of a droplet $i$ are
defined as the droplets that are within distance $R$ from the droplet $i$ (see
Fig. 4) and within the sector of the circle with an angle $2\pi\phi$. The
value of $R$ is chosen such that the circle contains the immediate neighbors
of the droplet $i$, while the values of $\phi$ are varied between $0$ and $1$.
We analyzed two cases, one in which the accounted neighbors are in front of
the moving droplets, and a further one in which the neighbors are behind the
moving droplets. This neighborhood definition allows us to see how these
droplets positioned at two different locations affect the motion of an
individual droplet $i$, an effect that can be inferred by counting the number
of instances in which a given droplet is aligned with the average direction of
its neighbors. In Fig. 5, we plot the fraction $F$ of points located within a
band indicated by the two black lines of each snapshot of Fig. 3. The points
within this band correspond to an instance in which an individual droplet’s
moving direction is within ten degrees of its neighbors, i.e.
$|\theta_{i}-\theta_{S}|<0.175$. Here $\theta_{i}$ is the moving direction of
droplet $i$, and $\theta_{S}$ is the moving direction of neighbors of droplet
$i$. The results show that the droplets are aligned for a comparatively longer
time with the droplets behind it rather than the droplets in front of it, thus
suggesting that the droplets pushing from behind have higher influence on the
individual droplet’s moving direction.

Figure 5: Fraction of instances in which a given droplet is closely aligned
with its neighbors as the neighborhood area is varied by changing the angle
$\phi$.

##  5 Conclusions

In this work, we revisit the deep learning-based algorithms to infer the
droplet trajectories by analyzing the output of lattice Boltzmann simulation
of dense emulsions. We measure various quantities like polar order parameter
and deviation of individual droplet’s moving direction with its neighbors. The
results suggest that the droplets are aligned for a comparatively longer time
with the droplets behind it rather than with the droplets in front of it,
meaning that the droplets pushing from behind have higher chance to affect the
direction of an individual droplet. It would be interesting to perform similar
analysis on data generated by active matter systems in order to further
investigate to which extent the analogy with moving emulsions holds and
provide further hints to write down a dynamical model of such emulsions in
confined systems.

##  6 Acknowledgments

The authors acknowledge funding from the European Research Council under the
European Union’s Horizon 2020 Framework Programme (No. FP/2014-2020) ERC Grant
Agreement No.739964 (COPMAT). We gratefully acknowledge the HPC infrastructure
and the Support Team at Fondazione Istituto Italiano di Tecnologia.

## References

  * [1] Andrew W. Senior, Richard Evans, John Jumper, James Kirkpatrick, Laurent Sifre, Tim Green, Chongli Qin, Augustin Žídek, Alexander W. R. Nelson, Alex Bridgland, Hugo Penedones, Stig Petersen, Karen Simonyan, Steve Crossan, Pushmeet Kohli, David T. Jones, David Silver, Koray Kavukcuoglu, and Demis Hassabis.  Improved protein structure prediction using potentials from deep learning.  Nature, 577(7792):706–710, Jan 2020.
  * [2] Giuseppe Carleo, Ignacio Cirac, Kyle Cranmer, Laurent Daudet, Maria Schuld, Naftali Tishby, Leslie Vogt-Maranto, and Lenka Zdeborová.  Machine learning and the physical sciences.  Rev. Mod. Phys., 91:045002, Dec 2019.
  * [3] Darmatasia and M. I. Fanany.  Handwriting recognition on form document using convolutional neural network and support vector machines (cnn-svm).  2017 5th International Conference on Information and Communication Technology (ICoIC7), pages 1–6, 2017.
  * [4] N. H. Tandel, H. B. Prajapati, and V. K. Dabhi.  Voice recognition and voice comparison using machine learning techniques: A survey.  2020 6th International Conference on Advanced Computing and Communication Systems (ICACCS), pages 459–465, 2020.
  * [5] Savita Ahlawat, Amit Choudhary, Anand Nayyar, Saurabh Singh, and Byungun Yoon.  Improved handwritten digit recognition using convolutional neural networks (cnn).  Sensors, 20(12), 2020.
  * [6] Kun Han, Dong Yu, and Ivan Tashev.  Speech emotion recognition using deep neural network and extreme learning machine.  In Interspeech 2014, September 2014.
  * [7] P. Hadikhani, N. Borhani, S.M.H. Hashemi, and D. Psaltis.  Learning from droplet flows in microfluidic channels using deep neural networks.  Scientific Reports, 9:8114, 2019.
  * [8] Yassine Mahdi and Kamel Daoud.  Microdroplet size prediction in microfluidic systems via artificial neural network modeling for water-in-oil emulsion formulation.  Journal of Dispersion Science and Technology, 38(10):1501–1508, 2017.
  * [9] Jian Wei Khor, Neal Jean, Eric S. Luxenberg, Stefano Ermon, and Sindy K. Y. Tang.  Using machine learning to discover shape descriptors for predicting emulsion stability in a microfluidic channel.  Soft Matter, 15:1361–1372, 2019.
  * [10] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [11] Tousif Osman, Shahreen Shahjahan Psyche, J. M. Shafi Ferdous, and Hasan U. Zaman.  Intelligent traffic management system for cross section of roads using computer vision.  2017 IEEE 7th Annual Computing and Communication Workshop and Conference (CCWC), pages 1–7, 2017.
  * [12] Andrea Montessori, Adriano Tiribocchi, Michał Bogdan, Fabio Bonaccorso, Marco Lauricella, Jan Guzowski, and Sauro Succi.  Translocation dynamics of high-internal phase double emulsions in narrow channels.  Langmuir, 37(30):9026–9033, Aug 2021.
  * [13] Andrea Montessori, Michele La Rocca, Pietro Prestininzi, Adriano Tiribocchi, and Sauro Succi.  Deformation and breakup dynamics of droplets within a tapered channel.  Physics of Fluids, 33(8):082008, 2021.
  * [14] Michał Bogdan, Andrea Montessori, Adriano Tiribocchi, Fabio Bonaccorso, Marco Lauricella, Leon Jurkiewicz, Sauro Succi, and Jan Guzowski.  Stochastic jetting and dripping in confined soft granular flows.  Phys. Rev. Lett., 128:128001, Mar 2022.
  * [15] M. Costantini, C. Colosi, J. Guzowski, A. Barbetta, J. Jaroszewicz, W. Swieszkowski, M. Dentini, and P. Garstecki.  Highly ordered and tunable polyhipes by using microfluidics.  J. Mater. Chem. B, 2:2290–2300, 2014.
  * [16] Durve, Mihir, Bonaccorso, Fabio, Montessori, Andrea, Lauricella, Marco, Tiribocchi, Adriano, and Succi, Sauro.  Tracking droplets in soft granular flows with deep learning techniques.  Eur. Phys. J. Plus, 136(8):864, 2021.
  * [17] A. S. Utada, E. L. Lorenceau, D. R. Link, P. D. Kaplan, H. A. Stone, and D. A. Weitz.  Monodisperse double emulsions generated from a microcapillary device.  Science, 308:537–541, 2005.
  * [18] A Montessori, P Prestininzi, M La Rocca, and S Succi.  Lattice boltzmann approach for complex nonequilibrium flows.  Physical Review E, 92(4):043308, 2015.
  * [19] Christophe Coreixas, Bastien Chopard, and Jonas Latt.  Comprehensive comparison of collision models in the lattice boltzmann framework: Theoretical investigations.  Physical Review E, 100(3):033305, 2019.
  * [20] S. Succi.  The Lattice Boltzmann Equation: For Complex States of Flowing Matter.  Oxford University Press, 2018.
  * [21] Mihir Durve, Fabio Bonaccorso, Andrea Montessori, Marco Lauricella, Adriano Tiribocchi, and Sauro Succi.  A fast and efficient deep learning procedure for tracking droplet motion in dense microfluidic emulsions.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 379(2208):20200400, 2021.
  * [22] The pascal visual object classes homepage.  http://host.robots.ox.ac.uk/pascal/VOC/.
  * [23] Coco dataset homepage.  http://cocodataset.org.
  * [24] Fangbo Zhou, Huailin Zhao, and Zhen Nie.  Safety helmet detection based on yolov5.  In 2021 IEEE International Conference on Power Electronics, Computer Applications (ICPECA), pages 6–11, 2021.
  * [25] Luiz Carlos M. Junior and JosÃ© Alfredo C. Ulson.  Real time weed detection using computer vision and deep learning.  In 2021 14th IEEE International Conference on Industry Applications (INDUSCON), pages 1131–1137, 2021.
  * [26] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi.  You only look once: Unified, real-time object detection.  2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 779–788, 2016.
  * [27] Joseph Redmon and Ali Farhadi.  Yolov3: An incremental improvement.  ArXiv:1804.02767v1, 2018.
  * [28] N. Wojke, A. Bewley, and D. Paulus.  Simple online and realtime tracking with a deep association metric.  2017 IEEE International Conference on Image Processing (ICIP), pages 3645–3649, 2017.
  * [29] Alex Bewley, Zongyuan Ge, Lionel Ott, Fabio Ramos, and Ben Upcroft.  Simple online and realtime tracking.  In 2016 IEEE International Conference on Image Processing (ICIP), pages 3464–3468, 2016.
  * [30] H. W. Kuhn.  The hungarian method for the assignment problem.  Naval Research Logistics Quarterly, 2(1-2):83–97, 1955.
  * [31] R. E. Kalman.  A New Approach to Linear Filtering and Prediction Problems.  Journal of Basic Engineering, 82(1):35–45, 03 1960.
  * [32] Andrea Cavagna, Lorenzo Del Castello, Irene Giardina, Tomas Grigera, Asja Jelic, Stefania Melillo, Thierry Mora, Leonardo Parisi, Edmondo Silvestri, Massimiliano Viale, and Aleksandra M. Walczak.  Flocking and turning: a new model for self-organized collective motion.  Journal of Statistical Physics, 158(3):601–627, Feb 2015.
  * [33] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic.  Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study.  Proceedings of the National Academy of Sciences, 105(4):1232–1237, 2008.
  * [34] Tamás Vicsek, András Czirók, Eshel Ben-Jacob, Inon Cohen, and Ofer Shochet.  Novel type of phase transition in a system of self-driven particles.  Phys. Rev. Lett., 75:1226–1229, Aug 1995.
  * [35] IAIN D. COUZIN, JENS KRAUSE, RICHARD JAMES, GRAEME D. RUXTON, and NIGEL R. FRANKS.  Collective memory and spatial sorting in animal groups.  Journal of Theoretical Biology, 218(1):1–11, 2002.
  * [36] Lucas Barberis and Fernando Peruani.  Large-scale patterns in a minimal cognitive flocking model: Incidental leaders, nematic patterns, and aggregates.  Phys. Rev. Lett., 117:248001, Dec 2016.


TTP23-042, P3H-23-070

# Phenomenology of axion-like particles with universal fermion couplings –
revisited

Giovani Dalla Valle Garcia<EMAIL_ADDRESS>Institut für
Astroteilchen Physik, Karlsruher Institut für Technologie (KIT), Hermann-von-
Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany  Felix Kahlhoefer
<EMAIL_ADDRESS>Institute for Theoretical Particle Physics (TTP), Karlsruhe
Institute of Technology (KIT), D-76131 Karlsruhe, Germany  Maksym Ovchynnikov
<EMAIL_ADDRESS>Institut für Astroteilchen Physik, Karlsruher
Institut für Technologie (KIT), Hermann-von-Helmholtz-Platz 1, 76344
Eggenstein-Leopoldshafen, Germany  Instituut-Lorentz, Leiden University, Niels
Bohrweg 2, 2333 CA Leiden, The Netherlands  Andrii Zaporozhchenko
<EMAIL_ADDRESS>Taras Shevchenko National University of Kyiv, 64
Volodymyrs’ka str., Kyiv 01601, Ukraine

###### Abstract

Axion-like particles (ALPs) emerge in many extensions of the Standard Model as
pseudo-Goldstone bosons of a spontaneously broken global symmetry.
Understanding their phenomenology in high-energy collisions is crucial for
optimizing experimental searches and understanding the exploration potential
of future experiments. In this paper, we revise the phenomenology of ALPs with
universal couplings to fermions. In particular, we analyze the hierarchy and
uncertainty of the various ALP production channels depending on the proton
collision energy and the placement of the experiment, and provide improved
calculations of the hadronic decay modes.

##  I Introduction

Axion-like particles (ALPs) $a$ are pseudoscalar particles that arise in
theories with spontaneously broken global chiral symmetries, generalizing the
idea of the QCD axion – a hypothetical light particle capable of solving the
strong CP problem Peccei and Quinn (1977); Weinberg (1978); Wilczek (1978).
While the QCD axion obtains its mass directly from its coupling to gluons, a
generic ALP may interact with various particles and have an arbitrary mass
Jaeckel and Ringwald (2010). The lowest-order gauge-invariant Lagrangian of
the ALP interaction takes the form Georgi et al. (1986); Bauer et al. (2017);
Brivio et al. (2017)

$\displaystyle\mathcal{L}=\frac{a}{f}\bigg{(}$ $\displaystyle
c_{G}\frac{\alpha_{s}}{4\pi}G_{\mu\nu}^{c}\tilde{G}^{\mu\nu,c}+c_{W}\frac{\alpha_{W}}{4\pi}W^{\mu\nu,c}\tilde{W}^{c}_{\mu\nu}$
$\displaystyle+c_{B}\frac{\alpha_{B}}{4\pi}B_{\mu\nu}\tilde{B}^{\mu\nu}\bigg{)}+\frac{\partial^{\mu}a}{\Lambda}\sum_{F}\bar{\Psi}_{F}\,C_{F}\,\gamma_{\mu}\,\Psi_{F}\;,$
(1)

where $G,W,B$ are field strengths of the Standard Model $SU_{c}(3)$,
$SU_{L}(2)$, and $U_{Y}(1)$ gauge groups,
$\alpha_{s},\alpha_{W},\alpha_{B}=g^{2}/4\pi$ are corresponding running
couplings,
$\tilde{G}_{\mu\nu}=\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}G^{\alpha\beta}$ is
the dual strength and $\Psi_{F}$ denote the SM fermion multiplets.
Furthermore, $f$ is a dimensional scale, $c_{i}$ are dimensionless parameters
and $C_{F}$ are (hermitian) matrices characterizing. For light ALPs with
$m_{a}\simeq 1\text{ GeV}$, past experiments have excluded combinations
$f/c_{i}\lesssim 1\text{ TeV}$ Beacham et al. (2020); Antel et al. (2023).
Therefore, ALPs belong to the class of the so-called Feebly-Interacting
Particles, or FIPs. Many studies have explored the phenomenology of the
individual terms in the Lagrangian above, in particular couplings to gluons
Aloni et al. (2019a, b); Chakraborty et al. (2021), photons Döbrich et al.
(2016); Dolan et al. (2017); Döbrich et al. (2019a), electroweak gauge bosons
Izaguirre et al. (2017); Alonso-Álvarez et al. (2019); Gavela et al. (2019),
fermions Dolan et al. (2015); Döbrich et al. (2019b); Carmona et al. (2021),
the effect of flavour-violation Cornella et al. (2020); Calibbi et al. (2021)
and renormalisation group evolution Chala et al. (2021); Bauer et al. (2021,
2022); Ferber et al. (2023), as well as the interplay between different
couplings Ertas and Kahlhoefer (2020); Jerhot et al. (2022); Liu et al.
(2023); Bruggisser et al. (2023).

One case of particular interest are ALPs that interact dominantly with
fermions with universal and flavour-diagonal couplings,

$\mathcal{L}_{\text{eff}}=\frac{\partial_{\mu}a}{f}\left(c_{\ell}\sum_{\ell}\bar{\ell}\gamma^{\mu}\gamma_{5}\ell+c_{q}\sum_{q}\bar{q}\gamma^{\mu}\gamma_{5}q\right),$
(2)

where $\ell$ are leptons and $q$ are quarks. The case $c_{\ell}=c_{q}$ has
been identified by the Physics Beyond Colliders (PBC) initiative Beacham et
al. (2020) as one of the benchmark models (called BC10) to demonstrate the FIP
exploration potential of future experiments. Due to a lack of in-depth
theoretical studies of this model, the description of the phenomenology of
such ALPs proposed in Ref. Beacham et al. (2020) suffers from two issues.
First, following Ref. Dolan et al. (2015), the contribution of the hadronic
decays to the total ALP decay width is assumed to be negligible. While this is
a reasonable assumption for light ALPs with $m_{a}\ll 1\text{ GeV}$, where the
only relevant decay mode is into three pions Bauer et al. (2021), the hadronic
decays may actually dominate the total width for heavier ALPs, as indicated by
calculations of the ALP decay width into quark and gluon pairs Bauer et al.
(2021); Ferber et al. (2023).

Another problem is that the production of such ALPs at beam dump experiments
has been approximated by considering decays $B\to X_{s}+a$, where
$X_{s}=K,K^{*}$ only Döbrich et al. (2019b); Ferber et al. (2023).
Nevertheless, there may be a sizable contribution from the $B$ decays into
other resonances, $X_{s}=K_{1},K^{*}_{2},K^{*}_{0}$. For the case of light
Higgs-like scalars Boiarska et al. (2019), which couple to the $b\to s$
operator in a way similar to ALPs, these decays have been shown to correspond
to $1/3$ of all possible decays. The same effect can be expected to be
relevant for the case of ALPs with interactions as in eqs. (1) and (2).
Second, additional production processes arise due to the mixing of the ALPs
with light pseudoscalar mesons $m^{0}=\pi^{0}/\eta/\eta^{\prime}$, similar to
the ALPs coupled to gluons.

In this paper, we address these issues by revising the phenomenology of the
model given in eq. (2). Our goal is to provide a detailed and comprehensive
description of the PBC BC10, which the community may easily implement to
consistently interpret ALP constraints from existing experiments and study the
projected sensitivities of proposed searches. The results summarized in this
paper are also accessible in a Mathematica notebook supplementing the paper
(see Appendix B).111Available on https://github.com/maksymovchynnikov/ALPs-
phenomenology We implement the revised phenomenology in SensCalc Ovchynnikov
et al. (2023a) – a public code that uses the semi-analytic approach to
calculate the number of events with decays of FIPs and sensitivities of
proposed experiments.

The remainder of this work is organized as follows. In Sec. II, we discuss the
details of the ALP model that we consider, and in particular, the choice of
the scale at which the underlying Lagrangian is defined. In Sec. III, we
discuss various contributions to the ALP production flux depending on the
proton collision energy for the given experiment and its geometric placement.
In Sec. IV, we study the decay palette of the ALP. We conclude in Sec. V.

##  II Model details

The phenomenology of GeV-scale ALPs depends on the scale $\Lambda$ at which
the interactions in eq. (2) are defined. This is because of a non-trivial
renormalization group (RG) flow, which generates additional effective
couplings absent in eq. (2) and modifies the initial couplings.

The RG evolution has been thoroughly studied for general models of ALPs (see
Refs. Bauer et al. (2022, 2021); Ferber et al. (2023) and references therein).
The evolution is split onto the flow from $\Lambda$ down to the electroweak
scale $\mu_{w}\equiv m_{t}$ and from $\mu_{w}$ down to the scale of the
process with ALPs, which is of the order of the ALP mass. For the processes
with hadrons, there is one additional scale $\sim 4\pi\Lambda_{\text{QCD}}$
where the perturbative QCD should be matched with Chiral Perturbation Theory
(ChPT).

The couplings $c_{G},c_{W},c_{B}$ entering the Lagrangian in eq. (1) are
scale-invariant at least up to the second order in the loop expansion.
Therefore, they will only be generated from the initial Lagrangian in eq. (2)
by integrating out heavy quarks. This is not the case for the fermion
couplings from eq. (2), which may evolve due to electroweak and strong
interaction loops. To study the dynamics of these couplings, we solve the RG
equations from Ref. Bauer et al. (2021). The behavior of the couplings
$c_{i},i=u,d,s,\ell$ to various SM fermions at the scale $\mu=2\text{ GeV}$ as
a function of $\Lambda$, assuming the universal coupling $c_{i}(\Lambda)=1$,
is shown in Fig. 1 (left panel).

Figure 1: The RG flow of the quark and lepton couplings from the Lagrangian
(2) (the left panel) and modulus of the flavor-changing $abs$ coupling from
eq. (8) at $f=1\text{ TeV}$ (the right panel), assuming the initial condition
$c_{q}(\Lambda)=c_{l}(\Lambda)=1$. The curves show the values of the couplings
at the EW scale $\mu=2\text{ GeV}$. The vertical dashed line shows the
reference scale $\Lambda=1\text{ TeV}$ used as our benchmark choice.

From the figure we see that all the lepton couplings $c_{l}$ evaluated at the
EW scale differ from $c_{l}(\Lambda)$ by less than $20\%$ for
$\Lambda\lesssim\mathcal{O}(10\text{ TeV})$; their dynamics is flavor-
independent. The situation is different for the quark couplings. The evolution
from the scale $\Lambda$ and down to $\mu_{w}$ is flavor-dependent because
different quarks have different weak isospins. The evolution from $\mu_{w}$
and down to the scale of interest is flavor-universal up to tiny EM
corrections, which may be neglected. As a result, the RG flow violates quark
coupling universality. For instance, assuming $\Lambda=1\text{ TeV}$, the
relative difference is $(c_{d}-c_{u})/c_{d}\simeq 10\%$. As we will see in
Sec. III, it cannot be neglected when studying the interactions of the ALPs
with neutral pions (see Sec. III).

The ALP-gluon interaction is another type of interaction important for ALP
production and decay. At the leading order in $\alpha_{s}$, the matrix element
of the type $GG\to a$ (the gluon fusion process) or $a\to GG$ (the ALP decay
into a pair of gluons) is generated by the following matrix element:

$\mathcal{M}_{GG\leftrightarrow a}\approx
c^{\text{eff}}_{G}(m_{a})\frac{\alpha_{s}}{2\pi
f}G^{a}_{\mu\nu}(p_{G,1})\tilde{G}^{\mu\nu,a}(p_{G,2}),$ (3)

where $G^{a}_{\mu\nu}(p)$ is the linear part of the gluon strength tensor with
the replacements $\partial_{\mu}\to ip_{\mu}$ and
$G_{\mu}^{a}\to\epsilon_{\mu}^{a}(p)$, with $\epsilon_{\mu}^{a}$ being the
polarization vector. The effective coupling $c^{\text{eff}}_{G}(m_{a})$ is
Bauer et al. (2021)

$c^{\text{eff}}_{G}\equiv\sum_{q}c_{q}B_{1}\left(\frac{4m_{q}^{2}}{m_{a}^{2}}\right)$
(4)

with

$B_{1}=1-\tau\times\begin{cases}\arcsin^{2}\left(\frac{1}{\sqrt{\tau}}\right),\quad\tau\geq
1\\\
\left(\frac{\pi}{2}+\frac{i}{2}\ln\left[\frac{1+\sqrt{1-\tau}}{1-\sqrt{1-\tau}}\right]\right)^{2},\quad\tau<1\end{cases}\,.$
(5)

and $\tau_{f}=4m_{f}^{2}/m_{a}^{2}$. In the regime $\tau\gg 1$, the function
$B_{1}$ behaves as $B_{1}(\tau\gg 1)\approx-(3\tau)^{-1}$. In the opposite
regime, $B_{1}(\tau\ll 1)\approx 1$. Therefore, the coupling to gluons is
mainly generated by the quarks lighter than the ALP.

Additionally, the RG flow (via loops involving top quarks and charged weak
gauge bosons) generates the flavor-changing neutral current (FCNC) coupling
$D_{i}D_{j}a$, where $D=d,s,b$ are down quarks and $i\neq j$ Bauer et al.
(2022):

$\mathcal{L}_{\text{FCNC}}=\frac{\partial_{\mu}a}{f}\sum_{i\neq
j}c_{ij}\,\bar{q}_{i}\gamma^{\mu}P_{L}q_{j}$ (6)

where $P_{L}=\frac{1}{2}(1-\gamma_{5})$ and $c_{ij}$ is the model-dependent
coupling:

$\displaystyle c_{ij}=-V_{ti}^{\ast}V_{tj}\Big{\\{}$
$\displaystyle-\frac{1}{6}I_{t}(\Lambda,\mu_{w})$
$\displaystyle+\dfrac{y_{t}^{2}}{16\pi^{2}}$
$\displaystyle\bigg{[}c_{q}\bigg{(}\ln{\dfrac{\mu_{w}^{2}}{m_{t}^{2}}}+\frac{1}{2}+3\dfrac{1-x_{t}+\ln{x_{t}}}{(1-x_{t})^{2}}\bigg{)}$
$\displaystyle+\dfrac{9\alpha_{\text{EM}}}{4\pi
s_{w}^{2}}\left(3c_{q}+c_{l}\right)\dfrac{1-x_{t}+x_{t}\ln{x_{t}}}{(1-x_{t})^{2}}\bigg{]}\Big{\\}}\,,$
(7)

with $x_{t}=m_{t}^{2}/m_{W}^{2}$. The term $I_{t}(\Lambda,\mu_{w})$ represents
the contribution of the RG flow from the scale $\Lambda$ down to
$\mu_{w}=m_{t}$ and vanishes if $\Lambda=\mu_{w}$.

Neglecting the mass of the lighter quark in eq. (6), the FCNC Lagrangian may
be rewritten as Batell et al. (2011)

$\mathcal{L}_{\text{FCNC}}=ia\sum_{i,j}C_{ij}\bar{q}_{i}(1+\gamma_{5})q_{j}+\text{h.c.},$
(8)

where $C_{ij}\equiv c_{ij}m_{q_{j}}/2f$, where $m_{q_{j}}$ is the mass of the
heavier quark among the pair $q_{i}q_{j}$.

The FCNC couplings for the transitions $s\to d$ and $b\to s$ have a huge
impact on the ALP phenomenology as they may dominate the production of the
ALPs depending on the amounts of kaons and $B$ mesons produced in the given
experiment. The value of this coupling is very sensitive to $\Lambda$, growing
by two orders of magnitude if increasing the scale from $\Lambda=\mu_{w}$ to
$\Lambda\sim 10\text{ TeV}$ (see Fig. 1, right panel).222As a cross-check of
the implementation, we reproduce the values of $c_{ij}$ for the scale
$\Lambda=4\pi\text{ TeV}$ reported in Bauer et al. (2022).

Taking this into account, we consider two representative choices of $\Lambda$:
the one with $\Lambda\equiv\mu_{w}$, and another one with $\Lambda\equiv
1\text{ TeV}$, which is the reference scale used for the PBC BC10 benchmark
Beacham et al. (2020).

##  III ALP production

The ALPs in eq. (2) may be produced by decays of kaons and $B$ mesons, via the
mixing with light pseudoscalar mesons $m^{0}$, or by deep inelastic
scatterings (DIS).

We describe these production channels in detail in the section below. The
amounts of the produced mesons and the DIS cross-section are collision energy
dependent. Therefore, to make an experiment-independent comparison, we will
consider three collision energies $\sqrt{s}\approx 16,28\text{ GeV}$, and $13$
TeV (corresponding to the collision energies at DUNE, the SPS beam and the
LHC). We took the meson production fractions from the SensCalc repository
Ovchynnikov et al. (2023a).

###  III.1 Decays of $B,K$ mesons

Scale $\Lambda$ | $|C_{bs}|$ | $|C_{bd}|$ | $|C_{sd}|$
---|---|---|---
$m_{t}$ | $2.9\cdot 10^{-5}$ | $5.6\cdot 10^{-6}$ | $4.1\cdot 10^{-9}$
$1\text{ TeV}$ | $1.8\cdot 10^{-3}$ | $3.4\cdot 10^{-4}$ | $2.8\cdot 10^{-7}$
Table 1: The values of the FCNC couplings from the Lagrangian (8) with
$f=1\text{ GeV}$, assuming the model (2) with $c_{q}=c_{\ell}=1$ and two
scales at which it is defined: $\Lambda=m_{t}$, and $\Lambda=1\text{ TeV}$.

We will consider the interactions $abs$, $abd$, and $asd$, for which the quark
running in the loop is the top quark; the other interactions are heavily
suppressed by the Yukawas of lighter quarks and/or CKM elements and are
irrelevant. The corresponding decay processes are $B\to a+X_{s}$, $B\to
a+\pi$, and $K\to a+\pi$. As we see from eq. (7), the values of the couplings
describing these transitions differ only by the CKM products
$V_{ti}^{*}V_{tj}$. This product is the largest for the $b\to s$ transition;
however, the other processes are also important. Namely, the $abd$ coupling is
suppressed by $|V_{t\to d}/V_{t\to s}|\approx 0.2$; however, the process $B\to
a+\pi$ is the only possible above the threshold $B\to K+a$. The relative
suppression of the $sd$ coupling is even larger, $|V_{t\to d}|\approx 8\cdot
10^{-3}$. However, depending on the experiment, the number of kaons may be
much larger than that of $B$ mesons, which may compensate for this
suppression. The values of the corresponding couplings for the two different
scales $\Lambda=1\text{ TeV}$ or $\Lambda=\mu_{w}$ are given in Table 1. In
particular, the value $C_{bs}(1\text{ TeV})$ matches with the one used for the
BC10 model Beacham et al. (2020); Batell et al. (2011).

Having the operator of the FCNC interaction (8), one may calculate the matrix
elements of the processes $B/K\to a+X$, where $X$ is a hadronic state
containing an $s$ quark or a $d$ quark. They have the form

$\mathcal{M}_{m\to a+m^{\prime}}=iC_{QQ^{\prime}}\left(\mathcal{M}^{(S)}_{m\to
m^{\prime}}+\mathcal{M}^{(P)}_{m\to m^{\prime}}\right),$ (9)

where the parity-even and parity-odd transition matrix elements are

$\displaystyle\mathcal{M}_{m\to m^{\prime}}^{(S)}$ $\displaystyle\equiv\langle
m^{\prime}|\bar{Q}^{\prime}Q|m\rangle,$ $\displaystyle\mathcal{M}^{(P)}_{m\to
m^{\prime}}$ $\displaystyle\equiv\langle
m^{\prime}|\bar{Q}^{\prime}\gamma_{5}Q|m\rangle$ (10)

Because of the parity conservation in QCD, if $m,m^{\prime}$ have the same
parity, only $\mathcal{M}_{m\to m^{\prime}}^{(S)}$ contributes, while for
$m^{\prime}$ having a different parity than $m$ only the
$\mathcal{M}^{(P)}_{m\to m^{\prime}}$ is non-zero.

The matrix elements (10) match with the matrix elements $M_{XX^{\prime}}$ from
eq. (B.7) from Boiarska et al. (2019), used to compute the production of the
Higgs-like scalars, which is caused by the similarity of the FCNC operator for
ALPs and Higgs-like scalars Grzadkowski and Krawczyk (1983); Leutwyler and
Shifman (1990); Haber et al. (1987); Chivukula and Manohar (1988). Therefore,
instead of computing the branching ratios using eq. (9) one may use the
results of Ref. Boiarska et al. (2019) after the rescaling of the branching
ratio with the proper coupling.

Ref. Boiarska et al. (2019) used the matrix elements computed using light-cone
QCD sum rules and considered the mesons $X_{s}=K,K^{*}(892)$, $K^{*}(1410)$,
$K^{*}(1680)$, $K^{0*}(700)$, $K^{0*}(1430)$, $K_{1}(1270)$,
$K_{1}(1430),K_{2}^{*}$, and $X_{d}=\pi$.

Figure 2: Branching ratios of the 2-body decays $B^{+}\rightarrow X+a$, where
$X$ is a hadron that contains an s or d quark. By the $K^{*}_{0}$ channel, we
denote the final states $K^{*}_{0}(700)$ and $K^{*}_{0}(1430)$, by $K_{1}$ —
$K_{1}(1270),K_{1}(1400)$. The dashed black lines correspond to the
probability of the process $B\to K/K^{*}(892)+a$ considered previously in the
literature.

The branching ratios of various decays are shown in Fig. 2. Compared to the
literature where only the decays $B\to a+K/K^{*}(892)$ have been considered
Aloni et al. (2019a); Beacham et al. (2020); Jerhot et al. (2022), we find
almost 4 times larger total production probability. In particular, in the
domain of masses $m_{a}\lesssim 1\text{ GeV}$, the dominant production channel
is into $K_{1}$ and $K^{0*}$ mesons.

###  III.2 Mixing with neutral mesons

####  III.2.1 Interaction

If the ALP is light ($m_{a}\lesssim 2\text{ GeV}$), the description of its
hadronic interactions in terms of the $qq$ and $GG$ operators from eq. (1)
becomes inadequate since the QCD enters the non-perturbative regime. Instead,
light mesons and their interactions represent the strongly interacting sector.

We follow the existing studies Bauer et al. (2021, 2022); Aloni et al. (2019a)
and obtain the Lagrangian of the ALP interactions with the pseudoscalar mesons
$P=\pi,\eta,K,\dots$ by using the matching of the operator (2) with the ChPT
Lagrangian.

Details are provided in Appendix A; we summarize the main features below. In
general, the interaction (2) leads to the kinetic mixings of the ALP with
neutral pseudoscalar mesons $m^{0}$. This contrasts with the case of the ALPs
coupled to gluons, where the gluon operator also induces the mass mixings
Aloni et al. (2019a). We need to diagonalize the kinetic term to find the
relevant interactions. The fields $m^{0}$ entering the Lagrangian are related
to the mass eigenstates $m^{0}_{\text{phys}},a_{\text{phys}}$ by

$m^{0}\approx
m^{0}_{\text{phys}}+\theta_{m^{0}a}a_{\text{phys}}+\sum_{m^{0^{\prime}}\neq
m^{0}}\theta_{m^{0}m^{0^{\prime}}}m^{0^{\prime}}_{\text{phys}}$ (11)

Here, the second term is due to the kinetic mixing with the ALPs, and the
third one appears from the mass mixing between the mesons emerging from the
minimal ChPT breaking term.

In the limit $m_{s}\gg m_{u,d}$, the mixing angles are

$\displaystyle\theta_{\pi^{0}a}$
$\displaystyle\approx\frac{f_{\pi}F(m_{a})}{f}\frac{m_{a}^{2}}{(m_{a}^{2}-m_{\pi^{0}}^{2})}\left((c_{d}-c_{u})+\delta\frac{m_{\pi^{0}}^{2}}{3}\left[\frac{m_{a}^{2}(c_{d}+2c_{s}+c_{u})}{m_{a}^{2}-m_{\eta^{\prime}}^{2}}+\frac{2m_{a}^{2}(-c_{s}+c_{u}+c_{d})}{m_{a}^{2}-m_{\eta}^{2}}\right]\right),$
$\displaystyle\theta_{\eta a}$
$\displaystyle\approx\frac{f_{\pi}F(m_{a})}{f}\sqrt{\frac{2}{3}}\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta}^{2}}\left((c_{u}+c_{d}-c_{s})-\delta\frac{m_{\pi^{0}}^{2}(c_{u}-c_{d})}{m_{a}^{2}-m_{\pi^{0}}^{2}}\right),$
(12) $\displaystyle\theta_{\eta^{\prime}a}$
$\displaystyle\approx\frac{f_{\pi}F(m_{a})}{f}\frac{1}{\sqrt{3}}\frac{m_{a}^{2}}{m_{a}^{2}-m_{\eta^{\prime}}^{2}}\left(-(c_{d}+2c_{s}+c_{u})-\delta\frac{m_{\pi^{0}}^{2}(c_{u}-c_{d})}{m_{\pi^{0}}^{2}-m_{a}^{2}}\right),$

where $\delta=(m_{d}-m_{u})/(m_{d}+m_{u})$ is the isospin symmetry breaking
parameter, $f_{\pi}\approx 93\text{ MeV}$ is the pion decay constant, and
$F(m_{a})$ is a phenomenological function ensuring the drop of the VMD
contribution according to quark counting sum rules Aloni et al. (2019a);
Jerhot et al. (2022). Similar to Aloni et al. (2019a); Cheng et al. (2022), in
eq. (12), we fix $\theta_{\eta}=\arcsin(-1/3)$, motivated by the fact that
various phenomenological studies of the effective Lagrangian of the decays of
light mesons consider this value. The expressions (12) are given in terms of
the couplings $c_{u},c_{d},c_{s}$ instead of a single $c_{q}$ to account for
the RG flow (remind Sec. II).

In the first order on the parameter $f_{\pi}/f\ll 1$ and far from resonance
domains $m_{a}=m_{m^{0}}$, the ALP and meson masses are left unchanged.
Therefore, the only impact of the diagonalization (11) is the appearance of
new interactions between the ALPs and mesons.

The ALP mixing with $\pi^{0}$ emerges either from the RG flow (the first term
in eq. (12)), or via the mixing of unphysical $\pi^{0}$ with $\eta$ and
$\eta^{\prime}$ (the second term).

The behavior of the squared mixing angles as a function of the ALP mass for
the two reference scales $\Lambda=m_{t}$ and $\Lambda=1\text{ TeV}$ is shown
in Fig. 3.

Figure 3: The ALP mass dependence of the square of the mixing angle of the ALP
with neutral mesons $\pi^{0}/\eta/\eta^{\prime}$ (eq. (12)). The solid lines
show the results assuming the scale $\Lambda=1\text{ TeV}$, while the dashed
line corresponds to $\Lambda=m_{t}$. The vertical gray bands correspond to the
vicinity of $m_{a}=m_{m^{0}}$ where the approximate description of the ALP
interactions via the small mixing with the neutral mesons breaks down (see
text for details).

####  III.2.2 Production

All relevant production processes with ALPs and light neutral mesons occur
through their mixings (12). We assume that the ALP production cross-section
due to the mixing is given by

$\sigma_{\text{prod,mixing}}=\sum_{m^{0}=\pi^{0}/\eta/\eta^{\prime}}\sigma_{\text{prod},m^{0}}\times|\theta_{m^{0}a}|^{2}$
(13)

However, depending on the ALP mass, its kinematics would be different from the
corresponding meson kinematics. We follow the procedure described in Jerhot et
al. (2022) to account for this.

We should stress that this description does not account for the ALP mass
dependence of the production cross-section. Clearly, the production of the
ALPs heavier than $m^{0}$ at the unit mixing angle must be kinematically
suppressed, which is not considered in eq. (13). This point is to be improved
in future works.

###  III.3 Deep inelastic scattering process

Another important process of ALP production is deep inelastic scattering. At
the parton level and at the leading order in $\alpha_{s}$, it is described by
the fusion333It is important to note that the higher-order processes which we
do not consider in this paper may contribute to the flux of the ALPs flying
off-axis.

$q+\bar{q}\to a,\quad G+G\to a,$ (14)

where for the second process, the matrix element is given by eq. (3). The
parton model applicability breaks down if the characteristic scale of the
process $\sqrt{s_{qq}}=m_{a}$ becomes comparable with $\Lambda_{\text{QCD}}$.
We “turn on” the process (14) at $m_{a}=1.5\text{ GeV}$.

The DIS process is the only relevant production channel for heavy ALPs with
$m_{a}>m_{B}-m_{K}$, given that its kinematic threshold is extended until the
center-of-mass energy at the experiment.

To estimate the cross-section of the process (14), we implement the fermionic
and gluonic matrix elements in MadGraph5 Alwall et al. (2014) using FeynRules
Alloul et al. (2014); Christensen and Duhr (2009) and generate the leading-
order processes of the gluon and quark fusion. For the parton distribution
function, we use NNPDF 3.1 NNLO set, which is a common choice for FIP
sensitivity studies Berlin and Kling (2019).

We have found that the quark fusion is strongly suppressed compared to the
gluon fusion. The reason for this is a large value $c_{G}^{\text{eff}}\approx
c_{u}+c_{d}+c_{s}\approx 3$ and the fact that the gluonic squared matrix
element is proportional to $m_{a}^{3}$ rather than to $m_{a}m_{q}^{2}$ as in
the case of the quark fusion, where $m_{q}\ll m_{a}$.

Figure 4: The probability of producing ALPs in the DIS process (14) at various
facilities. The bands show systematics uncertainties, which we obtained by
varying the factorization and renormalization scales
$\mu_{r}^{2}=\mu_{r}^{2}=m_{a}^{2}$ by a factor of two (see text for details).

The cross-section of the gluon fusion at various center-of-mass collision
energies is shown in Fig. 4. It suffers from a large systematical uncertainty,
which we illustrate by varying the renormalization and factorization scales
$\mu_{r}=\mu_{f}=m_{a}$ by a factor of 2. A huge fraction of this uncertainty
comes from the scaling $\sigma_{gg\to a}\propto\alpha_{s}^{2}$, which varies
rapidly for small scales $\mu=\mathcal{O}(1\text{ GeV})$. Another important
question relevant for the LHC is that the production of light ALPs sits at a
very small energy fraction $x=m_{a}^{2}/s_{\text{pp}}$, where the PDFs have a
large uncertainty. Including the uncertainties from PDFs would further
increase the overall error.

Unlike the production via mixing and decays of $B$ mesons, the DIS cross-
section is practically independent of the RG flow. Indeed, the effective gluon
coupling includes the summation over $u,d,s$ quarks. For the ALP masses of
interest, it is proportional to the combination $(c_{u}+c_{d}+c_{s})^{2}$,
which for the wide range of the choices of the scale $\Lambda$ changes
insignificantly.

###  III.4 Discussion

To compare the contribution of the particular production channels to the ALP
yield, we consider the production probabilities per proton collision:

$P_{\text{prod}}=\begin{cases}\chi_{m^{0}}\times|\theta_{am^{0}}|^{2},&\text{mixing}\,,\\\
\chi_{B}\times\text{Br}(B\to X_{s}a),&\text{B decays}\,,\\\
\frac{\sigma_{\text{DIS}}}{\sigma_{pp}},&\text{DIS}\;.\end{cases}$ (15)

Here, $\chi_{X}$ is the fraction of the produced $X$ particles per proton
collision (for $B$, we take both mesons and anti-mesons), and $\sigma_{pp}$ is
the total proton collision cross-section. The mixing angles $\theta_{m^{0}a}$
are taken from eqs. (12), and for the branching ratio $\text{Br}(B\to X_{s}a)$
we used Fig. 2 and Table 1.

We will not consider here the ALP production by decays of kaons, since the
latter are long-lived, and the ALP flux may be heavily affected by the
interaction of the kaons with the infrastructure surrounding the kaon
production point. The most important factor is their absorption by the
material: the absorption length of the kaons is typically smaller than their
decay length $3.7E_{K}/m_{K}\text{ m}$. For the beam dump experiments with
thick targets, kaons would already be heavily absorbed inside the target. For
the LHC-based experiments, the effective kaon decay volume is limited by the
detector, and only a small fraction of kaons would decay there due to their
huge boosts (see, e.g., Beltrán et al. (2023)).444For the experiments with a
thin target like DUNE Abi et al. (2020) and NA62 Hahn et al. (2010), whose
goal is to maximize the kaon flux within the detector acceptance, kaons bypass
the system of magnetic collimators. In such an experiment, one cannot directly
obtain the ALP yield based on the number of protons collisions. Estimating
sensitivities to FIPs, therefore, requires a dedicated study (see, e.g., Ref.
Ovchynnikov et al. (2023b).

Figure 5:  The ALP production probabilities (eq. (15) per proton collision at
the LHC, SPS, and FNAL (DUNE), with the collision energies of
$\sqrt{s}=13\text{ TeV},\approx 28\text{ GeV}$ and $\approx 16\text{ GeV}$
correspondingly. The probabilities are evaluated for the value of the ALP
coupling $f=1\text{ TeV}$ and assuming two scales $\Lambda$: $\Lambda=m_{t}$
(dashed lines) and $\Lambda=1\text{ TeV}$ (solid). The meaning of the vertical
gray bands is the same as in Fig. 3. See text for details.

The total production probabilities for these energies are shown in Fig. 5.
From the comparison, we see that the dominant production channels at the LHC
are decays of $B$ mesons (thanks to the large fraction of the produced
$b\bar{b}$ pairs) and Drell-Yan process, independently of the model scale
$\Lambda$ chosen. This finding differs strongly from the case of ALPs coupled
to gluons Chakraborty et al. (2021); Bauer et al. (2021), mostly related to
the loop suppression of the $bs$ coupling for the latter case.

The mixing with light mesons may become relevant for the experiments operating
at lower energies. Depending on the scale $\Lambda$, it may dominate the total
production at SPS at masses $m_{a}\lesssim 2\text{ GeV}$. It also dominates
the production at FNAL, even at large $\Lambda$, given the small center-of-
mass energy and the correspondingly tiny fraction of produced $B$ mesons.

To conclude this discussion, we emphasize that the hierarchy of the production
channels may change depending on the placement of the experiment with respect
to the proton beam axis. Generically, the angular distributions of the light
ALPs from $B$ decays and those produced by the Drell-Yan process are broader
than the distribution from the mixing with light mesons. Therefore, if these
production channels provide similar overall amounts of ALPs, the mixing with
the mesons would dominate for on-axis experiments with small angular coverage,
while the other channels dominate for off-axis experiments. To demonstrate
this point, in Fig. 6, we show the production probabilities for the ALPs
flying within the polar range $\theta<10\text{ mrad}$ and $\theta>10\text{
mrad}$ at SPS, assuming $\Lambda=300\text{ GeV}$.

Figure 6: The production probabilities of the ALPs flying in the polar range
$\theta<10\text{ mrad}$ (solid lines) and $\theta>10\text{ mrad}$ (dashed
lines) at SPS, assuming the scale $\Lambda=300\text{ GeV}$. By the curves
“mixing”, we summarize the production of the ALPs via the mixing angles with
the mesons $\pi^{0},\eta,$ and $\eta^{\prime}$.

##  IV Decay modes

###  IV.1 Decays into leptons and photons

The matrix element of the decay of ALPs into a pair of photons is given by

$\mathcal{M}_{a\to\gamma\gamma}=\frac{1}{f}\frac{\alpha_{\text{EM}}}{\pi}c^{\text{eff}}_{\gamma}F_{\mu\nu}(p_{\gamma
1})\tilde{F}^{\mu\nu}(p_{\gamma 2})$ (16)

Depending on the ALP mass, the effective coupling $c_{\gamma}^{\text{eff}}$ is

$c_{\gamma}^{\text{eff}}=c_{\gamma}^{\text{loop,l}}+\begin{cases}c_{\gamma}^{\text{VMD}},\quad
m_{a}<m_{a}^{\prime}\\\ c_{\gamma}^{\text{loop,q}},\quad m_{a}\geq
m_{a}^{\prime}\end{cases}$ (17)

where $m^{\prime}_{a}\simeq 2$ GeV is similar to the matching scale between
the ALP’s ChPT and QCD perturbative decays (see the next subsection).

The VMD contribution $c_{\gamma}^{\text{VMD}}$ originates from the mixing of
the vector mesons $\rho,\omega,\phi$ with photons. The corresponding
contributions are

$c_{\gamma}^{\text{VMD}}=-\frac{1}{9}F(m_{a})\big{(}4\sqrt{6}\theta_{\eta
a}+7\sqrt{3}\theta_{\eta^{\prime}a}+9\theta_{\pi^{0}a}\big{)}$ (18)

We have checked that this coupling approximately reproduces the widths of the
anomalous decays $m^{0}\to 2\gamma$ with $m^{0}=\pi^{0}/\eta/\eta^{\prime}$ in
the symbolic limit $\theta_{m^{0^{\prime}}a}=\delta_{m^{0}m^{0^{\prime}}}$,
$f=f_{\pi}$, and $m_{a}=m_{m^{0}}$.

The loop contribution is given by triangle diagrams with fermions $f=l,q$
running inside the loop Bauer et al. (2017, 2022):

$c_{\gamma}^{\text{loop,f}}=\sum_{f}2N_{c}^{f}Q_{f}^{2}c_{f}B_{1}(\tau_{f}),$
(19)

where $N_{c}^{q}=3$, $N_{c}^{l}=1$, $Q_{f}$ is the charge of the fermion, and
$B_{1}$ is from eq. (5).

Decay widths into lepton pairs
$\ell^{+}\ell^{-}=e^{+}e^{-},\mu^{+}\mu^{-},\tau^{+}\tau^{-}$ are described by
the formula

$\Gamma(a\to\ell^{+}\ell^{-})=\frac{c_{\ell}^{2}m_{a}m_{\ell}^{2}\sqrt{1-\frac{4m_{\ell}^{2}}{m_{a}^{2}}}}{2\pi
f^{2}},$ (20)

The total width into leptons and the width into photons are shown in Fig. 7.

###  IV.2 Hadronic decays

Let us now discuss the hadronic decays of ALPs. For
$m_{a}\gg\Lambda_{\text{QCD}}$, it is adequate to describe these decays by
decays into quarks and gluons Bauer et al. (2022):

$\displaystyle\Gamma(a\to q\bar{q})$
$\displaystyle=\frac{N_{c}c_{q}^{2}m_{a}m_{q}^{2}\sqrt{1-\frac{4m_{X_{q}}^{2}}{m_{a}^{2}}}}{2\pi
f^{2}},$ (21) $\displaystyle\Gamma(a\to GG)$
$\displaystyle=|c_{G}^{\text{eff}}|^{2}\frac{\alpha_{s}^{2}\left(1+\frac{83\alpha_{s}}{4\pi}\right)m_{a}^{3}}{8\pi
f^{2}},$ (22)

where $c_{G}^{\text{eff}}$ is given by eq. (4), and $N_{c}=3$ is the number of
colors. To approximately account for hadronization, when considering the
kinematics, we replace the quark’s mass with the mass of the lightest meson
containing the given quark Gunion et al. (2000); Winkler (2019). For instance,
for decays into $c\bar{c}$, $X_{q}$ is a $D$ meson.

Because of the same reason as it was discussed in the context of the Drell-Yan
production (Sec. III.3), decays into gluons dominate over decays into light
quarks $u\bar{u},d\bar{d},s\bar{s}$. However, above the $D$-meson pair
production threshold, the decay into $c\bar{c}$ contributes a sizable fraction
of the ALP total width because of a large $c$ mass.

For $m_{a}=\mathcal{O}(1\text{ GeV})$, perturbative QCD breaks down, and one
should use ChPT. To describe the palette of various decays (e.g.,
$\eta/\eta^{\prime}\to\pi\pi\gamma,\eta^{\prime}\to 4\pi$, or
$\eta^{\prime}\to\eta\pi\pi$) in agreement with the experimental data, the
minimal ChPT is supplemented by phenomenological Lagrangians of the
interactions of the pseudoscalar mesons $P$ with vector Fujiwara et al.
(1985); Guo et al. (2012), scalar Fariborz and Schechter (1999), and tensor
mesons Guo et al. (2012), with the operators and their couplings being fixed
by theoretical arguments (such as the chiral symmetry or anomaly matching
conditions) and to match the experimental data on interactions of $P$. These
mesons may contribute to the matrix elements as intermediate states; one
example is the mixing of neutral vector mesons $\rho^{0},\omega,\phi$ with
photons. The ChPT width should match the parton-level width at some mass
$m_{a}\simeq 1\text{ GeV}$.

Ref. Aloni et al. (2019a) followed this data-driven approach to describe the
decays of the ALPs coupled to gluons; however, their approach also applies to
the ALPs coupled to fermions. Ref. Cheng et al. (2022) repeated the analysis
of Aloni et al. (2019a) with some modifications for the ALPs with a non-
universal explicitly isospin-breaking coupling to quarks.

Figure 7: Summary of decays of the ALPs with the universal coupling to
fermions. Top panel: hadronic decay widths, assuming the scale
$\Lambda=1\text{ TeV}$. For the ALP mass $m_{a}\lesssim 2\text{ GeV}$, ChPT
describes decays of ALPs, while at higher masses, they may be approximated by
the perturbative QCD; the vertical dashed line shows the matching mass. Bottom
panel: the total non-hadronic (including di-$e$,$\mu$,$\tau$, and $\gamma$
processes) and hadronic widths. The solid lines correspond to the scale
$\Lambda=1\text{ TeV}$, while the dashed ones correspond to $\Lambda=m_{t}$.
The vertical gray bands show the vicinity of
$m_{a}=m_{\pi^{0}}/m_{\eta}/m_{\eta^{\prime}}$ where the description via the
ALP-meson mixing (and hence the results for the hadronic widths) breaks down.

In our analysis, we incorporate the ChPT Lagrangian (following the references
above) in the Mathematica notebook accompanying the paper and calculate the
matrix elements and decay widths for various processes (see Appendices B, A
for details). We include the decay channels $a\to\eta
2\pi,\eta^{\prime}2\pi,4\pi$, $3\pi,\gamma 2\pi$, $3\eta$, $\omega 2\pi$, and
$2V$, where $V=\omega,K^{*},\phi$. As a cross-check, we reproduce the results
of the SM decay widths of the mesons $\eta$ and $\eta^{\prime}$ in the limit
when the ALP matches them, i.e., when $\theta_{m^{0}a}=1$, $f=f_{\pi}$, and
$m_{a}=m_{m^{0}}$. We have also qualitatively reproduced the results from
Aloni et al. (2019a) for the model of the ALPs coupled to gluons (see a
discussion in Appendix A).

The summary of the hadronic widths for the ALPs is shown in Fig. 7. The decays
of low-mass ALPs $m_{a}\lesssim 1\text{ GeV}$ are saturated by $a\to 3\pi$,
$a\to\gamma\pi\pi$. At higher masses, decays into $\eta\pi\pi$, $4\pi$, and
$2V$ become the dominant channels. The ChPT width matches with the width in
the perturbative regime at $m_{\text{match}}\simeq 2\text{ GeV}$.

###  IV.3 Discussion

The leptonic, photonic, and hadronic widths are compared in Fig. 7. In Fig. 8,
we show the branching ratios of the ALP decays into various final states. From
the figures, we see that photonic decays are always sub-dominant, while
hadronic widths are irrelevant for the phenomenology of light ALPs with
$m_{a}\lesssim 1\text{ GeV}$, where leptonic decays dominate. For heavier
ALPs, however, decays into hadrons dominate, increasing the total width by up
to a factor of 100. This conclusion is in qualitative agreement with the paper
Domingo (2017), which studied a somewhat different model of a CP-odd scalar.
We emphasize that in the mass range around 1 GeV, the hadronic decay width is
significantly larger than the one obtained in Refs. Bauer et al. (2021);
Ferber et al. (2023) using perturbative QCD.

The decay palette above 1 GeV is also qualitatively similar to the case of the
ALPs coupled to gluons. This is because, for both of these models, the ALPs
have mixing with the three neutral pseudoscalar mesons
$\pi^{0}/\eta/\eta^{\prime}$.

Interestingly, the choice of the scale $\Lambda$ practically does not
influence the decay phenomenology. This is because it affects only the decays
where the mixing with pions dominates among the others. These are the decays
into $3\pi$ and $\gamma\pi\pi$, which are important only in the mass range
$m_{a}\lesssim 1\text{ GeV}$ where leptonic decay widths are much larger.

Figure 8: Left panel: branching ratios of the ALP decays into important final
states (see Fig. 7 for the description of the states). Right panel: the ALP
lifetime assuming $f=1\text{ TeV}$. The red and blue lines show,
correspondingly, the lifetime in the approximation of using only leptonic
width as in Beacham et al. (2020) and our result (both assuming
$\Lambda=1\text{ TeV}$), while the green line the results from Ferber et al.
(2023) evaluated for the scale $\Lambda=4\pi\text{ TeV}$

To further stress the importance of the hadronic decays and finalize the
discussion, we show the mass dependence of the ALP lifetime as computed in
this work and the one widely considered in the past Beacham et al. (2020),
when only leptonic decays have been included. In the mass range
$2m_{\mu}<m_{a}<2m_{\tau}$, the full width in the description from Beacham et
al. (2020) is saturated by dimuon decay. With hadronic decays being included,
the branching ratio $\text{Br}(a\to\mu\mu)$ is $<10\%$ for masses
$m_{a}\gtrsim 1\text{ GeV}$.

In Fig. 8, we show the branching ratios of the ALP decays into various final
states (left panel) and the ALP lifetime (right panel). For the lifetime, we
compare the predictions assuming the revised phenomenology and the description
from past works: Ref. Beacham et al. (2020), which neglects hadronic decay
modes and is widely used by the experiments community to derive constraints
and sensitivities, and Ref. Ferber et al. (2023), which, following Bauer et
al. (2021, 2022), approximates the hadronic decays by the decay $a\to 3\pi$ at
$m_{a}\lesssim 1\text{ GeV}$ and by the decays into a gluon and quark pairs
above this mass. The lifetime from Beacham et al. (2020) is always much larger
for $m_{a}\gtrsim 1\text{ GeV}$. The lifetime from Ferber et al. (2023)
coincides with our result for the ranges $m_{a}\lesssim 1\text{ GeV}$,
$m_{\text{match}}<m_{a}<2m_{c}$, and $m_{a}\gtrsim 2m_{D}$. The origin of the
discrepancy in the range $2m_{c}<m_{a}\lesssim 4\text{ GeV}$ is that Ref.
Ferber et al. (2023) turns on the decay $a\to c\bar{c}$ above $2m_{c}$, even
though this decay is kinematically impossible until the $2D$ threshold.

Figure 9: Re-interpretation of the model-independent LHCb constraints from the
searches $B\to K+(\text{FIP}\to)\mu\mu$ reported in Aaij et al. (2015, 2017)
for the model of ALPs with the fermion coupling (2), assuming $\Lambda=1\text{
TeV}$, for the plane $m_{a}-2v_{h}/f$, where $v_{h}=246\text{ GeV}$ is the
Higgs VEV. The blue solid line: constraints assuming the ALP phenomenology
obtained in this work. The red solid line: if assuming the phenomenology from
Beacham et al. (2020), which only includes leptonic decays. The light lines of
the same colors show the ALP decay lengths
$c\tau_{a}\langle\gamma_{a}\rangle=1\text{ cm}$ and $1\text{ m}$ assuming the
corresponding phenomenology.

The previously neglected production and decay modes are expected to
significantly change the landscape of the past constraints and future searches
for ALPs. For example, let us consider searches for $B\to K(a\to)\mu\mu$
performed at LHCb Aaij et al. (2015, 2017). It is sensitive to the total ALP
decay width as well as to the branching ratio $\text{Br}(a\to\mu\mu)$. The
constraints to ALPs are shown in Fig. 9, where, for comparison, we display the
bound obtained assuming the ALP phenomenology description from Ref. Beacham et
al. (2020) and the one obtained in this work. The updated constraints are
significantly weaker in the domain of large masses $m_{a}\gtrsim 1\text{
GeV}$. Interestingly, for the revised phenomenology, the lower bound of the
constraint lies in the regime where ALPs are short-lived and mainly decay
within the detector, whereas for the old phenomenology, it mainly belongs to
the parameter space of long-lived ALPs. We will revise further existing
constraints, consider the ones previously not accounted for in the literature,
and derive sensitivities for future experiments in a forthcoming work Garcia
et al. (2023).

##  V Conclusions

The model of ALPs universally coupled to fermions is considered by the physics
beyond colliders (PBC) group as one of the benchmark models to test the
potential of various experiments to explore the parameter space of feebly-
interacting particles. Therefore, understanding the phenomenology of such ALPs
is an important and timely question. In this work, we have revised the
production and decay modes of ALPs at hadronic accelerator experiments,
considering also the impact of the renormalization group (RG) flow of their
couplings depending on the scale $\Lambda$ at which the model is defined (see
Sec. II and in particular Fig. 1).

For the production (Sec. III), we have considered decays of kaons and $B$
mesons, the mixing with neutral mesons $\pi^{0}/\eta/\eta^{\prime}$, and the
Drell-Yan process. For the production via mixing, we have found that the RG
flow is very important, sizably changing the mixing angle squared between the
ALPs and $\pi^{0}$ (Fig. 3). For the production from $B$ mesons, we have
included the decays $B\to X_{s}+a$, with $X_{s}$ being heavy kaon resonances
$K^{*0},K_{1},K_{2},\dots$, which have not been considered previously in this
context in the literature and increase the total production branching ratio by
a factor of 3–4 for light ALPs (Fig. 2). Our results apply also to generic
ALPs, provided that the low-energy Lagrangian describing the decay has the
same operator expression. For the Drell-Yan production, we have considered the
leading-order fusion processes and shown that the cross-section suffers from a
large systematic uncertainty (Fig. 4).

Depending on the scale $\Lambda$ at which the ALP model (2) is defined and the
collision energy, we have found that any of these processes may dominate the
production (Fig. 5). In particular, at DUNE collision energy, the production
via mixing is the main production channel of the ALPs with mass below 2–3 GeV,
while at larger masses, decays of $B$ mesons are the main channel. At SPS
energies, the hierarchy of production channels depends heavily on the scale
$\Lambda$. Namely, if $\Lambda$ is close to the EW scale, the main channels
are mixing with mesons and the Drell-Yan process. Once $\Lambda$ departs from
the EW scale, decays via $B$ mesons dominate. The main production channel may
also depend on the geometric placement of the experiment (see Fig. 6).
Finally, at the LHC, $\Lambda$ practically does not influence the hierarchy
and affects only the magnitude of the ALP flux.

We have studied the full palette of the ALP decays (Sec. IV), including the
hadronic ones that were missing previously. In particular, we have found (Fig.
7 and also Fig. 8) that leptonic decays are the main channels for light ALPs
with $m_{a}\lesssim 1\text{ GeV}$, while the hadronic decays dominate at
higher masses, increasing the total ALP width by up to a factor of 100.
Contrary to the production case, the decay widths are only weakly sensitive to
the choice of $\Lambda$.

To simplify the use of our results by the community, we have implemented the
ALP phenomenology studied in this work in a Mathematica notebook accompanying
the paper. We have also implemented the model in SensCalc Ovchynnikov et al.
(2023a) – a public code evaluating sensitivities of different experiments.

## Acknowledgements

We thank Pilar Coloma, Torben Ferber, Joerg Jaeckel, Jan Jerhot, Felix Kling,
Thomas Schwetz, Yotam Soreq and Susanne Westhoff for helpful discussions. This
work was partially funded by the Deutsche Forschungsgemeinschaft (DFG) through
the Emmy Noether Grant No. KA 4662/1-2 and grant 396021762 – TRR 257. MO
received support from the European Union’s Horizon 2020 research and
innovation program under the Marie Sklodowska-Curie grant agreement No.
860881-HIDDeN. GG thanks the Doctoral School “Karlsruhe School of Elementary
and Astroparticle Physics: Science and Technology (KSETA)” for financial
support through the GSSP program of the German Academic Exchange Service
(DAAD).

## References

  * Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Jaeckel and Ringwald (2010) J. Jaeckel and A. Ringwald, Ann. Rev. Nucl. Part. Sci. 60, 405 (2010), eprint 1002.0329.
  * Georgi et al. (1986) H. Georgi, D. B. Kaplan, and L. Randall, Phys. Lett. B 169, 73 (1986).
  * Bauer et al. (2017) M. Bauer, M. Neubert, and A. Thamm, JHEP 12, 044 (2017), eprint 1708.00443.
  * Brivio et al. (2017) I. Brivio, M. B. Gavela, L. Merlo, K. Mimasu, J. M. No, R. del Rey, and V. Sanz, Eur. Phys. J. C 77, 572 (2017), eprint 1701.05379.
  * Beacham et al. (2020) J. Beacham et al., J. Phys. G 47, 010501 (2020), eprint 1901.09966.
  * Antel et al. (2023) C. Antel et al. (2023), eprint 2305.01715.
  * Aloni et al. (2019a) D. Aloni, Y. Soreq, and M. Williams, Phys. Rev. Lett. 123, 031803 (2019a), eprint 1811.03474.
  * Aloni et al. (2019b) D. Aloni, C. Fanelli, Y. Soreq, and M. Williams, Phys. Rev. Lett. 123, 071801 (2019b), eprint 1903.03586.
  * Chakraborty et al. (2021) S. Chakraborty, M. Kraus, V. Loladze, T. Okui, and K. Tobioka, Phys. Rev. D 104, 055036 (2021), eprint 2102.04474.
  * Döbrich et al. (2016) B. Döbrich, J. Jaeckel, F. Kahlhoefer, A. Ringwald, and K. Schmidt-Hoberg, JHEP 02, 018 (2016), eprint 1512.03069.
  * Dolan et al. (2017) M. J. Dolan, T. Ferber, C. Hearty, F. Kahlhoefer, and K. Schmidt-Hoberg, JHEP 12, 094 (2017), [Erratum: JHEP 03, 190 (2021)], eprint 1709.00009.
  * Döbrich et al. (2019a) B. Döbrich, J. Jaeckel, and T. Spadaro, JHEP 05, 213 (2019a), [Erratum: JHEP 10, 046 (2020)], eprint 1904.02091.
  * Izaguirre et al. (2017) E. Izaguirre, T. Lin, and B. Shuve, Phys. Rev. Lett. 118, 111802 (2017), eprint 1611.09355.
  * Alonso-Álvarez et al. (2019) G. Alonso-Álvarez, M. B. Gavela, and P. Quilez, Eur. Phys. J. C 79, 223 (2019), eprint 1811.05466.
  * Gavela et al. (2019) M. B. Gavela, R. Houtz, P. Quilez, R. Del Rey, and O. Sumensari, Eur. Phys. J. C 79, 369 (2019), eprint 1901.02031.
  * Dolan et al. (2015) M. J. Dolan, F. Kahlhoefer, C. McCabe, and K. Schmidt-Hoberg, JHEP 03, 171 (2015), [Erratum: JHEP 07, 103 (2015)], eprint 1412.5174.
  * Döbrich et al. (2019b) B. Döbrich, F. Ertas, F. Kahlhoefer, and T. Spadaro, Phys. Lett. B 790, 537 (2019b), eprint 1810.11336.
  * Carmona et al. (2021) A. Carmona, C. Scherb, and P. Schwaller, JHEP 08, 121 (2021), eprint 2101.07803.
  * Cornella et al. (2020) C. Cornella, P. Paradisi, and O. Sumensari, JHEP 01, 158 (2020), eprint 1911.06279.
  * Calibbi et al. (2021) L. Calibbi, D. Redigolo, R. Ziegler, and J. Zupan, JHEP 09, 173 (2021), eprint 2006.04795.
  * Chala et al. (2021) M. Chala, G. Guedes, M. Ramos, and J. Santiago, Eur. Phys. J. C 81, 181 (2021), eprint 2012.09017.
  * Bauer et al. (2021) M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, JHEP 04, 063 (2021), eprint 2012.12272.
  * Bauer et al. (2022) M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, JHEP 09, 056 (2022), eprint 2110.10698.
  * Ferber et al. (2023) T. Ferber, A. Filimonova, R. Schäfer, and S. Westhoff, JHEP 04, 131 (2023), eprint 2201.06580.
  * Ertas and Kahlhoefer (2020) F. Ertas and F. Kahlhoefer, JHEP 07, 050 (2020), eprint 2004.01193.
  * Jerhot et al. (2022) J. Jerhot, B. Döbrich, F. Ertas, F. Kahlhoefer, and T. Spadaro, JHEP 07, 094 (2022), eprint 2201.05170.
  * Liu et al. (2023) J. Liu, Y. Luo, and M. Song, JHEP 09, 104 (2023), eprint 2304.05435.
  * Bruggisser et al. (2023) S. Bruggisser, L. Grabitz, and S. Westhoff (2023), eprint 2308.11703.
  * Boiarska et al. (2019) I. Boiarska, K. Bondarenko, A. Boyarsky, V. Gorkavenko, M. Ovchynnikov, and A. Sokolenko, JHEP 11, 162 (2019), eprint 1904.10447.
  * Ovchynnikov et al. (2023a) M. Ovchynnikov, J.-L. Tastet, O. Mikulenko, and K. Bondarenko (2023a), eprint 2305.13383.
  * Batell et al. (2011) B. Batell, M. Pospelov, and A. Ritz, Phys. Rev. D 83, 054005 (2011), eprint 0911.4938.
  * Grzadkowski and Krawczyk (1983) B. Grzadkowski and P. Krawczyk, Z. Phys. C 18, 43 (1983).
  * Leutwyler and Shifman (1990) H. Leutwyler and M. A. Shifman, Nucl. Phys. B 343, 369 (1990).
  * Haber et al. (1987) H. E. Haber, A. S. Schwarz, and A. E. Snyder, Nucl. Phys. B 294, 301 (1987).
  * Chivukula and Manohar (1988) R. S. Chivukula and A. V. Manohar, Phys. Lett. B 207, 86 (1988), [Erratum: Phys.Lett.B 217, 568 (1989)].
  * Cheng et al. (2022) H.-C. Cheng, L. Li, and E. Salvioni, JHEP 01, 122 (2022), eprint 2110.10691.
  * Alwall et al. (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), eprint 1405.0301.
  * Alloul et al. (2014) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, Comput. Phys. Commun. 185, 2250 (2014), eprint 1310.1921.
  * Christensen and Duhr (2009) N. D. Christensen and C. Duhr, Comput. Phys. Commun. 180, 1614 (2009), eprint 0806.4194.
  * Berlin and Kling (2019) A. Berlin and F. Kling, Phys. Rev. D 99, 015021 (2019), eprint 1810.01879.
  * Beltrán et al. (2023) R. Beltrán, J. Günther, M. Hirsch, A. Titov, and Z. S. Wang (2023), eprint 2309.11546.
  * Abi et al. (2020) B. Abi et al. (DUNE), JINST 15, T08010 (2020), eprint 2002.03010.
  * Hahn et al. (2010) F. Hahn, F. Ambrosino, A. Ceccucci, H. Danielsson, N. Doble, F. Fantechi, A. Kluge, C. Lazzeroni, M. Lenti, G. Ruggiero, et al. (NA62), Tech. Rep., CERN, Geneva (2010), URL https://cds.cern.ch/record/1404985.
  * Ovchynnikov et al. (2023b) M. Ovchynnikov, T. Schwetz, and J.-Y. Zhu, Phys. Rev. D 107, 055029 (2023b), eprint 2210.13141.
  * Gunion et al. (2000) J. F. Gunion, H. E. Haber, G. L. Kane, and S. Dawson, _The Higgs Hunter’s Guide_ , vol. 80 (2000).
  * Winkler (2019) M. W. Winkler, Phys. Rev. D 99, 015018 (2019), eprint 1809.01876.
  * Fujiwara et al. (1985) T. Fujiwara, T. Kugo, H. Terao, S. Uehara, and K. Yamawaki, Prog. Theor. Phys. 73, 926 (1985).
  * Guo et al. (2012) F.-K. Guo, B. Kubis, and A. Wirzba, Phys. Rev. D 85, 014014 (2012), eprint 1111.5949.
  * Fariborz and Schechter (1999) A. H. Fariborz and J. Schechter, Phys. Rev. D 60, 034002 (1999), eprint hep-ph/9902238.
  * Domingo (2017) F. Domingo, JHEP 03, 052 (2017), eprint 1612.06538.
  * Aaij et al. (2015) R. Aaij et al. (LHCb), Phys. Rev. Lett. 115, 161802 (2015), eprint 1508.04094.
  * Aaij et al. (2017) R. Aaij et al. (LHCb), Phys. Rev. D 95, 071101 (2017), eprint 1612.07818.
  * Garcia et al. (2023) G. D. V. Garcia, F. Kahlhoefer, M. Ovchynnikov, and A. Zaporozhchenko (2023), eprint to appear.
  * Krauss and Wise (1986) L. M. Krauss and M. B. Wise, Phys. Lett. B 176, 483 (1986).
  * Bardeen et al. (1987) W. A. Bardeen, R. D. Peccei, and T. Yanagida, Nucl. Phys. B 279, 401 (1987).
  * Srednicki (1985) M. Srednicki, Nucl. Phys. B 260, 689 (1985).

##  Appendix A ChPT with ALPs

Let us, for completeness, assume that both the quark couplings $c_{q}$ from
eq. (2) and the gluon coupling $c_{G}$ in eq. (3) are present in the
Lagrangian. Both of them may contribute to the ChPT interactions. To this end,
let us first convert the gluon coupling in eq. (3) to the pure quark sector by
performing the following chiral rotation of the light quarks $q=(u,d,s)$
Krauss and Wise (1986); Bardeen et al. (1987); Srednicki (1985); Georgi et al.
(1986):

$q\to\mathcal{U}q,\quad\mathcal{U}=\exp\left[-ic_{G}\frac{a}{f}\kappa_{q}\gamma_{5}\right]\;,$
(23)

where $(\kappa_{q})_{ij}=\delta_{ij}m_{q_{i}}^{-1}/(m_{u}+m_{d}+m_{s})$ is
fixed in order to prevent any mass mixing between the mesons $\pi^{0}$ and
$\eta_{0}$ with $a$. The hadronic part of the Lagrangian becomes

$\langle\mathcal{L}_{\text{model
1,\text{hadr}}}\rangle_{c,b,t}=\bar{q}\,\mathcal{U}m\,\mathcal{U}q+\frac{\partial^{\mu}a}{f}(c_{q}+c_{G}\kappa_{q})\bar{q}\gamma_{\mu}\gamma_{5}q\;,$
(24)

where we have neglected the off-diagonal quark coupling $c_{sd}$ generated by
integrating out the top quark due to its strong CKM suppression. The relevant
ChPT Lagrangian then is Bauer et al. (2021); Aloni et al. (2019a)

$\mathcal{L}_{\text{ChPT,min}}=\frac{1}{2}(\partial_{\mu}a)^{2}-\frac{m_{a}^{2}}{2}a^{2}+\frac{f_{\pi}^{2}}{2}B_{0}\text{Tr}\left[\Sigma\hat{m}^{\dagger}_{q}+\hat{m}_{q}\Sigma^{\dagger}\right]+\frac{f_{\pi}^{2}}{4}\text{Tr}\left[D_{\mu}\Sigma
D^{\mu}\Sigma^{\dagger}\right]+\\\
\frac{f_{\pi}^{2}}{2}\frac{\partial_{\mu}a}{f}\text{Tr}[(\hat{c}_{q}+c_{G}\kappa_{q})(\Sigma
D^{\mu}\Sigma^{\dagger}-\Sigma^{\dagger}D^{\mu}\Sigma)]\;,$ (25)

where $c_{q}=\text{diag}(c_{u},c_{d},c_{s})$,

$\hat{m}_{q}=\exp\left[-ic_{G}\frac{a}{f}\kappa_{q}\right]m_{q}\exp\left[-ic_{G}\frac{a}{f}\kappa_{q}\right]\;,$
(26)

$\Sigma$ is the matrix of the pseudoscalar mesons

$\Sigma=\exp\left[\frac{2i\mathcal{P}}{f_{\pi}}\right],\quad\mathcal{P}=\frac{1}{\sqrt{2}}\begin{pmatrix}\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{3}}+\frac{\eta^{\prime}}{\sqrt{6}}&\pi^{+}&K^{+}\\\
\pi^{-}&-\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{3}}+\frac{\eta^{\prime}}{\sqrt{6}}&K^{0}\\\
K^{-}&\bar{K}^{0}&-\frac{\eta}{3}+2\frac{\eta^{\prime}}{6},\end{pmatrix}$ (27)

$D_{\mu}\Sigma=\partial_{\mu}\Sigma+ieA_{\mu}\Sigma$ is the covariant
derivative.

We also need to include the phenomenological Lagrangian of the interactions of
pseudoscalar mesons with other mesons: anomalous WZW interactions and
interactions with vector Fujiwara et al. (1985); Guo et al. (2012), scalar
Fariborz and Schechter (1999), and tensor meson $f_{2}$ Cheng et al. (2022)
(see also Guo et al. (2012)):

$\displaystyle\mathcal{L}_{\text{vec+an}}=$
$\displaystyle-\frac{3g^{2}}{8\pi^{2}f_{\pi}}\epsilon^{\mu\nu\alpha\beta}\text{Tr}[P(x)d_{\mu}V_{\nu}(x)d_{\alpha}V_{\beta}(x)]+\frac{7}{60\pi^{2}f_{\pi}^{5}}\epsilon^{\mu\nu\alpha\beta}\text{Tr}[P(x)\partial_{\mu}P\partial_{\nu}P\partial_{\alpha}P\partial_{\beta}P]$
(28)
$\displaystyle+2f_{\pi}^{2}\text{Tr}\left|gV_{\mu}-eA_{\mu}Q-\frac{i}{2f_{\pi}^{2}}[P,\partial_{\mu}P]\right|^{2}$
(29) $\displaystyle-
g_{f_{2}\pi\pi}\frac{f_{\pi}^{2}}{4}\text{Tr}\left[\left(\partial^{\mu}\Sigma^{\dagger}\partial^{\nu}\Sigma-\frac{1}{2}g^{\mu\nu}\partial^{\alpha}\Sigma^{\dagger}\partial_{\alpha}\Sigma\right)\bm{f}_{\mathbf{2}}\right]\phi_{\mu\nu}+\mathcal{L}_{\text{scalar}}$
(30)

Here, $g\approx m_{\rho}/\sqrt{2}f_{\pi}$, $Q=\text{diag}[2/3,-1/3,-1/3]$ is
the quark charge matrix, $V_{\mu}$ is the matrix of vector mesons,

$V_{\mu}=\frac{1}{\sqrt{2}}\begin{pmatrix}\frac{\rho^{0}+\omega}{\sqrt{2}}&\rho^{+}&K^{*+}\\\
\rho^{-}&\frac{-\rho^{0}+\omega}{\sqrt{2}}&K^{*0}\\\
K^{*-}&\bar{K}^{*0}&\phi\end{pmatrix}$ (31)

and $A_{\mu}$ is the EM field. Next, $\mathcal{L}_{\text{scalar}}$ is given by
eq. (A1) from Fariborz and Schechter (1999). The tensor meson $f_{2}$ is
denoted by $\phi_{\mu\nu}$, while $\bm{f_{2}}$ is the SU(3) generator of the
tensor meson. The coupling $g_{f_{2}\pi\pi}=13.1\text{ GeV}^{-1}$ Cheng et al.
(2022).

Having the Lagrangians (25), (30), we calculate the various contributions to
the matrix elements of the ALP production and decay (see Appendix B for the
impementation).

When calculating the decay matrix elements, we mostly follow the assumptions
considered in Cheng et al. (2022) based on observational data and unitarity
requirements. As an example, we artificially set the contact VMD terms
originating from the square of the last summand of the second line of Eq. (30)
to zero if $m_{a}>m_{\eta^{\prime}}$.

###  A.1 Comparison with ref. Aloni et al. (2019a)

There are some differences in the description of the ALP decays from Aloni et
al. (2019a) and Cheng et al. (2022) (and hence our approach). For instance,
unlike Cheng et al. (2022), ref. Aloni et al. (2019a) does not include the
contributions of $\mathcal{L}_{\text{vec+an}}$ to the decays $a\to 3\pi$,
which changes the corresponding width by orders of magnitude. It is crucial
since this width dominates the ALP decays in the mass range $m_{a}\lesssim
1\text{ GeV}$. Another difference is that the vector meson contribution to the
widths $a\to KK\pi$ is included in Cheng et al. (2022) but not in Aloni et al.
(2019a).

Yet another difference is in the sector of the interactions with scalar
mesons. We have used the interaction Lagrangian directly from Appendix A of
Fariborz and Schechter (1999), which assumes the SU(3) representation of the
scalar mesons given by eqs. (1.2), (1.3). Ref. Aloni et al. (2019a), using the
same reference, represented the ALP decay matrix elements in terms of SU(3)
representation of the scalar mesons that directly contradicts the definitions
(1.2), (1.3).

As a result of these differences, our prediction of the decays of ALPs coupled
solely to gluons differs from the one presented in Aloni et al. (2019a). In
particular, we have found a somewhat larger value of the mass where the ChPT
width matches with the perturbative QCD width, $m_{a}=2.3\text{ GeV}$.

##  Appendix B Mathematica notebook

To calculate and summarize the ALP production and decay rates, we implement
the Lagrangian (25), (30), as well as the RG flow for the couplings
$\\{c_{q}\\}=c_{u,d,s,c,b,t}$, $c_{l}=c_{e,\mu,\tau}$ in a Mathematica
notebook555Available on https://github.com/maksymovchynnikov/ALPs-
phenomenology. The structure of the notebook is as follows. First, we define
several ALP models at a scale $\Lambda$, such as the ALPs with universal
fermion and gluon couplings. Then, we solve the RG equations for the fermion
couplings (both the diagonal and FCNC couplings) at various scales $\Lambda$
following Bauer et al. (2021, 2022), and interpolate the solutions.

Next, we implement the ChPT Lagrangian (25), (30) keeping the arbitrary values
of the couplings $c_{u,d,s}$, and $c_{G}$. We diagonalize the quadratic ChPT
Lagrangian to get the mixing angles and the ALP interactions. Then, we define
the Feynman rules for the obtained Lagrangian and compute the matrix elements
and decay widths of the ALP decay processes listed in Sec. IV. In the last
step, we specify the model and the scale $\Lambda$ and insert the resulting
couplings into the decay widths. The resulting tabulated widths are then
exported.

Finally, we compute the total ALP production rates described in Sec. III. For
this, we again use the RG flow of the couplings and the pre-computed ALP
production cross-section in the gluon fusion for the unit value of $c_{G}$.

For ALPs coupled solely to $W$ and $B$ bosons, only the production rates are
currently evaluated. We will fully implement these models in future versions
of the notebook.