# Higgs Footprints of Hefty ALPs

Anisha<EMAIL_ADDRESS>School of Physics and Astronomy, University of
Glasgow, Glasgow G12 8QQ, United Kingdom
Supratim Das Bakshi<EMAIL_ADDRESS>CAFPE and Departamento de Física Teórica y del
Cosmos, Universidad de Granada, Campus de Fuentenueva, E–18071 Granada, Spain
Christoph Englert<EMAIL_ADDRESS>School of Physics and
Astronomy, University of Glasgow, Glasgow G12 8QQ, United Kingdom
Panagiotis Stylianou<EMAIL_ADDRESS>Deutsches Elektronen-
Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

###### Abstract

We discuss axion-like particles (ALPs) within the framework of Higgs Effective
Field Theory, targeting instances of close alignment of ALP physics with a
custodial singlet character of the Higgs boson. We tension constraints arising
from new contributions to Higgs boson decays against limits from high-momentum
transfer processes that become under increasing control at the LHC. Going
beyond leading-order approximations, we highlight the importance of multi-top
and multi-Higgs production for the pursuit of searches for physics beyond the
Standard Model extensions.

###### pacs:

††preprint: DESY-23-082

##  I Introduction

Searches for new interactions beyond the Standard Model of Particle Physics
have, so far, been unsuccessful. This is puzzling as the Standard Model
contains a plethora of flaws that are expected to be addressed by a more
comprehensive theory of microscopic interactions. A reconciliation of these
flaws can have direct phenomenological consequences for physics at or below
the weak scale $v\simeq 246~{}\text{GeV}$. This is particularly highlighted by
fine-tuning problems related to the Higgs mass or the neutron electric dipole
moment, both of which take small values due to cancellations which are not
protected by symmetries in the Standard Model (SM). Dynamical solutions to
these issues have a long history, leading to new interactions and states
around the TeV scale to address Higgs naturalness, or relaxing into and CP-
conserving QCD vacuum via a Peccei-Quinn-like mechanism Peccei and Quinn
(1977). Often such approaches yield an additional light pseudo-Nambu Goldstone
field, in the guise of a composite Higgs boson or the axion Wilczek (1978).

The search for a wider class of the latter states referred to as axion-like
particles (ALPs) bridges different areas of high energy physics. Efforts to
detect ALPs across different mass and coupling regimes have shaped the current
BSM programme in many experimental realms (see e.g. Graham et al. (2015);
Irastorza and Redondo (2018) for recent reviews). In particular, at the Large
Hadron Collider (LHC), ALP interactions have been discussed in relation to
their tell-tale signatures arising from $\sim F\tilde{F}$ coupling structures
Jaeckel and Spannowsky (2016), top quarks Esser et al. (2023), emerging
signatures Alonso-Álvarez et al. (2023), flavour physics Bauer et al. (2021,
2022), electroweak precision constraints Bauer et al. (2017), Higgs decays
Bauer et al. (2017); Biekötter et al. (2022), and mixing Bauer et al. (2023).
The methods of effective field theory Bonilla et al. (2021) naturally embed
ALP-related field theories into a broader framework of a more modern
perspective on renormalisability Buchalla et al. (2014, 2018, 2021).

Experimental searches for these states have been carried out using a variety
of techniques, including collider searches, precision measurements of atomic
and nuclear transitions (e.g. ACME ACME Collaboration (2018) and nEDM Abel et
al. (2020)), and searches from astrophysical events Raffelt (1996); Raffelt
and Stodolsky (1988), over a wide range of ALP mass O’Hare (2020). In
particular for the ALP mass range $M_{\mathcal{A}}\,\in$ [6, 100] GeV, the
most stringent exclusion limits for ALPs are derived from ultra-peripheral
lead nuclei collision data Aad et al. (2021a); Sirunyan et al. (2019). These
limits are from exclusive di-photon searches, and define SM–ALPs interactions
via electromagnetic interactions $(\sim F\widetilde{F})$. 111When the ALP mass
$M_{\mathcal{A}}<2m_{e}$, only the di-photon channel is the allowed decay
process via SM particles. In same manner, with greater ALP mass, the decay
modes to other leptons, quarks (jets), gauge bosons open up as well. The ATLAS
limits Aad et al. (2021a) on the ALP–photon cross sections when put in terms
of ALP–photon couplings is found in the range
$g_{\mathcal{A}\gamma}\in[0.05,1]$ TeV-1 d’Enterria (2021). However, in
general, the ALP–SM interactions can be defined via gauge bosons, fermions,
and scalars; although, its decays will depend on its mass. The limits on ALP
couplings to the SM fields (except photons) are less stringent. The exotic
decays of the SM Higgs and Z-bosons are promising channels for ALP searches
(particularly benefiting from the high-luminosity run of the LHC), e.g., with
the decay modes $h\rightarrow Z\mathcal{A}$ Sirunyan et al. (2023).

It is the latter perspective that we adopt in this note to focus on ALP
interactions with the Higgs boson, also beyond leading order. Adopting the
methodology of Higgs Effective Theory (HEFT), we can isolate particular
interactions of the ALP state and trace their importance (and thus the
potential for constraints) to representative collider processes that navigate
between the low energy precision and the large momentum transfer regions
accessible at the LHC. If the interactions of ALPs and Higgs particles is
predominantly related to a custodial singlet realisation of the Higgs boson,
these areas might well be the first phenomenological environments where BSM
could be unveiled as pointed out in, e.g., Ref. Brivio et al. (2017). In
parallel, our results demonstrate the further importance of multi-top and
multi-Higgs final states as promising candidates for the discovery of new
physics. With the LHC experiments closing in on both Higgs pair Tumasyan et
al. (2023); Aad et al. (2023a) and four-top production in the SM Hayrapetyan
et al. (2023); Aad et al. (2023b), such searches becoming increasingly
interesting for our better understanding of the BSM landscape.

This work is organised as follows: In Sec. II, we review the ALP-HEFT
framework that we use in this study to make this work self-contained (a
comprehensive discussion is presented in Brivio et al. (2017)). In Sec. III,
we focus on the decay phenomenology of the Higgs boson in the presence of ALP
interactions before we turn to discuss a priori sensitive processes that can
provide additional constraints due to their multi-scale nature and kinematic
coverage. Specifically, in Sec. IV we analyse ALP corrections to Higgs
propagation as accessible in four-top final states Englert et al. (2019);
Sirunyan et al. (2020); Aad et al. (2023b), which informs corrections to
multi-Higgs production. We conclude in Sec. V.

##  II ALP Chiral HEFT Lagrangian

The leading order ALP interactions with SM fields in the framework of chiral
(non-linear) electroweak theory are written as

$\mathcal{L}_{\text{LO}}=\mathcal{L}^{\text{HEFT}}_{\text{LO}}+\mathcal{L}^{\text{ALP}}_{\text{LO}}\,.$
(1)

$\mathcal{L}^{\text{HEFT}}_{\text{LO}}$ is the chiral dimension-2 HEFT
Lagrangian Alonso et al. (2013); Brivio et al. (2016); Herrero and Morales
(2020, 2021); Anisha et al. (2022). In this framework, the SM Higgs ($H$) is a
singlet field and the Goldstone bosons $\pi^{a}$ are parametrised non-linearly
using the matrix $U$

$U(\pi^{a})=\exp\left({i\pi^{a}\tau^{a}/v}\right)\,,$ (2)

with $\tau^{a}$ as the Pauli matrices with $a=1,2,3$ and $v\simeq
246~{}\text{GeV}$. The $U$ matrix transforms under $L\in
SU(2)_{L},U(1)_{Y}\subset SU(2)_{R}\ni R$ as $U\to LUR^{\dagger}$ and is
expanded as

$U(\pi^{a})={\mathbbm{1}}_{2}+i\frac{\pi^{a}}{v}\tau^{a}-\frac{2G^{+}G^{-}+G^{0}G^{0}}{2v^{2}}{\mathbbm{1}}_{2}+\dots\,,$
(3)

where $G^{\pm}=(\pi^{2}\pm i\pi^{1})/\sqrt{2}$ and $G^{0}=-\pi^{3}$. The
dynamics of the gauge bosons $W^{a}_{\mu}$ and $B_{\mu}$ are determined by the
usual $SU(2)_{L}\times U(1)_{Y}$ gauge symmetry. Weak gauging of
$SU(2)_{L}\times U(1)_{Y}$ is achieved through the standard covariant
derivative

$D_{\mu}U=\partial_{\mu}U+ig_{W}(W^{a}_{\mu}\tau^{a}/2)\;U-ig^{\prime}UB_{\mu}\tau^{3}/2\,.$
(4)

The gauge fields in the physical (mass and electromagnetic $U(1)_{\text{em}}$)
basis are are related to the gauge basis via the Weinberg angle
$s_{W}=\sin\theta_{W},c_{W}=\cos\theta_{W}$

$\begin{split}W^{\pm}_{\mu}&=\frac{1}{\sqrt{2}}(W^{1}_{\mu}\mp
W^{2}_{\mu})\,,\\\ \begin{pmatrix}Z_{\mu}\\\
A_{\mu}\end{pmatrix}&=\begin{pmatrix}c_{W}&s_{W}\\\
-s_{W}&c_{W}\end{pmatrix}\begin{pmatrix}W^{3}_{\mu}\\\
B_{\mu}\end{pmatrix}\,.\end{split}$ (5)

The leading order HEFT Lagrangian relevant for our discussion is then given by

$\mathcal{L}^{\text{HEFT}}_{\text{LO}}=-\frac{1}{4}W^{a}_{\mu\nu}W^{a\mu\nu}-\frac{1}{4}B_{\mu\nu}B^{\mu\nu}+{\cal{L}}_{\text{ferm}}+\mathcal{L}_{\text{Yuk}}\\\
+\frac{v^{2}}{4}\mathcal{F}_{H}\,\text{Tr}[D_{\mu}U^{\dagger}D^{\mu}U]\\\
+\frac{1}{2}\partial_{\mu}H\partial^{\mu}H-V(H)\,.$ (6a) The interactions of
the singlet Higgs field with gauge and Goldstone bosons are parametrised by
the flare function $\mathcal{F}_{H}$ given as
$\mathcal{F}_{H}=\Big{(}1+2(1+\zeta_{1})\frac{H}{v}+(1+\zeta_{2})\Big{(}\frac{H}{v}\Big{)}^{2}+...\Big{)}\,.$
(6b)

The couplings $\zeta_{i}$ denote the independent parameters that determine the
leading-order interactions of the Higgs boson with the gauge fields.
${\cal{L}}_{\text{ferm}}$ parametrises the fermion-gauge boson interactions,
which we take SM-like in the following. $V({H})$ is the Higgs potential, which
we relate to the SM expectation

$V({H})=\frac{1}{2}M_{H}^{2}H^{2}+\kappa_{3}H^{3}+\kappa_{4}H^{4}\,.$ (7)

with $\kappa_{3}\simeq 32~{}\text{GeV},\kappa_{4}\simeq 0.03$ in the SM.

In this work, we consider ALP interactions that particularly probe the singlet
character of the Higgs boson as a parametrised by the HEFT Lagrangian. The
interactions are given by

$\mathcal{L}^{\text{ALP}}_{\text{LO}}=\frac{1}{2}\partial_{\mu}\mathcal{A}\partial^{\mu}\mathcal{A}-\frac{1}{2}M_{\mathcal{A}}^{2}\mathcal{A}^{2}\\\
+a_{2D}\Big{(}i\,v^{2}\text{Tr}[U\tau^{3}U^{\dagger}\,\boldsymbol{\mathcal{V}}_{\mu}]\;\frac{\partial_{\mu}\mathcal{A}}{f_{A}}{\mathcal{F}}_{2D}\Big{)}\,$
(8a) with
${\mathcal{F}}_{2D}=\Big{(}1+2\zeta_{12D}\frac{H}{v}+\zeta_{22D}\Big{(}\frac{H}{v}\Big{)}^{2}+...\Big{)}\,,$
(8b) and $\boldsymbol{\mathcal{V}}_{\mu}=(D_{\mu}U)U^{\dagger}\,.$ (8c)

In Eq. (8), $f_{\mathcal{A}}$ denotes the scale linked with the ALP
interactions. The interactions specified by Eq. (8) are the leading order
chiral interactions of the ALP field with SM states. These couplings
specifically probe the custodial singlet nature of the Higgs boson Brivio et
al. (2017). Therefore, the phenomenology of these interactions provides
relevant insights into the mechanism of electroweak symmetry breaking and its
relation to axion-like states. The radiative imprints of these interactions on
SM correlations are then captured by the chiral dimension-4 interactions
contributing to $\mathcal{L}^{\text{HEFT}}_{\text{LO}}$ when HEFT parameters
coincide with the SM expectation.222All interactions detailed above are
implemented using the FeynRules package Christensen and Duhr (2009); Alloul et
al. (2014).

The aim of our analysis is to clarify the phenomenological reach to the
couplings involved in Eq. (8) from two different angles. Firstly, these
interactions are clear indicators of a singlet character of the Higgs boson in
HEFT. Secondly, the interactions $\sim\zeta_{12D}$ will introduce
modifications to the Higgs boson propagation and Higgs decay in HEFT,
$\zeta_{22D}$ will imply modifications to the Higgs pair production rate.
Although the ALP might be too light to be directly accessible at collider
experiments such as the LHC, its virtual imprint through specific predictions
between the correlations of four-top and Higgs pair production could reveal
its presence. We will turn to the expected constraints in the next section.
Throughout, we will identify the HEFT parameters with their corresponding
tree-level SM limit except for the deviations introduced by the ALP, which we
also detail below. We will focus on the interactions that are generated at
order $a_{2D}\zeta_{12D}$ etc.; fits against the ALP-less HEFT (or the SM as a
particular HEFT parameter choice) should be sensitive to these contributions
when data is consistent with the latter expectation. To reflect this we will
therefore also assume that HEFT operators coincide with their SM expectation.
Specifically this means that we will choose vanishing HEFT parameters arising
at chiral dimension-4. Departures from the SM correlations are then directly
related to (radiative) presence of the ALP.

##  III Direct constraints from Higgs decays

The interactions of an ALP with a Higgs boson via the $a_{2D}\zeta_{12D}$
coupling of Eq. (8) is tree-level mediated. The exotic decay of the Higgs
boson via $H\rightarrow\mathcal{A}\,Z$ at leading order is given by

$\parbox{76.82234pt}{\includegraphics[width=76.82234pt]{HaZ.pdf}}=-2i\frac{e}{c_{W}s_{W}}\frac{v}{f_{\mathcal{A}}}a_{2D}\,\zeta_{12D}\,q^{\mu}(H)\,,$
(9)

with $q^{\mu}(H)$ denoting the four-momentum carried by the Higgs leg. When
kinematically accessible, the decay width of the Higgs boson receives a non-SM
contribution

$\Gamma(H\rightarrow\mathcal{A}Z)=\frac{v^{2}a_{2D}^{2}\zeta_{12D}^{2}}{274\;s_{W}^{2}c_{W}^{2}f_{\mathcal{A}}^{2}M_{H}^{3}M_{Z}^{2}}\\\
\Big{(}M_{\mathcal{A}}^{4}-2M_{\mathcal{A}}^{2}(M_{H}^{2}+M_{Z}^{2})+(M_{H}^{2}-M_{Z}^{2})^{2}\Big{)}^{3/2}.$
(10)

Assuming this two-body process as the most dominating BSM decay involving the
ALP, the SM Higgs boson signal strengths get uniformly modified

Figure 1:  The allowed parameter region obtained from the di-photon signal
strength reported by ATLAS for their $139~{}\text{fb}^{-1}$ data set for
different values of $\zeta_{12D}$.

$\mu_{\text{SM},{\cal{A}}}=\frac{\Gamma(H)_{\text{SM}}}{\Gamma(H\rightarrow\mathcal{A}\,Z)+\Gamma(H)_{\text{SM}}}\,,$
(11)

with $\Gamma(H)_{\text{SM}}\simeq 4~{}\text{MeV}$ as the total Higgs boson
decay width in the SM Dittmaier et al. (2011). To constrain this BSM decay, we
use the well constrained and hence representative signal strength for
$H\rightarrow\gamma\gamma$. This has been measured

$\mu_{\gamma\gamma}=1.04^{+0.1}_{-0.09}~{}\hbox{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{ATLAS:2022tnm}{\@@citephrase{(}}{\@@citephrase{)}}} }$
(12)

for the representative ATLAS Run 2 dataset of $139~{}\text{fb}^{-1}$. For the
on-shell decay of $H\to\mathcal{A}\,Z$, the maximum value of ALP mass allowed
kinematically is $\simeq 34~{}\text{GeV}$ with $M_{H}=125~{}\text{GeV}$ and
$M_{Z}=91.18~{}\text{GeV}$. For heavier ALP masses, the branching ratio
quickly dies off due to the offshellness of the involved $Z$ boson.

The allowed parameter space in $a_{2D}/f_{\mathcal{A}}$ vs $M_{\mathcal{A}}$
plane is shown in Fig. 1 for three different values of $\zeta_{12D}$. The
above $95\%$ limit translates into the lower bound on
$\Gamma(H\rightarrow\mathcal{A}\,Z)<0.65~{}\text{MeV}$ using Eq. (11) for
$\zeta_{12D}=1$. The above bound on $\Gamma(H\rightarrow\mathcal{A}\,Z)$ is
reduced by half with the HL-LHC projections for $H\rightarrow\gamma\gamma$ at
$3~{}\text{ab}^{-1}$ Aad et al. (2018), i.e. we obtain
$\Gamma(H\rightarrow\mathcal{A}\,Z)<0.32~{}\text{MeV}$ for $\zeta_{12D}=1$.

##  IV Higgs signals of Virtual ALPs

### Propagation vs. on-shell properties: Four-top production

BSM corrections to the Higgs self energy $\Sigma_{H}$ can give rise to an
oblique correction

$\hat{H}=-\frac{M_{H}^{2}}{2}\Sigma_{H}^{\prime\prime}(M_{H}^{2})\;,$ (13)

analogously to the $\hat{W}$, $\hat{Y}$ parameters in the gauge sector, e.g.
Barbieri et al. (2004). Such a correction leads to a Higgs propagator
modification Englert et al. (2019)

$-i\Delta_{H}(q^{2})=\frac{1}{q^{2}-M_{H}^{2}}\left(1+\hat{H}\left(1-{q^{2}\over
M_{H}^{2}}\right)\right)\,,$ (14)

indicating a departure for large momentum transfers at unit pole residue.
Measurements of this parameter have by now been established by ATLAS and CMS
in Refs. Aad et al. (2023b); Sirunyan et al. (2020). The expected upper limit
is

$\hat{H}\leq 0.12\,,$ (15)

at $95$% CL from the recent four-top production results of Ref. Aad et al.
(2023b). We can re-interpret this in the framework that we consider. In
parallel, we can employ an extrapolation of four-top final states to estimate
sensitivity improvements that should become available in $\hat{H}$-specific
analyses at the high-luminosity LHC (ATLAS currently observe a small tension
in their $\hat{H}$ fit).

Explicit calculation in general $R_{\xi}$ gauge of the ALP insertion of Eq.
(9) into the Higgs two-point function yields the $\xi$-independent result (see
also remarks in Herrero and Morales (2020, 2021); Anisha et al. (2022))

$\Gamma(H(q)H(q))=\frac{a_{2D}^{2}\,\zeta_{12D}^{2}}{4\pi^{2}f_{\mathcal{A}}^{2}}(4M_{\mathcal{A}}^{4}-3M_{\mathcal{A}}^{2}\,q^{2}\\\
+(q^{2}-3M_{Z}^{2})q^{2})\Delta_{\text{UV}}+\dots\,,$ (16)

with $\overline{\text{MS}}$ factor

$\Delta_{\text{UV}}={\Gamma(1+\epsilon)\over\epsilon}\left({4\pi\mu^{2}\over
M_{H}^{2}}\right)^{\epsilon}$ (17)

in dimensional regularisation $D=4-2\epsilon$ with ‘t Hooft mass $\mu$. The
ellipses in Eq. (16) denote finite terms for $\epsilon\to 0$ (see below).

In the following we will adopt the on-shell scheme for field and mass
renormalisation (cf. Eq. (14)), and the $\overline{\text{MS}}$ scheme for HEFT
parameters (see also Delgado et al. (2014); Herrero and Morales (2020, 2021);
Asiáin et al. (2022); Anisha et al. (2022); Herrero and Morales (2022)). On
the one hand, part of the divergence of Eq. (16) are then cancelled by the
(divergent, div.) counterterms related to the Higgs wave function and mass
renormalisation

$\displaystyle\delta Z_{H}\big{|}_{\text{div.}}$
$\displaystyle=\frac{3\,a_{2D}^{2}\,\zeta_{12D}^{2}(M_{\mathcal{A}}^{2}+M_{Z}^{2})}{4\pi^{2}f_{\mathcal{A}}^{2}}\,,$
(18) $\displaystyle\delta M_{H}^{2}\big{|}_{\text{div.}}$
$\displaystyle=\frac{a_{2D}^{2}\,\zeta_{12D}^{2}(4M_{\mathcal{A}}^{4}-3M_{\mathcal{A}}^{2}M_{H}^{2}-3M_{H}^{2}M_{Z}^{2})}{4\pi^{2}f_{\mathcal{A}}^{2}}\,.$

On the other hand, the appearance of a $~{}q^{4}$ contribution signifies the
sourcing of the chiral dimension-4 operator $\mathcal{O}_{\square\square}$ of
the HEFT Lagrangian

$\mathcal{O}_{\square\square}=a_{\square\square}\frac{\square H\square
H}{v^{2}}\,.$ (19)

This operator is renormalised by the ALP interactions via

$\delta
a_{\square\square}=-\frac{a_{2D}^{2}\zeta_{12D}^{2}v^{2}}{8\pi^{2}f_{\mathcal{A}}^{2}}\Delta_{\text{UV}}\,.$
(20)

Together, the renormalised Higgs two-point function then links to the
$\hat{H}$ parameter as

$\begin{split}\hat{H}=-\frac{a_{2D}^{2}M_{H}^{2}\zeta_{12D}^{2}}{8\pi^{2}f_{\cal{A}}^{2}}\bigg{(}2B_{0}(M_{H}^{2},M_{\cal{A}}^{2},M_{Z}^{2})\big{|}_{\text{fin.}}-4(M_{\cal{A}}^{2}-M_{H}^{2}+M_{Z}^{2})B^{\prime}_{0}(M_{H}^{2},M_{\cal{A}}^{2},M_{Z}^{2})\bigg{.}\\\
\bigg{.}+\left[M_{\cal{A}}^{4}-2M_{\cal{A}}^{2}(M_{H}^{2}+M_{Z}^{2})+(M_{H}^{2}-M_{Z}^{2})^{2}\right]B^{\prime\prime}_{0}(M_{H}^{2},M_{\cal{A}}^{2},M_{Z}^{2})\bigg{)}\,,\end{split}$
(21)

where ‘fin.’ denotes the UV finite part of the Passarino-Veltman $B_{0}$
function after subtracting Eq. (17) and derivatives are taken with respect to
the first argument of the $B_{0}$ function (an explicit representation can be
found in Ref. Denner (1993)). $\hat{H}$ vanishes in the decoupling limit
$f_{\cal{A}}>M_{\cal{A}}\gg M_{H}$.

Equation (21) shows that propagator corrections that can be attributed to
$\hat{H}$ probe similar couplings as the Higgs decay of Eq. (10), however, in
a momentum transfer-enhanced way, at the price of a loop suppression. This way
the energy coverage of the LHC that becomes under increasing statistical
control provides additional sensitivity beyond the fixed scale Higgs decay.
Any enhanced sensitivity to the on vs. off-shell phenomenology that can be
gained from the combination of the processes discussed so far, can then break
the degeneracies between the different HEFT coefficients in Eq. (8b).

To obtain an extrapolation estimate from the current constraints on $\hat{H}$,
we implement the modifications from $\hat{H}$ in MadGraph5_aMC@NLO Alwall et
al. (2014) in order to estimate the changes caused in the four-top cross
section from different contributions to the Higgs self-energy, and extrapolate
the result of Eq. (15). Assuming a significance $S(\hat{H}=0.12)/\sqrt{B}=2$
from the constraint of Ref. Aad et al. (2023b) at $140$/fb, and then
subsequently rescaling the results to $3$/ab, we obtain the approximate
significance at HL-LHC. While using the more recent results yields improved
bounds compared to earlier projections of Ref. Sirunyan et al. (2018) that
include systematics (due to improvements in the analysis procedure utilising
ML techniques), our projections remain conservative compared to the previously
estimated significance with only statistical uncertainties, see Fig. 2.

In Fig. 3, we also see that if $M_{\cal{A}}$ is light, it will freely
propagate in the 2 point function thus imparting the characteristic $q^{4}$
dependence probed by $\hat{H}$. This also means that this behaviour is
essentially independent of the light ALP mass scale. Turning to heavier
states, this kinematic dependence is not sourced as efficiently anymore,
leading to a quick decoupling from the two-point Higgs function and reduced
sensitivity and larger theoretical uncertainty. We will return to the
relevance of $\hat{H}$ for the discussed scenario after discussing the
modifications to Higgs pair production in the next section.

Figure 2: Expected significance for different values of $\hat{H}$ at HL-LHC
using the projections of Ref. Sirunyan et al. (2018) (blue) and our estimates
from projections using Eq. (15) (red).  Figure 3: The solid line indicates the
$2\sigma$ contour on the $a_{2D}/f_{\cal{A}}$-$M_{\cal{A}}$ plane from the
$\hat{H}$ analysis for a scale $\mu=2M_{H}$. The shaded region indicates the
open $H\rightarrow{\cal{A}}Z$ decay channel.

### Higher terms of the ALP flare function: Higgs pair production

Corrections to Higgs pair production under the same assumptions as in the
previous section are contained in propagator corrections and corrections to
trilinear Higgs coupling. As with the chiral dimension-4 operator that leads
to new contributions to the Higgs-two point function, there are additional
operators that modify the Higgs trilinear interactions. The amputated off-
shell three-point function receives contributions (see also Herrero and
Morales (2022))

$v^{3}\Gamma_{1}(H(q)H(k_{1})H(k_{2}))=a_{\chi
1}(q^{4}+k_{1}^{4}+k_{2}^{4})\\\ +2a_{\chi
2}(q^{2}k_{1}^{2}+k_{1}^{2}k_{2}^{2}+q^{2}k_{2}^{2})+a_{\chi
3}v^{2}(q^{2}+k_{1}^{2}+k_{2}^{2})\,,$ (22)

which are renormalised in the $\overline{\text{MS}}$ scheme according to

$\begin{split}\delta a_{\chi
1}&=\frac{a_{2D}^{2}\zeta_{12D}v^{2}}{8\pi^{2}f_{\mathcal{A}}^{2}}(3(1+\zeta_{1})\zeta_{12D}+2\zeta_{22D})\Delta_{\text{UV}}\,,\\\
\delta a_{\chi
2}&=\frac{3a_{2D}^{2}\zeta_{12D}^{2}v^{2}}{8\pi^{2}f_{\mathcal{A}}^{2}}(1+\zeta_{1})\Delta_{\text{UV}}\,,\\\
\delta a_{\chi
3}&=\frac{3a_{2D}^{2}\zeta_{12D}v^{2}}{4\pi^{2}f_{\mathcal{A}}^{2}}\big{[}(M_{\cal{A}}^{2}+M_{Z}^{2})\zeta_{22D}\\\
&\hskip
34.14322pt-3(M_{\cal{A}}^{2}+2M_{Z}^{2})(1+\zeta_{1})\zeta_{12D}\big{]}\Delta_{\text{UV}}\,.\end{split}$
(23)

The remaining renormalisation of the chiral dimension-2 term follows from Eq.
(7)

$\delta\Gamma_{2}(H(q)H(k_{1})H(k_{2}))\big{|}_{\text{div}}=-\frac{9a_{2D}^{2}\zeta_{12D}^{2}}{8\pi^{2}f_{\mathcal{A}}^{2}}\kappa_{3}(M_{\cal{A}}^{2}+m_{Z}^{2})\\\
-\frac{a_{2D}^{2}\zeta_{12D}M_{\cal{A}}^{4}}{2\pi^{2}f_{\mathcal{A}}^{2}v}\left(2(1+\zeta_{1})\zeta_{12D}-\zeta_{22D}\right)\,.$
(24)

(a) (b)

Figure 4:  Invariant mass of the Higgs pair $M^{2}_{HH}=(p_{H,1}+p_{H,2})^{2}$
for different values of $\zeta_{12D,22D}$ and the SM in $gg\to HH$, and
$M_{\cal{A}}=1~{}\text{GeV}$ is shown on the left. Considering the final state
Higgs boson on-shell, the corrections to their decays according to Eq. (10)
are not included. The normalisation is chose to $\sigma(HH)^{\text{SM}}\simeq
32~{}\text{fb}$ de Florian et al. (2016) and we use SM K factors to
qualitatively include higher order QCD corrections. Squared contributions from
the renormalised $gg\to hh$ amplitude are included. On the right the ratio of
the new-physics cross section with respect to the SM is shown for different
values in the $\zeta_{12D}$-$\zeta_{22D}$ plane.

(a) (b)

Figure 5:  $95$% CL regions for the case of a light
$M_{\cal{A}}=1~{}\text{GeV}$ (left) and heavy $M_{\cal{A}}=40~{}\text{GeV}$
(right) ALP.

ATLAS (CMS) have set highly competitive expected 95% confidence level cross
section limits of $\sigma/\sigma_{\text{SM}}<3.9~{}(5.2)$ Tumasyan et al.
(2023); Aad et al. (2022b) in the $b\bar{b}\tau\tau$ channel Dolan et al.
(2012) alone. Slightly reduced sensitivity Aad et al. (2022c); Tumasyan et al.
(2022); Sirunyan et al. (2021) can be achieved in the $4b$ and $2b2\gamma$
modes Baur et al. (2004); Ferreira de Lima et al. (2014). ATLAS have combined
these channels to obtain a combined exclusion of $3.1\sigma_{\text{SM}}$ Aad
et al. (2021b) with the currently available data and forecast a sensitivity of
$\sigma/\sigma_{\text{SM}}\gtrsim 1.1$ at the HL-LHC Aad et al. (2022d). We
use the two latter result to gain a qualitative sensitivity reach of Higgs
pair production in the considered scenario.

In Fig. 4, we show representative invariant Higgs pair mass distributions for
13 TeV LHC collisions, which demonstrates the potential of multi-Higgs final
states’ sensitivity to the momentum-enhanced new physics contributions
characteristic to the ALP.333We have implemented these changes into an in-
house Monte Carlo event generate based on Vbfnlo Arnold et al. (2009); Baglio
et al. (2011) employing FeynArts, FormCalc, and LoopTools Hahn and Perez-
Victoria (1999); Hahn (2000, 2001); Hahn and Rauch (2006) and PackageX Patel
(2017) for numerical and analytical cross checks. Throughout this work we
chose a renormalisation scale of $\mu=2M_{h}$. The behaviour exhibited by the
invariant mass distribution is not sensitive to the mass of the ALP as long as
the latter is not close to the $\simeq 2M_{H}$ threshold that determines the
$gg\to HH$ phenomenology. In instances when hefty ALPs propagate freely, their
distinctive momentum enhancements will sculpt the Higgs-boson distributions.
In parallel, non-linear effects will be important away from the SM reference
point as shown in Fig. 4. This shows that the constraints that can be obtained
in the di-Higgs channel are relatively strongly coupled, which is motivation
for us to directly include “squared” BSM effects to our analysis in addition
to interference effects.

We combine the three representative analyses in a global $\chi^{2}$ to obtain
sensitivity estimates. In the case when the ALP is light, there are
significant modifications to Higgs physics, also at large momentum transfers,
see also Fig. 4. Of course, these large contributions in particular to the
Higgs pair rate are tamed by decreasing signal strengths into SM-like states,
which quickly result in tension with experimental observations for larger
couplings. As Higgs pair production observations need to rely on relatively
clean and high branching ratio final states, the prospects of Higgs pair
production (and four top) analyses to provide additional sensitivity is
relatively low. This is highlighted already in the combination of the Higgs
decay constraints with these processes in Fig. 5(a).

For parameter choices for which the ALP is above the Higgs decay threshold,
this picture changes. Multi-Higgs constraints remain relatively insensitive to
the ALP mass scale as long as these states are away from the $2M_{H}$
threshold. The cross section enhancement then translates directly into an
enhancement of the observable Higgs boson pair production rate. In turn,
constraints on the the higher order terms in the ALP flare function become
possible. It is important to note that these are independent of couplings (to
first order) that shape the ALP decay phenomenology.

As the large enhancements result from the tails of distributions there is a
question of validity. Nonetheless the momentum dependence introduced by Eq.
(9) leads to partial wave unitarity violation as, e.g. $HZ$ scattering
proceeds momentum-enhanced. A numerical investigation shows that for
${\cal{O}}(1)$ couplings in Eq. (8), conserved zeroth partial wave unitarity
up to scales $\sim 1.5~{}\text{TeV}$ sets a lower bound of $f_{a}\gtrsim
300~{}\text{GeV}$ for unsuppressed propagation ${M}_{\cal{A}}=1~{}\text{GeV}$.
These constraints are driven by the longitudinal $Z$ polarisations,
constraints from transverse modes are comparably weaker. This means that the
entire region that is shown in Fig. 4(a) is perturbative at tree-level. In
parallel, the HL-LHC is unlikely to probe Higgs pairs beyond invariant masses
$M_{HH}>600~{}\text{GeV}$ in the SM (for which the cross section drops to 10%
of the inclusive rate). Most sensitivity in HL-LHC searches results from the
threshold region. Therefore, the sensitivity expected by the HL extrapolation
of Aad et al. (2022d) will probe Eq. (8) in a perturbatively meaningful
regime.

The combined constraints are largely driven by Higgs pair constraints, Fig.
5(b). However, it is worth highlighting that the statistics-only extrapolation
does not include changes to the four top search methodology. Improvements of
the latter can be expected with increasing luminosity and the final verdict
from four top production might indeed be much more optimistic than our
$\sqrt{\text{luminosity}}$ extrapolation might suggest.

##  V Summary and Conclusions

Searches for new light propagating degrees of freedom such as axion-like
particles are cornerstones of the BSM programme in particle physics as
explored at, e.g., the Large Hadron Collider. The Higgs boson, since a global
picture of its interactions is still incomplete, provides a motivated avenue
for the potential discovery of new physics in the near future as the LHC
experiments gain increasingly phenomenological sensitivity in rare processes
that could be tell-tale signs of Higgs-related BSM physics.

We take recent experimental developments in multi-Higgs and multi-top analyses
as motivation to analyse effective Higgs-philic ALP interactions, also beyond
leading order. This enables us to tension constraints from different areas of
precision Higgs phenomenology, combining Higgs decay modifications with large-
momentum transfer processes that are becoming increasingly accessible at the
LHC. For light states and sizeable HEFT-like couplings, a large part of the
sensitivity is contained in Higgs signal strength measurements (see also
Brivio et al. (2017)), which, however, only provide limited insights into the
Higgs-ALP interactions. Higher terms of the Higgs-ALP flare function, still
have the phenomenological potential to sizeably modify Higgs pair final states
at a level that will be observable at the LHC in the near future. Our findings
therefore also highlight further the relevance of multi-top and multi-Higgs
final state for the quest for new physics.

###### Acknowledgements.

We thank Dave Sutherland for insightful discussions. C.E. thanks the high-
energy physics group (FTAE) at the University of Granada for their hospitality
during early stages of this work. A. is supported by the Leverhulme Trust
under grant RPG-2021-031. S.D.B is supported by SRA (Spain) under Grant No.
PID2019-106087GB-C21 (10.13039/501100011033), and
PID2021-128396NB-100/AEI/10.13039/501100011033; by the Junta de Andalucía
(Spain) under Grants No. FQM-101, A-FQM-467-UGR18, and P18-FR-4314 (FEDER).
C.E. is supported by the STFC under grant ST/T000945/1, the Leverhulme Trust
under grant RPG-2021-031, and the Institute for Particle Physics Phenomenology
Associateship Scheme. P.S. is supported by the Deutsche Forschungsgemeinschaft
under Germany’s Excellence strategy EXC2121 “Quantum Universe” - 390833306.
This work has been partially funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) - 491245950.

## References

  * Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Graham et al. (2015) P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, Ann. Rev. Nucl. Part. Sci. 65, 485 (2015), eprint 1602.00039.
  * Irastorza and Redondo (2018) I. G. Irastorza and J. Redondo, Prog. Part. Nucl. Phys. 102, 89 (2018), eprint 1801.08127.
  * Jaeckel and Spannowsky (2016) J. Jaeckel and M. Spannowsky, Phys. Lett. B 753, 482 (2016), eprint 1509.00476.
  * Esser et al. (2023) F. Esser, M. Madigan, V. Sanz, and M. Ubiali (2023), eprint 2303.17634.
  * Alonso-Álvarez et al. (2023) G. Alonso-Álvarez, J. Jaeckel, and D. D. Lopes (2023), eprint 2302.12262.
  * Bauer et al. (2021) M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, Phys. Rev. Lett. 127, 081803 (2021), eprint 2102.13112.
  * Bauer et al. (2022) M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, JHEP 09, 056 (2022), eprint 2110.10698.
  * Bauer et al. (2017) M. Bauer, M. Neubert, and A. Thamm, JHEP 12, 044 (2017), eprint 1708.00443.
  * Biekötter et al. (2022) A. Biekötter, M. Chala, and M. Spannowsky, Phys. Lett. B 834, 137465 (2022), eprint 2203.14984.
  * Bauer et al. (2023) M. Bauer, G. Rostagni, and J. Spinner, Phys. Rev. D 107, 015007 (2023), eprint 2207.05762.
  * Bonilla et al. (2021) J. Bonilla, I. Brivio, M. B. Gavela, and V. Sanz, JHEP 11, 168 (2021), eprint 2107.11392.
  * Buchalla et al. (2014) G. Buchalla, O. Catà, and C. Krause, Nucl. Phys. B 880, 552 (2014), [Erratum: Nucl.Phys.B 913, 475–478 (2016)], eprint 1307.5017.
  * Buchalla et al. (2018) G. Buchalla, O. Cata, A. Celis, M. Knecht, and C. Krause, Nucl. Phys. B 928, 93 (2018), eprint 1710.06412.
  * Buchalla et al. (2021) G. Buchalla, O. Catà, A. Celis, M. Knecht, and C. Krause, Phys. Rev. D 104, 076005 (2021), eprint 2004.11348.
  * ACME Collaboration (2018) ACME Collaboration, Nature 562, 355—360 (2018), ISSN 0028-0836, URL https://doi.org/10.1038/s41586-018-0599-8.
  * Abel et al. (2020) C. Abel, S. Afach, N. J. Ayres, C. A. Baker, G. Ban, G. Bison, K. Bodek, V. Bondar, M. Burghoff, E. Chanel, et al., Phys. Rev. Lett. 124, 081803 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.081803.
  * Raffelt (1996) G. G. Raffelt, _Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles_ (1996), ISBN 978-0-226-70272-8.
  * Raffelt and Stodolsky (1988) G. Raffelt and L. Stodolsky, Phys. Rev. D 37, 1237 (1988), URL https://link.aps.org/doi/10.1103/PhysRevD.37.1237.
  * O’Hare (2020) C. O’Hare, _cajohare/axionlimits: Axionlimits_ , https://cajohare.github.io/AxionLimits/ (2020).
  * Aad et al. (2021a) G. Aad et al. (ATLAS), JHEP 03, 243 (2021a), [Erratum: JHEP 11, 050 (2021)], eprint 2008.05355.
  * Sirunyan et al. (2019) A. M. Sirunyan et al. (CMS), Phys. Lett. B 797, 134826 (2019), eprint 1810.04602.
  * d’Enterria (2021) D. d’Enterria, in _Workshop on Feebly Interacting Particles_ (2021), eprint 2102.08971.
  * Sirunyan et al. (2023) A. M. Sirunyan et al. (CMS), Tech. Rep., CERN, Geneva (2023), URL http://cds.cern.ch/record/2853524.
  * 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.
  * Tumasyan et al. (2023) A. Tumasyan et al. (CMS), Phys. Lett. B 842, 137531 (2023), eprint 2206.09401.
  * Aad et al. (2023a) G. Aad et al. (ATLAS) (2023a), eprint 2301.03212.
  * Hayrapetyan et al. (2023) A. Hayrapetyan et al. (CMS) (2023), eprint 2305.13439.
  * Aad et al. (2023b) G. Aad et al. (ATLAS) (2023b), eprint 2303.15061.
  * Englert et al. (2019) C. Englert, G. F. Giudice, A. Greljo, and M. Mccullough, JHEP 09, 041 (2019), eprint 1903.07725.
  * Sirunyan et al. (2020) A. M. Sirunyan et al. (CMS), Eur. Phys. J. C 80, 75 (2020), eprint 1908.06463.
  * Alonso et al. (2013) R. Alonso, M. B. Gavela, L. Merlo, S. Rigolin, and J. Yepes, Phys. Lett. B 722, 330 (2013), [Erratum: Phys.Lett.B 726, 926 (2013)], eprint 1212.3305.
  * Brivio et al. (2016) I. Brivio, J. Gonzalez-Fraile, M. C. Gonzalez-Garcia, and L. Merlo, Eur. Phys. J. C 76, 416 (2016), eprint 1604.06801.
  * Herrero and Morales (2020) M. Herrero and R. A. Morales, Phys. Rev. D 102, 075040 (2020), eprint 2005.03537.
  * Herrero and Morales (2021) M. J. Herrero and R. A. Morales, Phys. Rev. D 104, 075013 (2021), eprint 2107.07890.
  * Anisha et al. (2022) Anisha, O. Atkinson, A. Bhardwaj, C. Englert, and P. Stylianou, JHEP 10, 172 (2022), eprint 2208.09334.
  * Christensen and Duhr (2009) N. D. Christensen and C. Duhr, Comput. Phys. Commun. 180, 1614 (2009), eprint 0806.4194.
  * 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.
  * Dittmaier et al. (2011) S. Dittmaier et al. (LHC Higgs Cross Section Working Group) (2011), eprint 1101.0593.
  * Aad et al. (2022a) G. Aad et al. (ATLAS) (2022a), eprint 2207.00348.
  * Aad et al. (2018) G. Aad et al. (ATLAS) (2018), eprint ATL-PHYS-PUB-2018-054.
  * Barbieri et al. (2004) R. Barbieri, A. Pomarol, R. Rattazzi, and A. Strumia, Nucl. Phys. B 703, 127 (2004), eprint hep-ph/0405040.
  * Delgado et al. (2014) R. L. Delgado, A. Dobado, and F. J. Llanes-Estrada, JHEP 02, 121 (2014), eprint 1311.5993.
  * Asiáin et al. (2022) I. n. Asiáin, D. Espriu, and F. Mescia, Phys. Rev. D 105, 015009 (2022), eprint 2109.02673.
  * Herrero and Morales (2022) M. J. Herrero and R. A. Morales, Phys. Rev. D 106, 073008 (2022), eprint 2208.05900.
  * Denner (1993) A. Denner, Fortsch. Phys. 41, 307 (1993), eprint 0709.1075.
  * 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.
  * Sirunyan et al. (2018) A. M. Sirunyan et al. (CMS) (2018), eprint CMS-PAS-FTR-18-031.
  * de Florian et al. (2016) D. de Florian et al. (LHC Higgs Cross Section Working Group), 2/2017 (2016), eprint 1610.07922.
  * Aad et al. (2022b) G. Aad et al. (ATLAS) (2022b), eprint 2209.10910.
  * Dolan et al. (2012) M. J. Dolan, C. Englert, and M. Spannowsky, JHEP 10, 112 (2012), eprint 1206.5001.
  * Aad et al. (2022c) G. Aad et al. (ATLAS), Phys. Rev. D 106, 052001 (2022c), eprint 2112.11876.
  * Tumasyan et al. (2022) A. Tumasyan et al. (CMS), Phys. Rev. Lett. 129, 081802 (2022), eprint 2202.09617.
  * Sirunyan et al. (2021) A. M. Sirunyan et al. (CMS), JHEP 03, 257 (2021), eprint 2011.12373.
  * Baur et al. (2004) U. Baur, T. Plehn, and D. L. Rainwater, Phys. Rev. D 69, 053004 (2004), eprint hep-ph/0310056.
  * Ferreira de Lima et al. (2014) D. E. Ferreira de Lima, A. Papaefstathiou, and M. Spannowsky, JHEP 08, 030 (2014), eprint 1404.7139.
  * Aad et al. (2021b) G. Aad et al. (ATLAS) (2021b), eprint ATLAS-CONF-2021-052.
  * Aad et al. (2022d) G. Aad et al. (ATLAS) (2022d), eprint ATL-PHYS-PUB-2022-053.
  * Arnold et al. (2009) K. Arnold et al., Comput. Phys. Commun. 180, 1661 (2009), eprint 0811.4559.
  * Baglio et al. (2011) J. Baglio et al. (2011), eprint 1107.4038.
  * Hahn and Perez-Victoria (1999) T. Hahn and M. Perez-Victoria, Comput. Phys. Commun. 118, 153 (1999), eprint hep-ph/9807565.
  * Hahn (2000) T. Hahn, Nucl. Phys. B Proc. Suppl. 89, 231 (2000), eprint hep-ph/0005029.
  * Hahn (2001) T. Hahn, Comput. Phys. Commun. 140, 418 (2001), eprint hep-ph/0012260.
  * Hahn and Rauch (2006) T. Hahn and M. Rauch, Nucl. Phys. B Proc. Suppl. 157, 236 (2006), eprint hep-ph/0601248.
  * Patel (2017) H. H. Patel, Comput. Phys. Commun. 218, 66 (2017), eprint 1612.00009.