Discovery prospects for photophobic axion-like particles in the WWjj final state at the High-Luminosity LHC

Quantum chromodynamics (QCD) allows for a CP-violating θ term, yet no corresponding effect has been observed in measurements of the neutron electric dipole moment [1, 2]. This long–standing "strong CP problem" can be resolved by enlarging the Standard Model (SM) with a global $U(1)_{\rm PQ}$ symmetry, known as the Peccei-Quinn symmetry, whose spontaneous breaking generates a pseudo-Nambu–Goldstone boson, the axion [3−6]. In traditional QCD-axion models, the axion mass and its couplings are tightly linked to the symmetry-breaking scale [5−12].

However, in many extensions of the SM, pseudo-scalar states with axion-like couplings appear, whose masses and interaction strengths are not fixed by the strong CP solution. These states are collectively referred to as axion-like particles (ALPs) [13−16]. The decoupling of mass and couplings opens a large parameter space and has led to extensive search programs, ranging from astrophysical and cosmological probes to intensity-frontier experiments and high-energy colliders [17−68].

At colliders, ALPs can interact with SM gauge bosons and/or fermions and can be probed either as on-shell resonances or through off-shell exchange in high-energy processes. A large body of experimental work has already constrained the ALP coupling to photons, $ g_{a\gamma\gamma} $, using $ e^+e^- $ data and heavy-ion or proton–proton collisions, particularly via diphoton or light-by-light signatures [69−76]. The resulting bounds are strong over a wide mass range and motivate scenarios in which $ g_{a\gamma\gamma} $ is tuned to be small while couplings to the electroweak gauge bosons remain sizeable. In such photophobic constructions, the $ a\gamma\gamma $ vertex is suppressed, but interactions with $ W^\pm $, Z, and $ \gamma Z $ can persist [77−79], leading to collider signatures that differ qualitatively from generic diphoton ALP searches.

Heavy photophobic ALPs at the LHC have been explored in several ways. One line of work has reinterpreted SM triboson measurements at $ \sqrt{s}=8\; \,{{\rm{TeV}}} $ as signals of ALP-mediated processes, deriving limits up to a few hundred GeV in mass [77]. Other studies have mapped nonresonant ALP exchange onto vector-boson-scattering topologies [80] or recast Run-2 analyses targeting $ W^\pm W^\pm W^\mp $ and $ Z\gamma $ final states to constrain resonant ALP production and decay into electroweak boson pairs [81]. An overview of these and related collider bounds can be found in Ref. [82]. Because the original analyses were optimized for SM processes rather than for on-shell ALP production, differences in acceptance, kinematic selections, and background control often make the recast limits conservative and non-uniform across the parameter space.

More recently, dedicated projections have been carried out specifically for photophobic ALPs at the High-Luminosity LHC (HL-LHC) [83]. In Ref. [82], a detector-level study of $ pp\to jj\,a(\to Z\gamma) $ was performed, including both s-channel and vector-boson–fusion (VBF)–like production mechanisms, reconstructing the ALP resonance and using a multivariate analysis (MVA) to obtain discovery reaches on both $ g_{aWW} $ and $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to Z\gamma) $ over $ m_a\simeq(100 $–$ 4000)\; \,{{\rm{GeV}}} $ at $ \sqrt{s}=14\; \,{{\rm{TeV}}} $. In parallel, Ref. [84] proposed a tri-W signature, $ pp\to W^\pm X(\to W^+W^-) $, using a same-sign dimuon plus hadronic-W final state and an MVA tailored for background rejection; with heavy photophobic ALPs used as a benchmark, that study obtained sensitivities well beyond those of earlier triboson recasts while remaining applicable to more general neutral resonances coupled to W bosons. These works highlight two important lessons: clean electroweak final states with sharp kinematic handles and multivariate discrimination can dramatically sharpen sensitivity, and photophobic ALP searches should exploit multiple channels with different couplings and systematics.

The present study extends this program to the $ WWjj $ topology. We consider photophobic ALPs that couple to the electroweak gauge bosons and are produced in association with two jets via both s-channel electroweak exchange and VBF-like topologies in $ pp $ collisions, $ pp\to jj\,a $, followed by $ a\to W^+W^- $ decay. Both production mechanisms are simulated and analyzed on an equal footing within a single selection-and-MVA framework. This strategy utilizes the full photophobic electroweak structure ($ aWW $, $ aZZ $, $ a\gamma Z $) and allows the analysis to adapt to whichever kinematic pattern dominates: s-channel–like configurations are more important at lower $ m_a $, while VBF-like events with forward jets and large rapidity gaps become increasingly relevant at high $ m_a $. As a clean and experimentally robust final state, we focus on the different-flavour dilepton mode $ W^+W^-\to e^\pm\mu^\mp\nu\bar{\nu} $ accompanied by at least two jets, which efficiently suppresses Drell–Yan and multijet backgrounds. Our study provides projected discovery sensitivities on both the coupling $ g_{aWW} $ and fiducial quantity $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $ at the HL-LHC with $ \sqrt{s}=14\; \,{{\rm{TeV}}} $ and $ {\cal{L}}=3\; \ {\rm{ab}}^{-1} $.

Throughout this study, unless stated otherwise, we impose the photophobic condition $ g_{a\gamma\gamma}=0 $, include both s-channel and VBF-like production of $ pp\to jj\,a $ with $ a\to W^+W^- $ decay, and optimize the boosted decision tree (BDT) threshold independently for each assumed ALP mass. The quoted significances are purely statistical (systematic uncertainties are not included; cf. Sec. VI).

We note that our HL-LHC projection in the $ pp \to jj\,a(\to W^+W^-) $ topology should be read in the context of several recent and complementary studies, including those on the $ jj\,a(\to Z\gamma) $ channel [82], tri-W strategy $ pp \to W^\pm X(\to W^+W^-) $ [84], and related Run-1/2 recasts [77, 80, 81]. Relative to the tri-W strategy [84], our $ WWjj $ analysis typically achieves a stronger reach across the full mass range, with a particularly pronounced advantage for $ m_a\gtrsim 700\;\ {\rm{GeV}} $, owing to the VBF-like production and kinematics. By contrast, $ jj\,a(\to Z\gamma) $ [82] attains a stronger reach overall—benefiting from a narrow $ Z\gamma $ resonance and cleaner backgrounds. Meanwhile, $ jj\,a(\to W^+W^-) $ remains complementary: it directly probes the $ aWW $ vertex, thus enabling nontrivial cross-checks of the photophobic relations among $ aWW $, $ aZ\gamma $, and $ aZZ $ couplings when combined with previous results [82, 84]; it is subject to largely orthogonal systematics and backgrounds relative to $ Z\gamma $; and our model-independent thresholds in $ \sigma(pp \to jj\,a)\times{\rm{Br}}(a \to W^+W^-) $ facilitate straightforward reinterpretation for other neutral resonances that couple to W bosons [84].

The remainder of this paper is organized as follows. Sec. II introduces the effective description of photophobic ALPs and their electroweak couplings. In Secs. III and IV, we define the signal process and summarize the relevant SM backgrounds. The event selection, observable construction, and multivariate strategy are described in Sec. V. Our main results are presented in Sec. VI, and we conclude in Sec. VII. Additional distributions, BDT response shapes, and efficiency tables are presented in the appendices.

We work with an ALP whose interactions with the SM gauge bosons arise solely through the ${S U}(2)_{L}$ and ${U}(1)_{Y}$ field strengths. At the level before electroweak symmetry breaking (EWSB), the relevant effective Lagrangian can be written as [85]

where a and $ m_a $ denote the ALP field and its mass, respectively, and $ f_a $ is the associated decay constant. The tensors $ W_{\mu\nu}^{b} $ (with $ b=1,2,3 $) and $ B_{\mu\nu} $ are the usual ${S U}(2)_{L}$ and ${U}(1)_{Y}$ field strengths, and the corresponding dual tensor is defined as $ \widetilde X^{\mu\nu}\equiv {1}/{2}\epsilon^{\mu\nu\alpha\beta}X_{\alpha\beta} $ with $ X=W^b,B $. The dimensionless coefficients $ c_{\widetilde{W}} $ and $ c_{\widetilde{B}} $ encode the strengths of the ALP couplings to the electroweak gauge sector.

After EWSB, it is convenient to express the interactions in the mass eigenstate basis for the gauge bosons. The terms involving neutral and charged electroweak bosons can be recast as

with the tree-level relations

Here, θ is the weak mixing angle, and we define $ s_\theta\equiv\sin\theta $, $ c_\theta\equiv\cos\theta $, and $ s_{2\theta}\equiv\sin 2\theta $ for simplicity. Guided by the strong experimental constraints on the diphoton channel, we restrict ourselves to the photophobic case in which the $ a\gamma\gamma $ coupling vanishes, $ g_{a\gamma\gamma}=0 $ (see, e.g., [77−82, 84]). Imposing

with $ t_\theta\equiv\tan\theta $, $ g_{a\gamma\gamma} $ is removed at the tree level; thus, the remaining electroweak couplings correlate with each other. In particular, one finds

In the following, we use $ m_a $ and $ g_{aWW} $ as the primary parameters, with $ g_{aZZ} $ and $ g_{aZ\gamma} $ fixed by these relations throughout our phenomenological analysis of the photophobic scenario.

Deviations from the strict photophobic limit can arise in realistic UV completions: $ g_{a\gamma\gamma}=0 $ is a matching condition and can receive threshold corrections; therefore, a small tree-level $ g_{a\gamma\gamma} $ at collider scales is possible [16, 77]. We parametrize a controlled departure from the strict relation via

for which $ g_{a\gamma\gamma}\propto\delta $ and thus $ {\rm{Br}}(a\to\gamma\gamma)\propto\delta^2 $ (at fixed $ g_{aWW} $). For $ m_a $ above the $ WW $ threshold, $ a\to W^+W^- $ remains dominant for modest departures, while larger δ enhances $ a\to\gamma\gamma $ and dilutes $ {\rm{Br}}(a\to W^+W^-) $ (see Appendix A.1). Independently, within the strict photophobic effective field theory (EFT), the photophobic condition is preserved by renormalization group (RG) evolution up to the two-loop level, and loop/RG-induced rare channels remain at the percent level (or below) in our mass range [82]. Therefore, our results are also reported as $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ to enable straightforward recasts.

We interpret the dimension-5 operators in Eq. (1) as part of an EFT valid below a model-dependent cutoff Λ; EFT-validity, unitarity, and width considerations become particularly relevant in the high-mass/high-coupling region. A weakly-coupled UV completion suggests $ \Lambda \lesssim f_a/|c_{\widetilde{W}}| $ up to coupling/loop factors, while naive dimensional analysis (NDA) for strong dynamics suggests $ \Lambda \lesssim 4\pi f_a/|c_{\widetilde{W}}| $ [86, 87]. Using $ g_{aWW}=4c_{\widetilde{W}}/f_a $, this corresponds to

We interpret our results in the region where $ m_a $ and representative hard scales (e.g., momentum transfer) are below Λ. Partial-wave unitarity provides an independent (process-dependent) guide: in the gauge-invariant ALP EFT the strongest bound arises from $ VV\to Va $ (dominant $ J=1 $ eigenvalue), primarily constraining the SU(2)$ _L $ coupling [88]; under our convention, this leads to Eq. (A2) in Appendix A.2, which we use as a conservative validity indicator.

Although $ m_a $ cannot be fully reconstructed in the $ \ell\ell+E_T^{\rm{miss}} $ final state, our signal modelling assumes on-shell $ pp\to jja $ and we present results in the factorised form $ \sigma(pp\to jja)\times{\rm{Br}}(a\to WW) $. We therefore compute the total width in the strict photophobic benchmark with $\tt{MadGraph5\_aMC@NLO} $ and report $ \Gamma_a/m_a $ along the $ 2\sigma/5\sigma $ reach curves in Appendix A.2. The width remains small in the sensitivity region for $ m_a\lesssim 2\;\ {\rm{TeV}} $, while the extreme high-mass/high-coupling corner becomes broad and is treated as an extrapolation ¹; We therefore quote our primary results also as model-independent thresholds on $ \sigma(pp\to jja)\times{\rm{Br}}(a\to WW) $ for reinterpretation; see Appendix A.2 for the EFT-validity, unitarity, and width checks.

We work with an ALP whose interactions with the SM gauge bosons arise solely through the ${S U}(2)_{L}$ and ${U}(1)_{Y}$ field strengths. At the level before electroweak symmetry breaking (EWSB), the relevant effective Lagrangian can be written as [85]

where a and $ m_a $ denote the ALP field and its mass, respectively, and $ f_a $ is the associated decay constant. The tensors $ W_{\mu\nu}^{b} $ (with $ b=1,2,3 $) and $ B_{\mu\nu} $ are the usual ${S U}(2)_{L}$ and ${U}(1)_{Y}$ field strengths, and the corresponding dual tensor is defined as $ \widetilde X^{\mu\nu}\equiv {1}/{2}\epsilon^{\mu\nu\alpha\beta}X_{\alpha\beta} $ with $ X=W^b,B $. The dimensionless coefficients $ c_{\widetilde{W}} $ and $ c_{\widetilde{B}} $ encode the strengths of the ALP couplings to the electroweak gauge sector.

After EWSB, it is convenient to express the interactions in the mass eigenstate basis for the gauge bosons. The terms involving neutral and charged electroweak bosons can be recast as

with the tree-level relations

Here, θ is the weak mixing angle, and we define $ s_\theta\equiv\sin\theta $, $ c_\theta\equiv\cos\theta $, and $ s_{2\theta}\equiv\sin 2\theta $ for simplicity. Guided by the strong experimental constraints on the diphoton channel, we restrict ourselves to the photophobic case in which the $ a\gamma\gamma $ coupling vanishes, $ g_{a\gamma\gamma}=0 $ (see, e.g., [77−82, 84]). Imposing

with $ t_\theta\equiv\tan\theta $, $ g_{a\gamma\gamma} $ is removed at the tree level; thus, the remaining electroweak couplings correlate with each other. In particular, one finds

In the following, we use $ m_a $ and $ g_{aWW} $ as the primary parameters, with $ g_{aZZ} $ and $ g_{aZ\gamma} $ fixed by these relations throughout our phenomenological analysis of the photophobic scenario.

Deviations from the strict photophobic limit can arise in realistic UV completions: $ g_{a\gamma\gamma}=0 $ is a matching condition and can receive threshold corrections; therefore, a small tree-level $ g_{a\gamma\gamma} $ at collider scales is possible [16, 77]. We parametrize a controlled departure from the strict relation via

for which $ g_{a\gamma\gamma}\propto\delta $ and thus $ {\rm{Br}}(a\to\gamma\gamma)\propto\delta^2 $ (at fixed $ g_{aWW} $). For $ m_a $ above the $ WW $ threshold, $ a\to W^+W^- $ remains dominant for modest departures, while larger δ enhances $ a\to\gamma\gamma $ and dilutes $ {\rm{Br}}(a\to W^+W^-) $ (see Appendix A.1). Independently, within the strict photophobic effective field theory (EFT), the photophobic condition is preserved by renormalization group (RG) evolution up to the two-loop level, and loop/RG-induced rare channels remain at the percent level (or below) in our mass range [82]. Therefore, our results are also reported as $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ to enable straightforward recasts.

We interpret the dimension-5 operators in Eq. (1) as part of an EFT valid below a model-dependent cutoff Λ; EFT-validity, unitarity, and width considerations become particularly relevant in the high-mass/high-coupling region. A weakly-coupled UV completion suggests $ \Lambda \lesssim f_a/|c_{\widetilde{W}}| $ up to coupling/loop factors, while naive dimensional analysis (NDA) for strong dynamics suggests $ \Lambda \lesssim 4\pi f_a/|c_{\widetilde{W}}| $ [86, 87]. Using $ g_{aWW}=4c_{\widetilde{W}}/f_a $, this corresponds to

We interpret our results in the region where $ m_a $ and representative hard scales (e.g., momentum transfer) are below Λ. Partial-wave unitarity provides an independent (process-dependent) guide: in the gauge-invariant ALP EFT the strongest bound arises from $ VV\to Va $ (dominant $ J=1 $ eigenvalue), primarily constraining the SU(2)$ _L $ coupling [88]; under our convention, this leads to Eq. (A2) in Appendix A.2, which we use as a conservative validity indicator.

Although $ m_a $ cannot be fully reconstructed in the $ \ell\ell+E_T^{\rm{miss}} $ final state, our signal modelling assumes on-shell $ pp\to jja $ and we present results in the factorised form $ \sigma(pp\to jja)\times{\rm{Br}}(a\to WW) $. We therefore compute the total width in the strict photophobic benchmark with $\tt{MadGraph5\_aMC@NLO} $ and report $ \Gamma_a/m_a $ along the $ 2\sigma/5\sigma $ reach curves in Appendix A.2. The width remains small in the sensitivity region for $ m_a\lesssim 2\;\ {\rm{TeV}} $, while the extreme high-mass/high-coupling corner becomes broad and is treated as an extrapolation ¹; We therefore quote our primary results also as model-independent thresholds on $ \sigma(pp\to jja)\times{\rm{Br}}(a\to WW) $ for reinterpretation; see Appendix A.2 for the EFT-validity, unitarity, and width checks.

We work with an ALP whose interactions with the SM gauge bosons arise solely through the ${S U}(2)_{L}$ and ${U}(1)_{Y}$ field strengths. At the level before electroweak symmetry breaking (EWSB), the relevant effective Lagrangian can be written as [85]

where a and $ m_a $ denote the ALP field and its mass, respectively, and $ f_a $ is the associated decay constant. The tensors $ W_{\mu\nu}^{b} $ (with $ b=1,2,3 $) and $ B_{\mu\nu} $ are the usual ${S U}(2)_{L}$ and ${U}(1)_{Y}$ field strengths, and the corresponding dual tensor is defined as $ \widetilde X^{\mu\nu}\equiv {1}/{2}\epsilon^{\mu\nu\alpha\beta}X_{\alpha\beta} $ with $ X=W^b,B $. The dimensionless coefficients $ c_{\widetilde{W}} $ and $ c_{\widetilde{B}} $ encode the strengths of the ALP couplings to the electroweak gauge sector.

After EWSB, it is convenient to express the interactions in the mass eigenstate basis for the gauge bosons. The terms involving neutral and charged electroweak bosons can be recast as

with the tree-level relations

Here, θ is the weak mixing angle, and we define $ s_\theta\equiv\sin\theta $, $ c_\theta\equiv\cos\theta $, and $ s_{2\theta}\equiv\sin 2\theta $ for simplicity. Guided by the strong experimental constraints on the diphoton channel, we restrict ourselves to the photophobic case in which the $ a\gamma\gamma $ coupling vanishes, $ g_{a\gamma\gamma}=0 $ (see, e.g., [77−82, 84]). Imposing

with $ t_\theta\equiv\tan\theta $, $ g_{a\gamma\gamma} $ is removed at the tree level; thus, the remaining electroweak couplings correlate with each other. In particular, one finds

In the following, we use $ m_a $ and $ g_{aWW} $ as the primary parameters, with $ g_{aZZ} $ and $ g_{aZ\gamma} $ fixed by these relations throughout our phenomenological analysis of the photophobic scenario.

Deviations from the strict photophobic limit can arise in realistic UV completions: $ g_{a\gamma\gamma}=0 $ is a matching condition and can receive threshold corrections; therefore, a small tree-level $ g_{a\gamma\gamma} $ at collider scales is possible [16, 77]. We parametrize a controlled departure from the strict relation via

for which $ g_{a\gamma\gamma}\propto\delta $ and thus $ {\rm{Br}}(a\to\gamma\gamma)\propto\delta^2 $ (at fixed $ g_{aWW} $). For $ m_a $ above the $ WW $ threshold, $ a\to W^+W^- $ remains dominant for modest departures, while larger δ enhances $ a\to\gamma\gamma $ and dilutes $ {\rm{Br}}(a\to W^+W^-) $ (see Appendix A.1). Independently, within the strict photophobic effective field theory (EFT), the photophobic condition is preserved by renormalization group (RG) evolution up to the two-loop level, and loop/RG-induced rare channels remain at the percent level (or below) in our mass range [82]. Therefore, our results are also reported as $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ to enable straightforward recasts.

We interpret the dimension-5 operators in Eq. (1) as part of an EFT valid below a model-dependent cutoff Λ; EFT-validity, unitarity, and width considerations become particularly relevant in the high-mass/high-coupling region. A weakly-coupled UV completion suggests $ \Lambda \lesssim f_a/|c_{\widetilde{W}}| $ up to coupling/loop factors, while naive dimensional analysis (NDA) for strong dynamics suggests $ \Lambda \lesssim 4\pi f_a/|c_{\widetilde{W}}| $ [86, 87]. Using $ g_{aWW}=4c_{\widetilde{W}}/f_a $, this corresponds to

We interpret our results in the region where $ m_a $ and representative hard scales (e.g., momentum transfer) are below Λ. Partial-wave unitarity provides an independent (process-dependent) guide: in the gauge-invariant ALP EFT the strongest bound arises from $ VV\to Va $ (dominant $ J=1 $ eigenvalue), primarily constraining the SU(2)$ _L $ coupling [88]; under our convention, this leads to Eq. (A2) in Appendix A.2, which we use as a conservative validity indicator.

Although $ m_a $ cannot be fully reconstructed in the $ \ell\ell+E_T^{\rm{miss}} $ final state, our signal modelling assumes on-shell $ pp\to jja $ and we present results in the factorised form $ \sigma(pp\to jja)\times{\rm{Br}}(a\to WW) $. We therefore compute the total width in the strict photophobic benchmark with $\tt{MadGraph5\_aMC@NLO} $ and report $ \Gamma_a/m_a $ along the $ 2\sigma/5\sigma $ reach curves in Appendix A.2. The width remains small in the sensitivity region for $ m_a\lesssim 2\;\ {\rm{TeV}} $, while the extreme high-mass/high-coupling corner becomes broad and is treated as an extrapolation ¹; We therefore quote our primary results also as model-independent thresholds on $ \sigma(pp\to jja)\times{\rm{Br}}(a\to WW) $ for reinterpretation; see Appendix A.2 for the EFT-validity, unitarity, and width checks.

In the photophobic limit $ g_{a\gamma\gamma}=0 $, the ALP retains anomalous couplings to electroweak gauge bosons, notably $ aWW $, $ aZZ $, and $ a\gamma Z $. We denote $ V\in\{W,Z,\gamma\} $ and $ V'\in\{W,Z\} $. Consequently, the ALP can be produced at proton colliders in association with one vector boson ($ pp\to V' \,a $) through s-channel vector-boson exchange or with two jets ($ pp\to jj \,a $) via vector-boson-fusion (VBF)– like topologies, and it can decay as $ a\to VV' $ with branching ratios varying with the ALP mass [82, 89]. For sufficiently large $ m_a $, the $ a\to W^+W^- $ mode dominates.

For decay channel choice, the strict photophobic condition guarantees the relations in Eqs. (8) and (9). Thus, for $ m_a\gg m_W $, one has $ {\rm{Br}}(a\to W^+W^-):{\rm{Br}}(a\to ZZ): {\rm{Br}}(a\to Z\gamma)\simeq 0.65:0.16:0.19 $ [82]. Including leptonic decays, $ {\rm{Br}}(WW\to e^\pm\mu^\mp+ \not{E}_T)\simeq2.3\times10^{-2} $ while $ {\rm{Br}}(ZZ\to4\ell)\simeq 4.5\times10^{-3} $ ($ \ell=e,\mu $) [90], implying $ [{\rm{Br}}(a\to WW) {\rm{Br}}(WW\to e\mu+ \not{E}_T)]/[{\rm{Br}}(a\to ZZ){\rm{Br}}(ZZ\to4\ell)]\sim{\cal{O}}(20) $ in the heavy-mass regime. We thus focus on $ W^+W^-\to e^\pm\mu^\mp+ \not{E}_T $ as the rate-favored channel, noting that $ ZZ\to 4\ell $ is highly complementary due to its excellent mass resolution and low backgrounds.

As illustrated in Fig. 1, we consider two representative production processes: (i) an s-channel electroweak-boson exchange yielding $ pp\to jj\,a $, with the two jets predominantly arising from the decay of the associated vector boson $ V' $, and (ii) an electroweak VBF-like topology in which the ALP is radiated off the t-channel gauge bosons. In both cases, the effective a–V–$ V' $ vertices include $ aWW $, $ aZZ $, and $ a\gamma Z $. We focus on the different-flavour dilepton final state, $ W^+W^-\to e^\pm\mu^\mp\nu \bar{\nu} $, which suppresses multijet backgrounds and Drell–Yan $ Z/\gamma^\ast\to \ell^+\ell^- $ contamination while retaining high signal efficiency.

Figure 1.  Representative production mechanisms for heavy photophobic ALPs at $ pp $ colliders: s-channel exchange (left) and VBF-like topology (right), yielding $ pp\to jj\,a $ with a–V–$ V' $ vertices ($ aWW $, $ aZZ $, $ a\gamma Z $). The ALP subsequently decays to two W bosons, defining the signal $ pp\to jj\,a(\to W^\pm W^\mp) $ with $ W^\pm\to \ell^\pm\bar{\nu}/\nu $.

The VBF-like component is characterized by two forward tagging jets with large invariant mass $ m_{jj} $ and sizeable pseudo-rapidity difference $ |\Delta\eta_{jj}| $, together with reduced central hadronic activity; meanwhile, the s-channel contribution typically yields a more central dijet event with smaller rapidity separation and an invariant mass near the parent-boson mass when the jets arise from $ V' \to jj $. Accordingly, the MVA in Sec. V is trained on the inclusive (combined) signal and uses variables that are simultaneously sensitive to both topologies and kinematic observables that reconstruct the common decay $ a\to W^+W^-\to e^\pm\mu^\mp\nu \bar{\nu} $ in the $ e^\pm\mu^\mp+\not{E}_T $ final state.

Signal and background samples are generated with MadGraph5_aMC@NLO [91] using the NNPDF2.3 parton distribution functions (PDFs) [92]. The ALP interactions are implemented via a UFO model based on the linear Lagrangian [87, 93]. Parton showering and hadronization are performed with PYTHIA 8.3 [94], and the detector response is simulated with Dᴇʟᴘʜᴇs 3 [95] using HL-LHC CMS detector cards. The HL-LHC center-of-mass energy is set to $ \sqrt{s} = 14 $ TeV (7 TeV per beam). The hard-scattering process $ pp\to jj\,a $ is generated with MadGraph5_aMC@NLO, while the subsequent decays $ a\to W^+W^- $ and $ W^\pm\to \ell^\pm \nu / \bar{\nu} $ ($ \ell=e,\mu $), yielding $ a\to W^+(\to \ell^+\nu)\,W^-(\to \ell^-\bar{\nu})\,jj $, are simulated with PYTHIA 8.3.

We generate signal samples at representative ALP masses

At the HL-LHC, we generate $ 4.8 $, $ 0.8 $, and $ 4.0 $ million events for the benchmark points $ m_a=200,\,260, $ and $ 350\;\ {\rm{GeV}} $, respectively: $ 0.4 $ million events for each point with $ m_a=170,\,185,\,230\;\ {\rm{GeV}} $; $ 2.0 $ million events for each point with $ m_a=500,\,750,\,4000\;\ {\rm{GeV}} $; and $ 1.0 $ million events for each point with $ m_a=1000,\,1500,\,2000,\,2500 $, $ 3000\;\ {\rm{GeV}} $. To remain inclusive at the generator level and to match experimental conditions, we employ the following loose parton-level cuts in the MadGraph5 configuration card:

(i) Minimal transverse momentum of jets, photons, and charged leptons:

$ p_{\rm{T}}(j/\gamma/\ell^\pm) \gt 0.5\;\ {\rm{GeV}} $;

(ii) Maximal pseudorapidity:

$ |\eta(j)|<10 $ for jets and $ |\eta(\gamma/\ell^\pm)|<5 $ for photons and charged leptons;

(iii) Minimal angular separation $ \Delta R\equiv \sqrt{(\Delta\eta)^2+(\Delta\phi)^2} $ between any two objects:

$ \Delta R>0.1 $.

These generator-level requirements are deliberately much looser than the requirements for the analysis selection and are never applied as analysis cuts. They avoid generator-induced inefficiencies (i.e., an artificial suppression of soft/forward configurations) and only regulate infrared-sensitive regions for technical stability [91]. Event realism is achieved downstream: parton showering and hadronization model QCD radiation. The detector simulation emulates experimental acceptance and resolution [94, 95]. Accordingly, the fiducial phase space and physics acceptance are defined by the detector simulation together with the selection in Sec. V, in line with LHC fiducial-measurement practice.

Figure 2 shows the production cross section for $ pp\to jj\,a(\to W^\pm W^\mp) $ at $ \sqrt{s}=14\;\ {\rm{TeV}} $ under the benchmark couplings $ g_{a\gamma\gamma}=0 $ and $ g_{aWW}=1\;\ {\rm{TeV}}^{-1} $. The two-body $ WW $ decay channel opens rapidly once $ m_a\gtrsim 2m_W\simeq 160 \ {\rm{GeV}} $, causing a steep rise in $ {\rm{Br}}(a\to W^+W^-) $; for example, at $ m_a=185\;\ {\rm{GeV}} $, the branching fraction reaches ~49% [82, 89]. Meanwhile, the inclusive production cross section $ \sigma(pp\to jj\,a) $ decreases with increasing $ m_a $ due to diminishing PDFs, such that $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $ peaks around $ m_a\sim200\;\ {\rm{GeV}} $.

Figure 2.  Production cross section of $ pp\to jj\,a(\to W^\pm W^\mp) $ at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $) for $ g_{a\gamma\gamma}=0 $ and $ g_{aWW}=1\;\ {\rm{TeV}}^{-1} $ as a function of the ALP mass $ m_a $.

In the photophobic limit $ g_{a\gamma\gamma}=0 $, the ALP retains anomalous couplings to electroweak gauge bosons, notably $ aWW $, $ aZZ $, and $ a\gamma Z $. We denote $ V\in\{W,Z,\gamma\} $ and $ V'\in\{W,Z\} $. Consequently, the ALP can be produced at proton colliders in association with one vector boson ($ pp\to V' \,a $) through s-channel vector-boson exchange or with two jets ($ pp\to jj \,a $) via vector-boson-fusion (VBF)– like topologies, and it can decay as $ a\to VV' $ with branching ratios varying with the ALP mass [82, 89]. For sufficiently large $ m_a $, the $ a\to W^+W^- $ mode dominates.

For decay channel choice, the strict photophobic condition guarantees the relations in Eqs. (8) and (9). Thus, for $ m_a\gg m_W $, one has $ {\rm{Br}}(a\to W^+W^-):{\rm{Br}}(a\to ZZ): {\rm{Br}}(a\to Z\gamma)\simeq 0.65:0.16:0.19 $ [82]. Including leptonic decays, $ {\rm{Br}}(WW\to e^\pm\mu^\mp+ \not{E}_T)\simeq2.3\times10^{-2} $ while $ {\rm{Br}}(ZZ\to4\ell)\simeq 4.5\times10^{-3} $ ($ \ell=e,\mu $) [90], implying $ [{\rm{Br}}(a\to WW) {\rm{Br}}(WW\to e\mu+ \not{E}_T)]/[{\rm{Br}}(a\to ZZ){\rm{Br}}(ZZ\to4\ell)]\sim{\cal{O}}(20) $ in the heavy-mass regime. We thus focus on $ W^+W^-\to e^\pm\mu^\mp+ \not{E}_T $ as the rate-favored channel, noting that $ ZZ\to 4\ell $ is highly complementary due to its excellent mass resolution and low backgrounds.

As illustrated in Fig. 1, we consider two representative production processes: (i) an s-channel electroweak-boson exchange yielding $ pp\to jj\,a $, with the two jets predominantly arising from the decay of the associated vector boson $ V' $, and (ii) an electroweak VBF-like topology in which the ALP is radiated off the t-channel gauge bosons. In both cases, the effective a–V–$ V' $ vertices include $ aWW $, $ aZZ $, and $ a\gamma Z $. We focus on the different-flavour dilepton final state, $ W^+W^-\to e^\pm\mu^\mp\nu \bar{\nu} $, which suppresses multijet backgrounds and Drell–Yan $ Z/\gamma^\ast\to \ell^+\ell^- $ contamination while retaining high signal efficiency.

Figure 1.  Representative production mechanisms for heavy photophobic ALPs at $ pp $ colliders: s-channel exchange (left) and VBF-like topology (right), yielding $ pp\to jj\,a $ with a–V–$ V' $ vertices ($ aWW $, $ aZZ $, $ a\gamma Z $). The ALP subsequently decays to two W bosons, defining the signal $ pp\to jj\,a(\to W^\pm W^\mp) $ with $ W^\pm\to \ell^\pm\bar{\nu}/\nu $.

The VBF-like component is characterized by two forward tagging jets with large invariant mass $ m_{jj} $ and sizeable pseudo-rapidity difference $ |\Delta\eta_{jj}| $, together with reduced central hadronic activity; meanwhile, the s-channel contribution typically yields a more central dijet event with smaller rapidity separation and an invariant mass near the parent-boson mass when the jets arise from $ V' \to jj $. Accordingly, the MVA in Sec. V is trained on the inclusive (combined) signal and uses variables that are simultaneously sensitive to both topologies and kinematic observables that reconstruct the common decay $ a\to W^+W^-\to e^\pm\mu^\mp\nu \bar{\nu} $ in the $ e^\pm\mu^\mp+\not{E}_T $ final state.

Signal and background samples are generated with MadGraph5_aMC@NLO [91] using the NNPDF2.3 parton distribution functions (PDFs) [92]. The ALP interactions are implemented via a UFO model based on the linear Lagrangian [87, 93]. Parton showering and hadronization are performed with PYTHIA 8.3 [94], and the detector response is simulated with Dᴇʟᴘʜᴇs 3 [95] using HL-LHC CMS detector cards. The HL-LHC center-of-mass energy is set to $ \sqrt{s} = 14 $ TeV (7 TeV per beam). The hard-scattering process $ pp\to jj\,a $ is generated with MadGraph5_aMC@NLO, while the subsequent decays $ a\to W^+W^- $ and $ W^\pm\to \ell^\pm \nu / \bar{\nu} $ ($ \ell=e,\mu $), yielding $ a\to W^+(\to \ell^+\nu)\,W^-(\to \ell^-\bar{\nu})\,jj $, are simulated with PYTHIA 8.3.

We generate signal samples at representative ALP masses

At the HL-LHC, we generate $ 4.8 $, $ 0.8 $, and $ 4.0 $ million events for the benchmark points $ m_a=200,\,260, $ and $ 350\;\ {\rm{GeV}} $, respectively: $ 0.4 $ million events for each point with $ m_a=170,\,185,\,230\;\ {\rm{GeV}} $; $ 2.0 $ million events for each point with $ m_a=500,\,750,\,4000\;\ {\rm{GeV}} $; and $ 1.0 $ million events for each point with $ m_a=1000,\,1500,\,2000,\,2500 $, $ 3000\;\ {\rm{GeV}} $. To remain inclusive at the generator level and to match experimental conditions, we employ the following loose parton-level cuts in the MadGraph5 configuration card:

(i) Minimal transverse momentum of jets, photons, and charged leptons:

$ p_{\rm{T}}(j/\gamma/\ell^\pm) \gt 0.5\;\ {\rm{GeV}} $;

(ii) Maximal pseudorapidity:

$ |\eta(j)|<10 $ for jets and $ |\eta(\gamma/\ell^\pm)|<5 $ for photons and charged leptons;

(iii) Minimal angular separation $ \Delta R\equiv \sqrt{(\Delta\eta)^2+(\Delta\phi)^2} $ between any two objects:

$ \Delta R>0.1 $.

These generator-level requirements are deliberately much looser than the requirements for the analysis selection and are never applied as analysis cuts. They avoid generator-induced inefficiencies (i.e., an artificial suppression of soft/forward configurations) and only regulate infrared-sensitive regions for technical stability [91]. Event realism is achieved downstream: parton showering and hadronization model QCD radiation. The detector simulation emulates experimental acceptance and resolution [94, 95]. Accordingly, the fiducial phase space and physics acceptance are defined by the detector simulation together with the selection in Sec. V, in line with LHC fiducial-measurement practice.

Figure 2 shows the production cross section for $ pp\to jj\,a(\to W^\pm W^\mp) $ at $ \sqrt{s}=14\;\ {\rm{TeV}} $ under the benchmark couplings $ g_{a\gamma\gamma}=0 $ and $ g_{aWW}=1\;\ {\rm{TeV}}^{-1} $. The two-body $ WW $ decay channel opens rapidly once $ m_a\gtrsim 2m_W\simeq 160 \ {\rm{GeV}} $, causing a steep rise in $ {\rm{Br}}(a\to W^+W^-) $; for example, at $ m_a=185\;\ {\rm{GeV}} $, the branching fraction reaches ~49% [82, 89]. Meanwhile, the inclusive production cross section $ \sigma(pp\to jj\,a) $ decreases with increasing $ m_a $ due to diminishing PDFs, such that $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $ peaks around $ m_a\sim200\;\ {\rm{GeV}} $.

Figure 2.  Production cross section of $ pp\to jj\,a(\to W^\pm W^\mp) $ at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $) for $ g_{a\gamma\gamma}=0 $ and $ g_{aWW}=1\;\ {\rm{TeV}}^{-1} $ as a function of the ALP mass $ m_a $.

In the photophobic limit $ g_{a\gamma\gamma}=0 $, the ALP retains anomalous couplings to electroweak gauge bosons, notably $ aWW $, $ aZZ $, and $ a\gamma Z $. We denote $ V\in\{W,Z,\gamma\} $ and $ V'\in\{W,Z\} $. Consequently, the ALP can be produced at proton colliders in association with one vector boson ($ pp\to V' \,a $) through s-channel vector-boson exchange or with two jets ($ pp\to jj \,a $) via vector-boson-fusion (VBF)– like topologies, and it can decay as $ a\to VV' $ with branching ratios varying with the ALP mass [82, 89]. For sufficiently large $ m_a $, the $ a\to W^+W^- $ mode dominates.

For decay channel choice, the strict photophobic condition guarantees the relations in Eqs. (8) and (9). Thus, for $ m_a\gg m_W $, one has $ {\rm{Br}}(a\to W^+W^-):{\rm{Br}}(a\to ZZ): {\rm{Br}}(a\to Z\gamma)\simeq 0.65:0.16:0.19 $ [82]. Including leptonic decays, $ {\rm{Br}}(WW\to e^\pm\mu^\mp+ \not{E}_T)\simeq2.3\times10^{-2} $ while $ {\rm{Br}}(ZZ\to4\ell)\simeq 4.5\times10^{-3} $ ($ \ell=e,\mu $) [90], implying $ [{\rm{Br}}(a\to WW) {\rm{Br}}(WW\to e\mu+ \not{E}_T)]/[{\rm{Br}}(a\to ZZ){\rm{Br}}(ZZ\to4\ell)]\sim{\cal{O}}(20) $ in the heavy-mass regime. We thus focus on $ W^+W^-\to e^\pm\mu^\mp+ \not{E}_T $ as the rate-favored channel, noting that $ ZZ\to 4\ell $ is highly complementary due to its excellent mass resolution and low backgrounds.

As illustrated in Fig. 1, we consider two representative production processes: (i) an s-channel electroweak-boson exchange yielding $ pp\to jj\,a $, with the two jets predominantly arising from the decay of the associated vector boson $ V' $, and (ii) an electroweak VBF-like topology in which the ALP is radiated off the t-channel gauge bosons. In both cases, the effective a–V–$ V' $ vertices include $ aWW $, $ aZZ $, and $ a\gamma Z $. We focus on the different-flavour dilepton final state, $ W^+W^-\to e^\pm\mu^\mp\nu \bar{\nu} $, which suppresses multijet backgrounds and Drell–Yan $ Z/\gamma^\ast\to \ell^+\ell^- $ contamination while retaining high signal efficiency.

Figure 1.  Representative production mechanisms for heavy photophobic ALPs at $ pp $ colliders: s-channel exchange (left) and VBF-like topology (right), yielding $ pp\to jj\,a $ with a–V–$ V' $ vertices ($ aWW $, $ aZZ $, $ a\gamma Z $). The ALP subsequently decays to two W bosons, defining the signal $ pp\to jj\,a(\to W^\pm W^\mp) $ with $ W^\pm\to \ell^\pm\bar{\nu}/\nu $.

The VBF-like component is characterized by two forward tagging jets with large invariant mass $ m_{jj} $ and sizeable pseudo-rapidity difference $ |\Delta\eta_{jj}| $, together with reduced central hadronic activity; meanwhile, the s-channel contribution typically yields a more central dijet event with smaller rapidity separation and an invariant mass near the parent-boson mass when the jets arise from $ V' \to jj $. Accordingly, the MVA in Sec. V is trained on the inclusive (combined) signal and uses variables that are simultaneously sensitive to both topologies and kinematic observables that reconstruct the common decay $ a\to W^+W^-\to e^\pm\mu^\mp\nu \bar{\nu} $ in the $ e^\pm\mu^\mp+\not{E}_T $ final state.

Signal and background samples are generated with MadGraph5_aMC@NLO [91] using the NNPDF2.3 parton distribution functions (PDFs) [92]. The ALP interactions are implemented via a UFO model based on the linear Lagrangian [87, 93]. Parton showering and hadronization are performed with PYTHIA 8.3 [94], and the detector response is simulated with Dᴇʟᴘʜᴇs 3 [95] using HL-LHC CMS detector cards. The HL-LHC center-of-mass energy is set to $ \sqrt{s} = 14 $ TeV (7 TeV per beam). The hard-scattering process $ pp\to jj\,a $ is generated with MadGraph5_aMC@NLO, while the subsequent decays $ a\to W^+W^- $ and $ W^\pm\to \ell^\pm \nu / \bar{\nu} $ ($ \ell=e,\mu $), yielding $ a\to W^+(\to \ell^+\nu)\,W^-(\to \ell^-\bar{\nu})\,jj $, are simulated with PYTHIA 8.3.

We generate signal samples at representative ALP masses

At the HL-LHC, we generate $ 4.8 $, $ 0.8 $, and $ 4.0 $ million events for the benchmark points $ m_a=200,\,260, $ and $ 350\;\ {\rm{GeV}} $, respectively: $ 0.4 $ million events for each point with $ m_a=170,\,185,\,230\;\ {\rm{GeV}} $; $ 2.0 $ million events for each point with $ m_a=500,\,750,\,4000\;\ {\rm{GeV}} $; and $ 1.0 $ million events for each point with $ m_a=1000,\,1500,\,2000,\,2500 $, $ 3000\;\ {\rm{GeV}} $. To remain inclusive at the generator level and to match experimental conditions, we employ the following loose parton-level cuts in the MadGraph5 configuration card:

(i) Minimal transverse momentum of jets, photons, and charged leptons:

$ p_{\rm{T}}(j/\gamma/\ell^\pm) \gt 0.5\;\ {\rm{GeV}} $;

(ii) Maximal pseudorapidity:

$ |\eta(j)|<10 $ for jets and $ |\eta(\gamma/\ell^\pm)|<5 $ for photons and charged leptons;

(iii) Minimal angular separation $ \Delta R\equiv \sqrt{(\Delta\eta)^2+(\Delta\phi)^2} $ between any two objects:

$ \Delta R>0.1 $.

These generator-level requirements are deliberately much looser than the requirements for the analysis selection and are never applied as analysis cuts. They avoid generator-induced inefficiencies (i.e., an artificial suppression of soft/forward configurations) and only regulate infrared-sensitive regions for technical stability [91]. Event realism is achieved downstream: parton showering and hadronization model QCD radiation. The detector simulation emulates experimental acceptance and resolution [94, 95]. Accordingly, the fiducial phase space and physics acceptance are defined by the detector simulation together with the selection in Sec. V, in line with LHC fiducial-measurement practice.

Figure 2 shows the production cross section for $ pp\to jj\,a(\to W^\pm W^\mp) $ at $ \sqrt{s}=14\;\ {\rm{TeV}} $ under the benchmark couplings $ g_{a\gamma\gamma}=0 $ and $ g_{aWW}=1\;\ {\rm{TeV}}^{-1} $. The two-body $ WW $ decay channel opens rapidly once $ m_a\gtrsim 2m_W\simeq 160 \ {\rm{GeV}} $, causing a steep rise in $ {\rm{Br}}(a\to W^+W^-) $; for example, at $ m_a=185\;\ {\rm{GeV}} $, the branching fraction reaches ~49% [82, 89]. Meanwhile, the inclusive production cross section $ \sigma(pp\to jj\,a) $ decreases with increasing $ m_a $ due to diminishing PDFs, such that $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $ peaks around $ m_a\sim200\;\ {\rm{GeV}} $.

Figure 2.  Production cross section of $ pp\to jj\,a(\to W^\pm W^\mp) $ at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $) for $ g_{a\gamma\gamma}=0 $ and $ g_{aWW}=1\;\ {\rm{TeV}}^{-1} $ as a function of the ALP mass $ m_a $.

The signal targets two oppositely charged, different-flavour leptons, with at least two jets and moderate missing transverse momentum, i.e., $ e^\pm\mu^\mp+jj+\not{E}_T $. In the $ e\mu $ channel, Drell–Yan contamination is strongly suppressed; the dominant backgrounds are well known to include diboson and top-quark production, as well as processes with nonprompt or misidentified leptons. We group them by physics origin and reducibility at the detector level.

Irreducible (same visible final state)

● $ W^+W^-jj $ (QCD+EW). Inclusive $ WW $ production accompanied by jets, dominated by $ q\bar q \to WW $ with additional QCD radiation [96]. After leptonic W decays, the visible final state coincides with the signal ($ e^\pm\mu^\mp+ jj+\not{E}_T $).

● $ t\bar t $ (dileptonic). When both W bosons decay leptonically, $ t\bar t \to \ell^+\nu\,\ell^-\bar{\nu}\,b\bar b $ can be indistinguishable from the signal up to global kinematics if b jets are not correctly tagged; thus, it is effectively irreducible in this channel.

Reducible (lost/misidentified objects)

● $ WZjj $, $ ZZjj $. Either when one prompt lepton is not reconstructed or fails identification/isolation, or when $ Z \to \tau\tau $ produces a pair $ e\mu $.

● $ W^\pm jj $. Contributes via one prompt lepton from the W and one nonprompt or misidentified lepton (heavy-flavour decays, hadrons mis-identification, photon conversions).

● $ Zjj $. Suppressed in $ e\mu $ but nonzero through a misidentified lepton, a lost lepton plus a nonprompt lepton, or $ Z \to \tau\tau \to e\mu+X $.

The $ pp \to W^+ W^- jj $ background is generated at matrix-element level with MadGraph5_aMC@NLO, while the subsequent decays $ W \to \ell^\pm \nu_\ell $ ($ \ell = e,\mu $) and parton showering are modeled with PYTHIA 8.3. The $ t\bar{t} $, $ WZjj $, and $ ZZjj $ samples are likewise produced with MadGraph5_aMC@NLO and interfaced to PYTHIA 8.3, which handles all kinematically allowed W- and Z-boson decays. Single-boson processes $ W^\pm(\to \ell^\pm \nu_\ell)jj $ and $ Z(\to \ell^+\ell^-)jj $ with $ \ell = e,\mu $ are generated fully differentially in MadGraph5_aMC@NLO, keeping the leptonic W and Z decays explicit in the matrix element.

The expected number of events scales with luminosity as

where $ \sigma_{\rm{pro}} $ is the production cross section, $ {\cal{L}} $ is the integrated luminosity, and $ \varepsilon_{\rm{sel}} $ is the net acceptance efficiency multiplied by the efficiencies from object reconstruction and analysis selections. Table 1 lists the generator-level inclusive production cross sections of background processes at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $). In practice, although the $ Wjj $ and $ Zjj $ processes dominate at the generator level (see Table 1), their contribution after the dilepton+ $ \not{E}_T $ selections is strongly reduced, while $ WWjj $ and $ t\bar t $ processes remain the leading backgrounds in the $ e\mu $ channel [96].

The signal targets two oppositely charged, different-flavour leptons, with at least two jets and moderate missing transverse momentum, i.e., $ e^\pm\mu^\mp+jj+\not{E}_T $. In the $ e\mu $ channel, Drell–Yan contamination is strongly suppressed; the dominant backgrounds are well known to include diboson and top-quark production, as well as processes with nonprompt or misidentified leptons. We group them by physics origin and reducibility at the detector level.

Irreducible (same visible final state)

● $ W^+W^-jj $ (QCD+EW). Inclusive $ WW $ production accompanied by jets, dominated by $ q\bar q \to WW $ with additional QCD radiation [96]. After leptonic W decays, the visible final state coincides with the signal ($ e^\pm\mu^\mp+ jj+\not{E}_T $).

● $ t\bar t $ (dileptonic). When both W bosons decay leptonically, $ t\bar t \to \ell^+\nu\,\ell^-\bar{\nu}\,b\bar b $ can be indistinguishable from the signal up to global kinematics if b jets are not correctly tagged; thus, it is effectively irreducible in this channel.

Reducible (lost/misidentified objects)

● $ WZjj $, $ ZZjj $. Either when one prompt lepton is not reconstructed or fails identification/isolation, or when $ Z \to \tau\tau $ produces a pair $ e\mu $.

● $ W^\pm jj $. Contributes via one prompt lepton from the W and one nonprompt or misidentified lepton (heavy-flavour decays, hadrons mis-identification, photon conversions).

● $ Zjj $. Suppressed in $ e\mu $ but nonzero through a misidentified lepton, a lost lepton plus a nonprompt lepton, or $ Z \to \tau\tau \to e\mu+X $.

The $ pp \to W^+ W^- jj $ background is generated at matrix-element level with MadGraph5_aMC@NLO, while the subsequent decays $ W \to \ell^\pm \nu_\ell $ ($ \ell = e,\mu $) and parton showering are modeled with PYTHIA 8.3. The $ t\bar{t} $, $ WZjj $, and $ ZZjj $ samples are likewise produced with MadGraph5_aMC@NLO and interfaced to PYTHIA 8.3, which handles all kinematically allowed W- and Z-boson decays. Single-boson processes $ W^\pm(\to \ell^\pm \nu_\ell)jj $ and $ Z(\to \ell^+\ell^-)jj $ with $ \ell = e,\mu $ are generated fully differentially in MadGraph5_aMC@NLO, keeping the leptonic W and Z decays explicit in the matrix element.

The expected number of events scales with luminosity as

where $ \sigma_{\rm{pro}} $ is the production cross section, $ {\cal{L}} $ is the integrated luminosity, and $ \varepsilon_{\rm{sel}} $ is the net acceptance efficiency multiplied by the efficiencies from object reconstruction and analysis selections. Table 1 lists the generator-level inclusive production cross sections of background processes at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $). In practice, although the $ Wjj $ and $ Zjj $ processes dominate at the generator level (see Table 1), their contribution after the dilepton+ $ \not{E}_T $ selections is strongly reduced, while $ WWjj $ and $ t\bar t $ processes remain the leading backgrounds in the $ e\mu $ channel [96].

The signal targets two oppositely charged, different-flavour leptons, with at least two jets and moderate missing transverse momentum, i.e., $ e^\pm\mu^\mp+jj+\not{E}_T $. In the $ e\mu $ channel, Drell–Yan contamination is strongly suppressed; the dominant backgrounds are well known to include diboson and top-quark production, as well as processes with nonprompt or misidentified leptons. We group them by physics origin and reducibility at the detector level.

Irreducible (same visible final state)

● $ W^+W^-jj $ (QCD+EW). Inclusive $ WW $ production accompanied by jets, dominated by $ q\bar q \to WW $ with additional QCD radiation [96]. After leptonic W decays, the visible final state coincides with the signal ($ e^\pm\mu^\mp+ jj+\not{E}_T $).

● $ t\bar t $ (dileptonic). When both W bosons decay leptonically, $ t\bar t \to \ell^+\nu\,\ell^-\bar{\nu}\,b\bar b $ can be indistinguishable from the signal up to global kinematics if b jets are not correctly tagged; thus, it is effectively irreducible in this channel.

Reducible (lost/misidentified objects)

● $ WZjj $, $ ZZjj $. Either when one prompt lepton is not reconstructed or fails identification/isolation, or when $ Z \to \tau\tau $ produces a pair $ e\mu $.

● $ W^\pm jj $. Contributes via one prompt lepton from the W and one nonprompt or misidentified lepton (heavy-flavour decays, hadrons mis-identification, photon conversions).

● $ Zjj $. Suppressed in $ e\mu $ but nonzero through a misidentified lepton, a lost lepton plus a nonprompt lepton, or $ Z \to \tau\tau \to e\mu+X $.

The $ pp \to W^+ W^- jj $ background is generated at matrix-element level with MadGraph5_aMC@NLO, while the subsequent decays $ W \to \ell^\pm \nu_\ell $ ($ \ell = e,\mu $) and parton showering are modeled with PYTHIA 8.3. The $ t\bar{t} $, $ WZjj $, and $ ZZjj $ samples are likewise produced with MadGraph5_aMC@NLO and interfaced to PYTHIA 8.3, which handles all kinematically allowed W- and Z-boson decays. Single-boson processes $ W^\pm(\to \ell^\pm \nu_\ell)jj $ and $ Z(\to \ell^+\ell^-)jj $ with $ \ell = e,\mu $ are generated fully differentially in MadGraph5_aMC@NLO, keeping the leptonic W and Z decays explicit in the matrix element.

The expected number of events scales with luminosity as

where $ \sigma_{\rm{pro}} $ is the production cross section, $ {\cal{L}} $ is the integrated luminosity, and $ \varepsilon_{\rm{sel}} $ is the net acceptance efficiency multiplied by the efficiencies from object reconstruction and analysis selections. Table 1 lists the generator-level inclusive production cross sections of background processes at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $). In practice, although the $ Wjj $ and $ Zjj $ processes dominate at the generator level (see Table 1), their contribution after the dilepton+ $ \not{E}_T $ selections is strongly reduced, while $ WWjj $ and $ t\bar t $ processes remain the leading backgrounds in the $ e\mu $ channel [96].

We target final states with two oppositely charged, different-flavour leptons, with at least two jets and moderate missing transverse momentum. Throughout the analysis, objects are ordered by transverse momentum, and the leading (subleading) jet is denoted by $ j_1 $ ($ j_2 $). The following preselection is applied to select the $ e^\pm\mu^\mp+jj+ \not{E}_T $ final state and suppress backgrounds at the first stage:

(i) Exactly two leptons: $ N(\ell)=2 $;

(ii) Lepton transverse momenta: $ p_T(\ell)\ge 10\;\ {\rm{GeV}} $ for both leptons;

(iii) Different flavour: exactly one electron and one muon, $ N(e)=1 $ and $ N(\mu)=1 $;

(iv) Opposite charge: $ Q(e)+Q(\mu)=0 $;

(v) Jet multiplicity: at least two jets, $ N(j)\ge 2 $;

(vi) Leading-jet thresholds: $ p_T(j_1)\ge 30\;\ {\rm{GeV}} $ and $ p_T(j_2)\ge 30\;\ {\rm{GeV}} $;

(vii) Top-quark suppression: zero b-tagged jets, $ N(j_b)=0 $.

This opposite-sign $ e\mu $ requirement strongly reduces Drell–Yan backgrounds, while the b-jet veto suppresses top-quark processes, in line with LHC $ WW\to e^\pm\nu\mu^\mp \bar{\nu} $ measurements [96].

As summarized in Table 2, the $ Wjj $ and $ Zjj $ processes have the lowest preselection efficiencies because of the tight dilepton requirement. Diboson $ WZjj $ and $ ZZjj $ remain as subleading reducible backgrounds primarily through lost-lepton or τ-lepton decays, while $ WWjj $ and $ t\bar t $ are effectively irreducible up to global kinematics; the latter is also strongly reduced by the b-tag veto.

Table 3 lists the event yields after sequentially applying preselection criteria (i)–(vii) for the signal with benchmark $ m_a=750\;\ {\rm{GeV}} $ and for background processes at the HL-LHC with $ \sqrt{s}=14\;\ {\rm{TeV}} $ and $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $. With the opposite-sign $ e\mu $ requirement, the large inclusive $ Wjj $ and $ Zjj $ rates are greatly reduced, while $ WWjj $ and $ t\bar t $ remain the leading backgrounds in the $ e\mu $ channel [96].

We target final states with two oppositely charged, different-flavour leptons, with at least two jets and moderate missing transverse momentum. Throughout the analysis, objects are ordered by transverse momentum, and the leading (subleading) jet is denoted by $ j_1 $ ($ j_2 $). The following preselection is applied to select the $ e^\pm\mu^\mp+jj+ \not{E}_T $ final state and suppress backgrounds at the first stage:

(i) Exactly two leptons: $ N(\ell)=2 $;

(ii) Lepton transverse momenta: $ p_T(\ell)\ge 10\;\ {\rm{GeV}} $ for both leptons;

(iii) Different flavour: exactly one electron and one muon, $ N(e)=1 $ and $ N(\mu)=1 $;

(iv) Opposite charge: $ Q(e)+Q(\mu)=0 $;

(v) Jet multiplicity: at least two jets, $ N(j)\ge 2 $;

(vi) Leading-jet thresholds: $ p_T(j_1)\ge 30\;\ {\rm{GeV}} $ and $ p_T(j_2)\ge 30\;\ {\rm{GeV}} $;

(vii) Top-quark suppression: zero b-tagged jets, $ N(j_b)=0 $.

This opposite-sign $ e\mu $ requirement strongly reduces Drell–Yan backgrounds, while the b-jet veto suppresses top-quark processes, in line with LHC $ WW\to e^\pm\nu\mu^\mp \bar{\nu} $ measurements [96].

As summarized in Table 2, the $ Wjj $ and $ Zjj $ processes have the lowest preselection efficiencies because of the tight dilepton requirement. Diboson $ WZjj $ and $ ZZjj $ remain as subleading reducible backgrounds primarily through lost-lepton or τ-lepton decays, while $ WWjj $ and $ t\bar t $ are effectively irreducible up to global kinematics; the latter is also strongly reduced by the b-tag veto.

Table 3 lists the event yields after sequentially applying preselection criteria (i)–(vii) for the signal with benchmark $ m_a=750\;\ {\rm{GeV}} $ and for background processes at the HL-LHC with $ \sqrt{s}=14\;\ {\rm{TeV}} $ and $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $. With the opposite-sign $ e\mu $ requirement, the large inclusive $ Wjj $ and $ Zjj $ rates are greatly reduced, while $ WWjj $ and $ t\bar t $ remain the leading backgrounds in the $ e\mu $ channel [96].

We target final states with two oppositely charged, different-flavour leptons, with at least two jets and moderate missing transverse momentum. Throughout the analysis, objects are ordered by transverse momentum, and the leading (subleading) jet is denoted by $ j_1 $ ($ j_2 $). The following preselection is applied to select the $ e^\pm\mu^\mp+jj+ \not{E}_T $ final state and suppress backgrounds at the first stage:

(i) Exactly two leptons: $ N(\ell)=2 $;

(ii) Lepton transverse momenta: $ p_T(\ell)\ge 10\;\ {\rm{GeV}} $ for both leptons;

(iii) Different flavour: exactly one electron and one muon, $ N(e)=1 $ and $ N(\mu)=1 $;

(iv) Opposite charge: $ Q(e)+Q(\mu)=0 $;

(v) Jet multiplicity: at least two jets, $ N(j)\ge 2 $;

(vi) Leading-jet thresholds: $ p_T(j_1)\ge 30\;\ {\rm{GeV}} $ and $ p_T(j_2)\ge 30\;\ {\rm{GeV}} $;

(vii) Top-quark suppression: zero b-tagged jets, $ N(j_b)=0 $.

This opposite-sign $ e\mu $ requirement strongly reduces Drell–Yan backgrounds, while the b-jet veto suppresses top-quark processes, in line with LHC $ WW\to e^\pm\nu\mu^\mp \bar{\nu} $ measurements [96].

As summarized in Table 2, the $ Wjj $ and $ Zjj $ processes have the lowest preselection efficiencies because of the tight dilepton requirement. Diboson $ WZjj $ and $ ZZjj $ remain as subleading reducible backgrounds primarily through lost-lepton or τ-lepton decays, while $ WWjj $ and $ t\bar t $ are effectively irreducible up to global kinematics; the latter is also strongly reduced by the b-tag veto.

Table 3 lists the event yields after sequentially applying preselection criteria (i)–(vii) for the signal with benchmark $ m_a=750\;\ {\rm{GeV}} $ and for background processes at the HL-LHC with $ \sqrt{s}=14\;\ {\rm{TeV}} $ and $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $. With the opposite-sign $ e\mu $ requirement, the large inclusive $ Wjj $ and $ Zjj $ rates are greatly reduced, while $ WWjj $ and $ t\bar t $ remain the leading backgrounds in the $ e\mu $ channel [96].

After the preselection, we employ a multivariate analysis (MVA) to further separate signal from background. The input observables are as follows:

(i) Object kinematics: energy, transverse momentum, pseudorapidity, and azimuth of the final-state objects, $ E(O) $, $ p_T(O) $, $ \eta(O) $, $ \phi(O) $ with $ O=\mu,\,e,\,j_1,\,j_2 $.

(ii) Missing transverse momentum (magnitude and azimuth): $ \not{E}_T $ and $ \phi( \not{E}_T) $.

(iii) Among all jet pairs, the pair with the minimal angular separation: $ \Delta R (j, j')_{\rm{min}} $ and the corresponding dijet mass $ m(j+j')_{ {\rm{min}} \Delta R} $, where $ \Delta R = \sqrt{(\Delta \eta)^2 + (\Delta \varphi)^2} $.

(iv) Among all jet pairs, the pair with the maximal pseudorapidity separation: $ \Delta \eta(j, j')_{\rm{max}} $ and the corresponding dijet mass $ m(j + j')_{ {\rm{max}} \Delta \eta} $.

(v) The invariant mass of the first two leading jets and the jet multiplicity: $ m(j_{1}+j_{2}) $ and $ N(j) $.

(vi) Angular separations between key objects: $ \Delta R(j_1,j_2) $, $ \Delta\eta(j_1,j_2) $, $ \Delta R(e,\mu) $, $ \Delta R(\mu,j_1) $, $ \Delta R(\mu,j_2) $, $ \Delta R(e,j_1) $, and $ \Delta R(e,j_2) $.

(vii) Transverse masses of leptons and $ \not{E}_T $: $ m_{\rm{T}}(\mu + \not{E}_T) $, $ m_{\rm{T}}(e + \not{E}_T) $, and $ m_{\rm{T}}(\mu + e + \not{E}_T) $. Here, the transverse mass $ m_{\mathrm T} $ of a visible system and the missing transverse momentum is defined as $ m_{\mathrm T}^2 = \big(E_{\mathrm T}^{\rm{vis}} + \not{E}_T \big)^2 - \bigl|\vec{p}_{\mathrm T}^{\,\rm vis} + \vec{p}_{\mathrm T}^{\,\rm miss}\bigr|^2 , $ where $ E_{\mathrm T}^{\rm{vis}} = \sqrt{|\vec{p}_{\mathrm T}^{\,\rm vis}|^2 + (m^{\rm{vis}})^2} $ is the transverse energy of the visible system and $ \not{E}_T = |\vec{p}_{\mathrm T}^{\,\rm miss}| $ assumes a massless invisible system.

Observables (iv)–(vi) exploit the VBF-like topology: two forward tagging jets, large $ \Delta\eta_{jj} $, and high $ m_{jj} $ with reduced central hadronic activity, while (iii) captures s-channel topologies where jets tend to be closer in phase space. This division follows standard LHC practices for VBF/di-jet tagging. Observables (vii) are related to the kinematics of the decay $ a\to W^+W^-\to e^\pm\mu^\mp\nu \bar{\nu} $ in the $ e^\pm\mu^\mp+ \not{E}_T $ final state.

We implement the MVA using a BDT via the TMVA package integrated in ROOT, adopting the default BDT configuration for classification [97]. This setup is used to maximize background rejection while retaining high signal efficiency. Figure 3 shows the BDT responses for two benchmark masses, $ m_a=750\;\ {\rm{GeV}} $ and $ m_a=2000\;\ {\rm{GeV}} $. As $ m_a $ increases, the underlying kinematics becomes more boosted and the separation between signal and background events improves, as shown by the shift of the signal (background) distribution towards higher (lower) BDT scores.

Figure 3.  (color online) BDT response distributions for signal (shaded) and total background (outline) at the HL-LHC with $ \sqrt{s}= $$ 14\;\ {\rm{TeV}} $, for $ m_a=750\;\ {\rm{GeV}} $ (left) and $ m_a=2000\;\ {\rm{GeV}} $ (right).

The BDT score threshold is optimized independently at each $ m_a $ by maximizing the statistical significance

where $ N_s $ and $ N_b $ denote the expected signal and total background yields after preselection and the BDT cut, respectively [98].

For completeness, Appendix B presents representative input-observable distributions after preselection for benchmark $ m_a=750\;\ {\rm{GeV}} $, while Appendix C shows BDT responses across the full $ m_a $ range for various $ m_a $ cases. Appendix D reports the preselection and BDT efficiencies for signals and backgrounds at $ \sqrt{s}=14\;\ {\rm{TeV}} $. In that table, a dash ("−") indicates that the surviving yield of a given background is negligible at $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $. Despite their large inclusive cross sections, the $ W^\pm(\to\ell^\pm\nu)jj $ and $ Z(\to\ell^+\ell^-)jj $ background events are almost completely removed by the opposite-sign $ e\mu $ requirement and multivariate selection, so they are omitted from the appendix table for brevity.

After the preselection, we employ a multivariate analysis (MVA) to further separate signal from background. The input observables are as follows:

(i) Object kinematics: energy, transverse momentum, pseudorapidity, and azimuth of the final-state objects, $ E(O) $, $ p_T(O) $, $ \eta(O) $, $ \phi(O) $ with $ O=\mu,\,e,\,j_1,\,j_2 $.

(ii) Missing transverse momentum (magnitude and azimuth): $ \not{E}_T $ and $ \phi( \not{E}_T) $.

(iii) Among all jet pairs, the pair with the minimal angular separation: $ \Delta R (j, j')_{\rm{min}} $ and the corresponding dijet mass $ m(j+j')_{ {\rm{min}} \Delta R} $, where $ \Delta R = \sqrt{(\Delta \eta)^2 + (\Delta \varphi)^2} $.

(iv) Among all jet pairs, the pair with the maximal pseudorapidity separation: $ \Delta \eta(j, j')_{\rm{max}} $ and the corresponding dijet mass $ m(j + j')_{ {\rm{max}} \Delta \eta} $.

(v) The invariant mass of the first two leading jets and the jet multiplicity: $ m(j_{1}+j_{2}) $ and $ N(j) $.

(vi) Angular separations between key objects: $ \Delta R(j_1,j_2) $, $ \Delta\eta(j_1,j_2) $, $ \Delta R(e,\mu) $, $ \Delta R(\mu,j_1) $, $ \Delta R(\mu,j_2) $, $ \Delta R(e,j_1) $, and $ \Delta R(e,j_2) $.

(vii) Transverse masses of leptons and $ \not{E}_T $: $ m_{\rm{T}}(\mu + \not{E}_T) $, $ m_{\rm{T}}(e + \not{E}_T) $, and $ m_{\rm{T}}(\mu + e + \not{E}_T) $. Here, the transverse mass $ m_{\mathrm T} $ of a visible system and the missing transverse momentum is defined as $ m_{\mathrm T}^2 = \big(E_{\mathrm T}^{\rm{vis}} + \not{E}_T \big)^2 - \bigl|\vec{p}_{\mathrm T}^{\,\rm vis} + \vec{p}_{\mathrm T}^{\,\rm miss}\bigr|^2 , $ where $ E_{\mathrm T}^{\rm{vis}} = \sqrt{|\vec{p}_{\mathrm T}^{\,\rm vis}|^2 + (m^{\rm{vis}})^2} $ is the transverse energy of the visible system and $ \not{E}_T = |\vec{p}_{\mathrm T}^{\,\rm miss}| $ assumes a massless invisible system.

Observables (iv)–(vi) exploit the VBF-like topology: two forward tagging jets, large $ \Delta\eta_{jj} $, and high $ m_{jj} $ with reduced central hadronic activity, while (iii) captures s-channel topologies where jets tend to be closer in phase space. This division follows standard LHC practices for VBF/di-jet tagging. Observables (vii) are related to the kinematics of the decay $ a\to W^+W^-\to e^\pm\mu^\mp\nu \bar{\nu} $ in the $ e^\pm\mu^\mp+ \not{E}_T $ final state.

We implement the MVA using a BDT via the TMVA package integrated in ROOT, adopting the default BDT configuration for classification [97]. This setup is used to maximize background rejection while retaining high signal efficiency. Figure 3 shows the BDT responses for two benchmark masses, $ m_a=750\;\ {\rm{GeV}} $ and $ m_a=2000\;\ {\rm{GeV}} $. As $ m_a $ increases, the underlying kinematics becomes more boosted and the separation between signal and background events improves, as shown by the shift of the signal (background) distribution towards higher (lower) BDT scores.

Figure 3.  (color online) BDT response distributions for signal (shaded) and total background (outline) at the HL-LHC with $ \sqrt{s}= $$ 14\;\ {\rm{TeV}} $, for $ m_a=750\;\ {\rm{GeV}} $ (left) and $ m_a=2000\;\ {\rm{GeV}} $ (right).

The BDT score threshold is optimized independently at each $ m_a $ by maximizing the statistical significance

where $ N_s $ and $ N_b $ denote the expected signal and total background yields after preselection and the BDT cut, respectively [98].

For completeness, Appendix B presents representative input-observable distributions after preselection for benchmark $ m_a=750\;\ {\rm{GeV}} $, while Appendix C shows BDT responses across the full $ m_a $ range for various $ m_a $ cases. Appendix D reports the preselection and BDT efficiencies for signals and backgrounds at $ \sqrt{s}=14\;\ {\rm{TeV}} $. In that table, a dash ("−") indicates that the surviving yield of a given background is negligible at $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $. Despite their large inclusive cross sections, the $ W^\pm(\to\ell^\pm\nu)jj $ and $ Z(\to\ell^+\ell^-)jj $ background events are almost completely removed by the opposite-sign $ e\mu $ requirement and multivariate selection, so they are omitted from the appendix table for brevity.

After the preselection, we employ a multivariate analysis (MVA) to further separate signal from background. The input observables are as follows:

(i) Object kinematics: energy, transverse momentum, pseudorapidity, and azimuth of the final-state objects, $ E(O) $, $ p_T(O) $, $ \eta(O) $, $ \phi(O) $ with $ O=\mu,\,e,\,j_1,\,j_2 $.

(ii) Missing transverse momentum (magnitude and azimuth): $ \not{E}_T $ and $ \phi( \not{E}_T) $.

(iii) Among all jet pairs, the pair with the minimal angular separation: $ \Delta R (j, j')_{\rm{min}} $ and the corresponding dijet mass $ m(j+j')_{ {\rm{min}} \Delta R} $, where $ \Delta R = \sqrt{(\Delta \eta)^2 + (\Delta \varphi)^2} $.

(iv) Among all jet pairs, the pair with the maximal pseudorapidity separation: $ \Delta \eta(j, j')_{\rm{max}} $ and the corresponding dijet mass $ m(j + j')_{ {\rm{max}} \Delta \eta} $.

(v) The invariant mass of the first two leading jets and the jet multiplicity: $ m(j_{1}+j_{2}) $ and $ N(j) $.

(vi) Angular separations between key objects: $ \Delta R(j_1,j_2) $, $ \Delta\eta(j_1,j_2) $, $ \Delta R(e,\mu) $, $ \Delta R(\mu,j_1) $, $ \Delta R(\mu,j_2) $, $ \Delta R(e,j_1) $, and $ \Delta R(e,j_2) $.

(vii) Transverse masses of leptons and $ \not{E}_T $: $ m_{\rm{T}}(\mu + \not{E}_T) $, $ m_{\rm{T}}(e + \not{E}_T) $, and $ m_{\rm{T}}(\mu + e + \not{E}_T) $. Here, the transverse mass $ m_{\mathrm T} $ of a visible system and the missing transverse momentum is defined as $ m_{\mathrm T}^2 = \big(E_{\mathrm T}^{\rm{vis}} + \not{E}_T \big)^2 - \bigl|\vec{p}_{\mathrm T}^{\,\rm vis} + \vec{p}_{\mathrm T}^{\,\rm miss}\bigr|^2 , $ where $ E_{\mathrm T}^{\rm{vis}} = \sqrt{|\vec{p}_{\mathrm T}^{\,\rm vis}|^2 + (m^{\rm{vis}})^2} $ is the transverse energy of the visible system and $ \not{E}_T = |\vec{p}_{\mathrm T}^{\,\rm miss}| $ assumes a massless invisible system.

Observables (iv)–(vi) exploit the VBF-like topology: two forward tagging jets, large $ \Delta\eta_{jj} $, and high $ m_{jj} $ with reduced central hadronic activity, while (iii) captures s-channel topologies where jets tend to be closer in phase space. This division follows standard LHC practices for VBF/di-jet tagging. Observables (vii) are related to the kinematics of the decay $ a\to W^+W^-\to e^\pm\mu^\mp\nu \bar{\nu} $ in the $ e^\pm\mu^\mp+ \not{E}_T $ final state.

We implement the MVA using a BDT via the TMVA package integrated in ROOT, adopting the default BDT configuration for classification [97]. This setup is used to maximize background rejection while retaining high signal efficiency. Figure 3 shows the BDT responses for two benchmark masses, $ m_a=750\;\ {\rm{GeV}} $ and $ m_a=2000\;\ {\rm{GeV}} $. As $ m_a $ increases, the underlying kinematics becomes more boosted and the separation between signal and background events improves, as shown by the shift of the signal (background) distribution towards higher (lower) BDT scores.

Figure 3.  (color online) BDT response distributions for signal (shaded) and total background (outline) at the HL-LHC with $ \sqrt{s}= $$ 14\;\ {\rm{TeV}} $, for $ m_a=750\;\ {\rm{GeV}} $ (left) and $ m_a=2000\;\ {\rm{GeV}} $ (right).

The BDT score threshold is optimized independently at each $ m_a $ by maximizing the statistical significance

where $ N_s $ and $ N_b $ denote the expected signal and total background yields after preselection and the BDT cut, respectively [98].

For completeness, Appendix B presents representative input-observable distributions after preselection for benchmark $ m_a=750\;\ {\rm{GeV}} $, while Appendix C shows BDT responses across the full $ m_a $ range for various $ m_a $ cases. Appendix D reports the preselection and BDT efficiencies for signals and backgrounds at $ \sqrt{s}=14\;\ {\rm{TeV}} $. In that table, a dash ("−") indicates that the surviving yield of a given background is negligible at $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $. Despite their large inclusive cross sections, the $ W^\pm(\to\ell^\pm\nu)jj $ and $ Z(\to\ell^+\ell^-)jj $ background events are almost completely removed by the opposite-sign $ e\mu $ requirement and multivariate selection, so they are omitted from the appendix table for brevity.

Building on the strategy above, we present the projected discovery sensitivities at the HL-LHC with $ \sqrt{s}= 14\;\ {\rm{TeV}} $ and $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $. We show both the discovery reach on the coupling $ g_{aWW} $ and corresponding model-independent reach on $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $ across $ m_a\in[170,4000]\;\ {\rm{GeV}} $. Statistical significances are evaluated with Eq. (13). All results assume the photophobic limit $ g_{a\gamma\gamma}=0 $, as defined in Ref. [77]; the BDT score threshold is optimized independently at each mass point (see Sec. V.B). Unless stated otherwise, the quoted sensitivities are statistical only; the impact of systematics can be estimated in reinterpretations by introducing background-normalization nuisances in a profile-likelihood fit following Ref. [98].

Figure 4 reports the corresponding $ 2\sigma $ and $ 5\sigma $ sensitivities on the model-independent fiducial production cross section $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $, facilitating reinterpretations in alternative theory models. Near the threshold ($ m_a\sim170 $–$ 200\;\ {\rm{GeV}} $), the sensitivity varies quickly as the $ W^+W^- $ mode turns on and the event kinematics evolve; meanwhile, in the larger $ m_a $ region, the signal–background separation improves (cf. Sec. V.B), which offsets the falling production rate and allows competitiveness up to multi-TeV masses.

Figure 4.  (color online) Projected $ 2\sigma $ and $ 5\sigma $ discovery sensitivities on the model-independent fiducial production cross-section $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $ versus $ m_a $ at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $, $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $). The projection assumes $ g_{a\gamma\gamma}=0 $ and uses per-mass BDT optimization. Red and blue curves correspond to $ 2\sigma $ and $ 5\sigma $, respectively.

Figure 5 summarizes the $ 2\sigma $ and $ 5\sigma $ discovery reaches on $ g_{aWW} $ for heavy photophobic ALPs. At $ m_a=170\;\ {\rm{GeV}} $, the $ 2\sigma\,(5\sigma) $ discovery threshold is $g_{aWW}= $ $ 0.60\,(0.76)\;\ {\rm{TeV}}^{-1} $. As the $ a\to W^+W^- $ channel opens and production cross-section of $ pp\to jj\,a(\to W^\pm W^\mp) $ grows rapidly between $ 170 $ and $ 200\;\ {\rm{GeV}} $ (cf. Fig. 2), the required coupling decreases; over $ 260 \le m_a \le 1500\;\ {\rm{GeV}} $, it remains near $ 0.61\,(0.76)\;\ {\rm{TeV}}^{-1} $. At very high masses, the reach weakens due to falling PDFs and reduced acceptance, with the $ 2\sigma\,(5\sigma) $ discovery threshold $ g_{aWW}= 2.75\,(3.48)\;\ {\rm{TeV}}^{-1} $ at $ m_a=4\;\ {\rm{TeV}} $.

Figure 5.  (color online) Projected discovery reach on $ g_{aWW} $ for photophobic ALPs as a function of $ m_a $ (170–4000 GeV) at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $, $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $). Solid (dashed) black curves denote $ 2\sigma $ ($ 5\sigma $) significances. The projection assumes $ g_{a\gamma\gamma}=0 $ and uses per–mass-point BDT optimization; significances are statistical only. For comparison, we overlay 95% CL limits recast from triboson $ pp \to Z (\to \nu\bar{\nu}) \gamma \gamma $ at 8 TeV with 20.3 fb⁻¹ [77] and from Run-2 analyses at 13 TeV: $ pp \to W^\pm W^\pm W^\mp $ with 35.9 fb^-1 and $ pp \to Z(\to\ell^+\ell^-)\gamma jj $ with 140 fb^-1 [81]. Also shown are HL-LHC ($ 14\;\ {\rm{TeV}} $, $ 3\;\ {\rm{ab}}^{-1} $) $ 2\sigma $ projections for $ pp \to W^\pm a(\to W^+W^-) $ [84] and $ pp \to jj\,a(\to Z\gamma) $ [82].

Figure 5 also displays previous constraints/projections for context. The magenta curve corresponds to the 95% C.L. limit obtained by reinterpreting the Run-1 triboson channel $ pp \to Z(\to\nu\bar{\nu})\gamma\gamma $ at $ \sqrt{s}=8 $ TeV with $ 20.3\; {\rm{fb}}^{-1} $ [77]. The blue and red curves show Run-2 recasts at $ \sqrt{s}=13 $ TeV: $ pp \to W^\pm W^\pm W^\mp $ with $ 35.9\; {\rm{fb}}^{-1} $ and $ pp \to Z(\to\ell^+\ell^-)\gamma jj $ with $ 140\; {\rm{fb}}^{-1} $, respectively, as reported in Ref. [81]. The green curve is an HL-LHC $ 2\sigma $ projection for the tri-W topology $ pp \to W^\pm a(\to W^+W^-) $ using a same-sign dimuon plus hadronic-W selection and an MVA [84]. The yellow curve shows the HL-LHC $ 2\sigma $ projection for $ pp \to jj\,a(\to Z\gamma) $, which exploits a clean $ Z(\to\ell^+\ell^-) $+γ resonance with both the VBF-like and s-channel productions and detector-level MVA [82].

Relative to the tri-W strategy (green), our dedicated $ jj\, a(\to W^+ W^-) $ analysis achieves a stronger reach (lower $ g_{aWW} $) over all mass ranges, with a particularly visible improvement in the high mass region ($ m_a \gtrsim $ 700 GeV), where the VBF-like kinematics sharpen signal–background separation. Compared with that of the $ jj\,a(\to Z\gamma) $ projection (yellow), our reach is weaker across the mass range. The reason is clear: the $ jj\,a(\to Z\gamma) $ final state benefits from a narrow mass peak, excellent photon/lepton resolution, and lower irreducible backgrounds [82]; by contrast, although the branching ratio of $ a \to W^+ W^- $ is larger than that of $ a \to Z \gamma $, the dileptonic-$ WW $ channel involves missing energy, larger SM backgrounds ($ WW $, $ t\bar t $), and no fully reconstructible resonance, which intrinsically reduces sensitivity. Despite this, the $ jj\, a(\to W^+ W^-) $ topology remains important and complementary. First, it directly probes the $ aWW $ vertex via $ a \to W^+W^- $, while $ a \to Z\gamma $ tests the $ aZ\gamma $ interaction; a joint observation/non-observation across channels would enable a nontrivial check of the photophobic relations $ g_{aZ\gamma}=t_\theta\,g_{aWW} $ and $ g_{aZZ}=(1-t_\theta^2)\,g_{aWW} $ and thus of the underlying electroweak structure [77, 81]. Second, systematics and backgrounds are largely orthogonal between $ jj\, a(\to W^+ W^-) $ and $ jj\,a(\to Z\gamma) $, providing robustness against analysis-specific effects and a richer basis for global fits. Third, we provide model-independent discovery thresholds for $ \sigma(pp \to jj\,a)\times{\rm{Br}}(a \to W^+W^-) $, which enable straightforward reinterpretation for other neutral resonances decaying to $ WW $ (including scenarios where $ aZ\gamma $ is suppressed but $ aWW $ remains sizable) [84]. Taken together, these considerations warrant a dedicated $ jj\, a(\to W^+ W^-) $ search parallel to $ Z\gamma $ analyses.

To gauge the impact of background systematics, we take the statistical-only reach as a baseline and estimate a background-normalization uncertainty using an Asimov profile-likelihood treatment (Appendix A.3) [98, 99]. For the benchmark $ m_a=750\;\ {\rm{GeV}} $, Table A2 indicates that the required $ \sigma\times{\rm{Br}} $ corresponding to 2σ (5σ) significance increases by $ \sim1.82\; (1.86) $ for $ \delta_b = 5\% $ and by $ \sim 3.45\; (3.58) $ for $ \delta_b$ = 10% compared with those in the statistical-only case. In realistic analyses, dominant backgrounds are typically constrained with dedicated control regions, which can reduce the effective $ \delta_b $; our $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ thresholds therefore facilitate straightforward rescaling under alternative systematic assumptions.

Regarding EFT validity and interpretation, our $ g_{aWW} $ reach is derived in the minimal photophobic ALP EFT truncated at dimension-5. At large masses/couplings, the relevant partonic energies in the selected $ jjWW $ events can approach the regime where perturbative unitarity of this truncation is lost; therefore, we regard the $ J=1 $ partial-wave bound from Ref. [88] (Appendix A.2) as an indication of where additional UV dynamics or unitarisation effects must enter. For clarity, we do not translate Eq. (A2) into a mass-dependent line in Fig. 5, as the relevant $ \sqrt{\hat{s}} $ is process- and selection-dependent and cannot be identified with $ m_a $; we use it only as a conservative validity guide (Appendix A.2). In this region, the mapping between an experimental rate and a single EFT parameter $ g_{aWW} $ becomes UV-sensitive (e.g., via higher-dimensional operators, resonances, or form factors); therefore, we emphasize our model-agnostic sensitivity in $ \sigma(pp\to jja)\times{\rm{Br}}(a\to WW) $ as the primary benchmark for reinterpretation.

Building on the strategy above, we present the projected discovery sensitivities at the HL-LHC with $ \sqrt{s}= 14\;\ {\rm{TeV}} $ and $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $. We show both the discovery reach on the coupling $ g_{aWW} $ and corresponding model-independent reach on $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $ across $ m_a\in[170,4000]\;\ {\rm{GeV}} $. Statistical significances are evaluated with Eq. (13). All results assume the photophobic limit $ g_{a\gamma\gamma}=0 $, as defined in Ref. [77]; the BDT score threshold is optimized independently at each mass point (see Sec. V.B). Unless stated otherwise, the quoted sensitivities are statistical only; the impact of systematics can be estimated in reinterpretations by introducing background-normalization nuisances in a profile-likelihood fit following Ref. [98].

Figure 4 reports the corresponding $ 2\sigma $ and $ 5\sigma $ sensitivities on the model-independent fiducial production cross section $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $, facilitating reinterpretations in alternative theory models. Near the threshold ($ m_a\sim170 $–$ 200\;\ {\rm{GeV}} $), the sensitivity varies quickly as the $ W^+W^- $ mode turns on and the event kinematics evolve; meanwhile, in the larger $ m_a $ region, the signal–background separation improves (cf. Sec. V.B), which offsets the falling production rate and allows competitiveness up to multi-TeV masses.

Figure 4.  (color online) Projected $ 2\sigma $ and $ 5\sigma $ discovery sensitivities on the model-independent fiducial production cross-section $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $ versus $ m_a $ at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $, $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $). The projection assumes $ g_{a\gamma\gamma}=0 $ and uses per-mass BDT optimization. Red and blue curves correspond to $ 2\sigma $ and $ 5\sigma $, respectively.

Figure 5 summarizes the $ 2\sigma $ and $ 5\sigma $ discovery reaches on $ g_{aWW} $ for heavy photophobic ALPs. At $ m_a=170\;\ {\rm{GeV}} $, the $ 2\sigma\,(5\sigma) $ discovery threshold is $g_{aWW}= $ $ 0.60\,(0.76)\;\ {\rm{TeV}}^{-1} $. As the $ a\to W^+W^- $ channel opens and production cross-section of $ pp\to jj\,a(\to W^\pm W^\mp) $ grows rapidly between $ 170 $ and $ 200\;\ {\rm{GeV}} $ (cf. Fig. 2), the required coupling decreases; over $ 260 \le m_a \le 1500\;\ {\rm{GeV}} $, it remains near $ 0.61\,(0.76)\;\ {\rm{TeV}}^{-1} $. At very high masses, the reach weakens due to falling PDFs and reduced acceptance, with the $ 2\sigma\,(5\sigma) $ discovery threshold $ g_{aWW}= 2.75\,(3.48)\;\ {\rm{TeV}}^{-1} $ at $ m_a=4\;\ {\rm{TeV}} $.

Figure 5.  (color online) Projected discovery reach on $ g_{aWW} $ for photophobic ALPs as a function of $ m_a $ (170–4000 GeV) at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $, $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $). Solid (dashed) black curves denote $ 2\sigma $ ($ 5\sigma $) significances. The projection assumes $ g_{a\gamma\gamma}=0 $ and uses per–mass-point BDT optimization; significances are statistical only. For comparison, we overlay 95% CL limits recast from triboson $ pp \to Z (\to \nu\bar{\nu}) \gamma \gamma $ at 8 TeV with 20.3 fb⁻¹ [77] and from Run-2 analyses at 13 TeV: $ pp \to W^\pm W^\pm W^\mp $ with 35.9 fb^-1 and $ pp \to Z(\to\ell^+\ell^-)\gamma jj $ with 140 fb^-1 [81]. Also shown are HL-LHC ($ 14\;\ {\rm{TeV}} $, $ 3\;\ {\rm{ab}}^{-1} $) $ 2\sigma $ projections for $ pp \to W^\pm a(\to W^+W^-) $ [84] and $ pp \to jj\,a(\to Z\gamma) $ [82].

Figure 5 also displays previous constraints/projections for context. The magenta curve corresponds to the 95% C.L. limit obtained by reinterpreting the Run-1 triboson channel $ pp \to Z(\to\nu\bar{\nu})\gamma\gamma $ at $ \sqrt{s}=8 $ TeV with $ 20.3\; {\rm{fb}}^{-1} $ [77]. The blue and red curves show Run-2 recasts at $ \sqrt{s}=13 $ TeV: $ pp \to W^\pm W^\pm W^\mp $ with $ 35.9\; {\rm{fb}}^{-1} $ and $ pp \to Z(\to\ell^+\ell^-)\gamma jj $ with $ 140\; {\rm{fb}}^{-1} $, respectively, as reported in Ref. [81]. The green curve is an HL-LHC $ 2\sigma $ projection for the tri-W topology $ pp \to W^\pm a(\to W^+W^-) $ using a same-sign dimuon plus hadronic-W selection and an MVA [84]. The yellow curve shows the HL-LHC $ 2\sigma $ projection for $ pp \to jj\,a(\to Z\gamma) $, which exploits a clean $ Z(\to\ell^+\ell^-) $+γ resonance with both the VBF-like and s-channel productions and detector-level MVA [82].

Relative to the tri-W strategy (green), our dedicated $ jj\, a(\to W^+ W^-) $ analysis achieves a stronger reach (lower $ g_{aWW} $) over all mass ranges, with a particularly visible improvement in the high mass region ($ m_a \gtrsim $ 700 GeV), where the VBF-like kinematics sharpen signal–background separation. Compared with that of the $ jj\,a(\to Z\gamma) $ projection (yellow), our reach is weaker across the mass range. The reason is clear: the $ jj\,a(\to Z\gamma) $ final state benefits from a narrow mass peak, excellent photon/lepton resolution, and lower irreducible backgrounds [82]; by contrast, although the branching ratio of $ a \to W^+ W^- $ is larger than that of $ a \to Z \gamma $, the dileptonic-$ WW $ channel involves missing energy, larger SM backgrounds ($ WW $, $ t\bar t $), and no fully reconstructible resonance, which intrinsically reduces sensitivity. Despite this, the $ jj\, a(\to W^+ W^-) $ topology remains important and complementary. First, it directly probes the $ aWW $ vertex via $ a \to W^+W^- $, while $ a \to Z\gamma $ tests the $ aZ\gamma $ interaction; a joint observation/non-observation across channels would enable a nontrivial check of the photophobic relations $ g_{aZ\gamma}=t_\theta\,g_{aWW} $ and $ g_{aZZ}=(1-t_\theta^2)\,g_{aWW} $ and thus of the underlying electroweak structure [77, 81]. Second, systematics and backgrounds are largely orthogonal between $ jj\, a(\to W^+ W^-) $ and $ jj\,a(\to Z\gamma) $, providing robustness against analysis-specific effects and a richer basis for global fits. Third, we provide model-independent discovery thresholds for $ \sigma(pp \to jj\,a)\times{\rm{Br}}(a \to W^+W^-) $, which enable straightforward reinterpretation for other neutral resonances decaying to $ WW $ (including scenarios where $ aZ\gamma $ is suppressed but $ aWW $ remains sizable) [84]. Taken together, these considerations warrant a dedicated $ jj\, a(\to W^+ W^-) $ search parallel to $ Z\gamma $ analyses.

To gauge the impact of background systematics, we take the statistical-only reach as a baseline and estimate a background-normalization uncertainty using an Asimov profile-likelihood treatment (Appendix A.3) [98, 99]. For the benchmark $ m_a=750\;\ {\rm{GeV}} $, Table A2 indicates that the required $ \sigma\times{\rm{Br}} $ corresponding to 2σ (5σ) significance increases by $ \sim1.82\; (1.86) $ for $ \delta_b = 5\% $ and by $ \sim 3.45\; (3.58) $ for $ \delta_b$ = 10% compared with those in the statistical-only case. In realistic analyses, dominant backgrounds are typically constrained with dedicated control regions, which can reduce the effective $ \delta_b $; our $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ thresholds therefore facilitate straightforward rescaling under alternative systematic assumptions.

Regarding EFT validity and interpretation, our $ g_{aWW} $ reach is derived in the minimal photophobic ALP EFT truncated at dimension-5. At large masses/couplings, the relevant partonic energies in the selected $ jjWW $ events can approach the regime where perturbative unitarity of this truncation is lost; therefore, we regard the $ J=1 $ partial-wave bound from Ref. [88] (Appendix A.2) as an indication of where additional UV dynamics or unitarisation effects must enter. For clarity, we do not translate Eq. (A2) into a mass-dependent line in Fig. 5, as the relevant $ \sqrt{\hat{s}} $ is process- and selection-dependent and cannot be identified with $ m_a $; we use it only as a conservative validity guide (Appendix A.2). In this region, the mapping between an experimental rate and a single EFT parameter $ g_{aWW} $ becomes UV-sensitive (e.g., via higher-dimensional operators, resonances, or form factors); therefore, we emphasize our model-agnostic sensitivity in $ \sigma(pp\to jja)\times{\rm{Br}}(a\to WW) $ as the primary benchmark for reinterpretation.

Building on the strategy above, we present the projected discovery sensitivities at the HL-LHC with $ \sqrt{s}= 14\;\ {\rm{TeV}} $ and $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $. We show both the discovery reach on the coupling $ g_{aWW} $ and corresponding model-independent reach on $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $ across $ m_a\in[170,4000]\;\ {\rm{GeV}} $. Statistical significances are evaluated with Eq. (13). All results assume the photophobic limit $ g_{a\gamma\gamma}=0 $, as defined in Ref. [77]; the BDT score threshold is optimized independently at each mass point (see Sec. V.B). Unless stated otherwise, the quoted sensitivities are statistical only; the impact of systematics can be estimated in reinterpretations by introducing background-normalization nuisances in a profile-likelihood fit following Ref. [98].

Figure 4 reports the corresponding $ 2\sigma $ and $ 5\sigma $ sensitivities on the model-independent fiducial production cross section $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $, facilitating reinterpretations in alternative theory models. Near the threshold ($ m_a\sim170 $–$ 200\;\ {\rm{GeV}} $), the sensitivity varies quickly as the $ W^+W^- $ mode turns on and the event kinematics evolve; meanwhile, in the larger $ m_a $ region, the signal–background separation improves (cf. Sec. V.B), which offsets the falling production rate and allows competitiveness up to multi-TeV masses.

Figure 4.  (color online) Projected $ 2\sigma $ and $ 5\sigma $ discovery sensitivities on the model-independent fiducial production cross-section $ \sigma(pp\to jj\,a)\times {\rm{Br}}(a\to W^+W^-) $ versus $ m_a $ at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $, $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $). The projection assumes $ g_{a\gamma\gamma}=0 $ and uses per-mass BDT optimization. Red and blue curves correspond to $ 2\sigma $ and $ 5\sigma $, respectively.

Figure 5 summarizes the $ 2\sigma $ and $ 5\sigma $ discovery reaches on $ g_{aWW} $ for heavy photophobic ALPs. At $ m_a=170\;\ {\rm{GeV}} $, the $ 2\sigma\,(5\sigma) $ discovery threshold is $g_{aWW}= $ $ 0.60\,(0.76)\;\ {\rm{TeV}}^{-1} $. As the $ a\to W^+W^- $ channel opens and production cross-section of $ pp\to jj\,a(\to W^\pm W^\mp) $ grows rapidly between $ 170 $ and $ 200\;\ {\rm{GeV}} $ (cf. Fig. 2), the required coupling decreases; over $ 260 \le m_a \le 1500\;\ {\rm{GeV}} $, it remains near $ 0.61\,(0.76)\;\ {\rm{TeV}}^{-1} $. At very high masses, the reach weakens due to falling PDFs and reduced acceptance, with the $ 2\sigma\,(5\sigma) $ discovery threshold $ g_{aWW}= 2.75\,(3.48)\;\ {\rm{TeV}}^{-1} $ at $ m_a=4\;\ {\rm{TeV}} $.

Figure 5.  (color online) Projected discovery reach on $ g_{aWW} $ for photophobic ALPs as a function of $ m_a $ (170–4000 GeV) at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $, $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $). Solid (dashed) black curves denote $ 2\sigma $ ($ 5\sigma $) significances. The projection assumes $ g_{a\gamma\gamma}=0 $ and uses per–mass-point BDT optimization; significances are statistical only. For comparison, we overlay 95% CL limits recast from triboson $ pp \to Z (\to \nu\bar{\nu}) \gamma \gamma $ at 8 TeV with 20.3 fb⁻¹ [77] and from Run-2 analyses at 13 TeV: $ pp \to W^\pm W^\pm W^\mp $ with 35.9 fb^-1 and $ pp \to Z(\to\ell^+\ell^-)\gamma jj $ with 140 fb^-1 [81]. Also shown are HL-LHC ($ 14\;\ {\rm{TeV}} $, $ 3\;\ {\rm{ab}}^{-1} $) $ 2\sigma $ projections for $ pp \to W^\pm a(\to W^+W^-) $ [84] and $ pp \to jj\,a(\to Z\gamma) $ [82].

Figure 5 also displays previous constraints/projections for context. The magenta curve corresponds to the 95% C.L. limit obtained by reinterpreting the Run-1 triboson channel $ pp \to Z(\to\nu\bar{\nu})\gamma\gamma $ at $ \sqrt{s}=8 $ TeV with $ 20.3\; {\rm{fb}}^{-1} $ [77]. The blue and red curves show Run-2 recasts at $ \sqrt{s}=13 $ TeV: $ pp \to W^\pm W^\pm W^\mp $ with $ 35.9\; {\rm{fb}}^{-1} $ and $ pp \to Z(\to\ell^+\ell^-)\gamma jj $ with $ 140\; {\rm{fb}}^{-1} $, respectively, as reported in Ref. [81]. The green curve is an HL-LHC $ 2\sigma $ projection for the tri-W topology $ pp \to W^\pm a(\to W^+W^-) $ using a same-sign dimuon plus hadronic-W selection and an MVA [84]. The yellow curve shows the HL-LHC $ 2\sigma $ projection for $ pp \to jj\,a(\to Z\gamma) $, which exploits a clean $ Z(\to\ell^+\ell^-) $+γ resonance with both the VBF-like and s-channel productions and detector-level MVA [82].

Relative to the tri-W strategy (green), our dedicated $ jj\, a(\to W^+ W^-) $ analysis achieves a stronger reach (lower $ g_{aWW} $) over all mass ranges, with a particularly visible improvement in the high mass region ($ m_a \gtrsim $ 700 GeV), where the VBF-like kinematics sharpen signal–background separation. Compared with that of the $ jj\,a(\to Z\gamma) $ projection (yellow), our reach is weaker across the mass range. The reason is clear: the $ jj\,a(\to Z\gamma) $ final state benefits from a narrow mass peak, excellent photon/lepton resolution, and lower irreducible backgrounds [82]; by contrast, although the branching ratio of $ a \to W^+ W^- $ is larger than that of $ a \to Z \gamma $, the dileptonic-$ WW $ channel involves missing energy, larger SM backgrounds ($ WW $, $ t\bar t $), and no fully reconstructible resonance, which intrinsically reduces sensitivity. Despite this, the $ jj\, a(\to W^+ W^-) $ topology remains important and complementary. First, it directly probes the $ aWW $ vertex via $ a \to W^+W^- $, while $ a \to Z\gamma $ tests the $ aZ\gamma $ interaction; a joint observation/non-observation across channels would enable a nontrivial check of the photophobic relations $ g_{aZ\gamma}=t_\theta\,g_{aWW} $ and $ g_{aZZ}=(1-t_\theta^2)\,g_{aWW} $ and thus of the underlying electroweak structure [77, 81]. Second, systematics and backgrounds are largely orthogonal between $ jj\, a(\to W^+ W^-) $ and $ jj\,a(\to Z\gamma) $, providing robustness against analysis-specific effects and a richer basis for global fits. Third, we provide model-independent discovery thresholds for $ \sigma(pp \to jj\,a)\times{\rm{Br}}(a \to W^+W^-) $, which enable straightforward reinterpretation for other neutral resonances decaying to $ WW $ (including scenarios where $ aZ\gamma $ is suppressed but $ aWW $ remains sizable) [84]. Taken together, these considerations warrant a dedicated $ jj\, a(\to W^+ W^-) $ search parallel to $ Z\gamma $ analyses.

To gauge the impact of background systematics, we take the statistical-only reach as a baseline and estimate a background-normalization uncertainty using an Asimov profile-likelihood treatment (Appendix A.3) [98, 99]. For the benchmark $ m_a=750\;\ {\rm{GeV}} $, Table A2 indicates that the required $ \sigma\times{\rm{Br}} $ corresponding to 2σ (5σ) significance increases by $ \sim1.82\; (1.86) $ for $ \delta_b = 5\% $ and by $ \sim 3.45\; (3.58) $ for $ \delta_b$ = 10% compared with those in the statistical-only case. In realistic analyses, dominant backgrounds are typically constrained with dedicated control regions, which can reduce the effective $ \delta_b $; our $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ thresholds therefore facilitate straightforward rescaling under alternative systematic assumptions.

Regarding EFT validity and interpretation, our $ g_{aWW} $ reach is derived in the minimal photophobic ALP EFT truncated at dimension-5. At large masses/couplings, the relevant partonic energies in the selected $ jjWW $ events can approach the regime where perturbative unitarity of this truncation is lost; therefore, we regard the $ J=1 $ partial-wave bound from Ref. [88] (Appendix A.2) as an indication of where additional UV dynamics or unitarisation effects must enter. For clarity, we do not translate Eq. (A2) into a mass-dependent line in Fig. 5, as the relevant $ \sqrt{\hat{s}} $ is process- and selection-dependent and cannot be identified with $ m_a $; we use it only as a conservative validity guide (Appendix A.2). In this region, the mapping between an experimental rate and a single EFT parameter $ g_{aWW} $ becomes UV-sensitive (e.g., via higher-dimensional operators, resonances, or form factors); therefore, we emphasize our model-agnostic sensitivity in $ \sigma(pp\to jja)\times{\rm{Br}}(a\to WW) $ as the primary benchmark for reinterpretation.

We assessed the discovery potential for heavy photophobic ALPs at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $), where the ALP–diphoton coupling vanishes, i.e., $ g_{a\gamma\gamma}=0 $, and interactions with electroweak gauge bosons [77]. In this limit, the ALP is produced in association with two jets, $ pp\to jj\,a $, via either s-channel electroweak-boson exchange or VBF-like topologies; for $ m_a\gtrsim 2\,m_W $, the decay channel $ a\to W^+W^- $ dominates. We target the opposite-sign, different-flavour dilepton final state $ e^\pm\mu^\mp+jj+ \not{E}_T $, which strongly suppresses Drell–Yan contamination and multijet backgrounds, and we adopt the HL-LHC baseline of $ \sqrt{s}=14\;\ {\rm{TeV}} $ and $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $ [83].

Events are generated with MadGraph5_aMC@NLO, showered and hadronised with PYTHIA 8.3, and passed through a Delphes 3 fast detector simulation [91, 94, 95]. The analysis applies a physics-motivated preselection for the $ e\mu+jj+ \not{E}_T $ topology and then a multivariate classification based on a BDT implemented in TMVA [97]. The BDT input features capture both VBF-like jet–gap kinematics (large $ m_{jj} $ and $ \Delta\eta_{jj} $, reduced central activity) and s-channel configurations (small $ \Delta R_{jj} $ and reconstructed $ m_{jj} \sim m_V $), and the BDT score threshold is optimised independently at each $ m_a $.

Our results provide the projected discovery reaches at $ 2\sigma $ and $ 5\sigma $ statistical significances on the ALP–W coupling $ g_{aWW} $ and, in a model-independent form, on $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $, over the mass range $ m_a\in [170,4000]\;\ {\rm{GeV}} $. Specifically, for $ 260\leqslant m_a\leqslant 1500 \ {\rm{GeV}} $, the sensitivity remains approximately constant at approximately $ 0.61\,(0.76)\;\ {\rm{TeV}}^{-1} $ at the $ 2\sigma\,(5\sigma) $ level. We also report discovery sensitivities on the fiducial production cross sections $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ across $ m_a\in[170,4000]\;\ {\rm{GeV}} $, enabling straightforward reinterpretations in alternative theory models.

For orientation, we compare our $ p p \to j j a\,(\to W^+ W^-) $ projection with recent HL-LHC studies and Run-1/2 recasts: the tri-W strategy $ pp \to W^\pm a(\to W^+W^-) $ [84], $ jj\,a(\to Z\gamma) $ channel [82], and related recasts [77, 81]. In short, our reach improves on the tri-W strategy (notably at high $ m_a $) but is weaker than the clean $ Z\gamma $ resonance channel (see Fig. 5). The $ WWjj $ topology remains complementary: it directly probes $ aWW $, carries largely orthogonal systematics to $ Z\gamma $, and provides model-independent thresholds in $ \sigma(pp \to jj\,a)\times{\rm{Br}}(a \to W^+W^-) $; full details and a quantitative comparison are deferred to Sec. VI.

In summary, the $ e\mu+jj+ \not{E}_T $ topology at the HL-LHC offers competitive sensitivity to heavy photophobic ALPs, with robust performance in the multi-TeV regime owing to improved signal–background separation at large $ m_a $. The model-independent $ \sigma\times{\rm{Br}} $ projections are intended to facilitate recasts beyond the specific benchmark ALP model investigated in this study.

We assessed the discovery potential for heavy photophobic ALPs at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $), where the ALP–diphoton coupling vanishes, i.e., $ g_{a\gamma\gamma}=0 $, and interactions with electroweak gauge bosons [77]. In this limit, the ALP is produced in association with two jets, $ pp\to jj\,a $, via either s-channel electroweak-boson exchange or VBF-like topologies; for $ m_a\gtrsim 2\,m_W $, the decay channel $ a\to W^+W^- $ dominates. We target the opposite-sign, different-flavour dilepton final state $ e^\pm\mu^\mp+jj+ \not{E}_T $, which strongly suppresses Drell–Yan contamination and multijet backgrounds, and we adopt the HL-LHC baseline of $ \sqrt{s}=14\;\ {\rm{TeV}} $ and $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $ [83].

Events are generated with MadGraph5_aMC@NLO, showered and hadronised with PYTHIA 8.3, and passed through a Delphes 3 fast detector simulation [91, 94, 95]. The analysis applies a physics-motivated preselection for the $ e\mu+jj+ \not{E}_T $ topology and then a multivariate classification based on a BDT implemented in TMVA [97]. The BDT input features capture both VBF-like jet–gap kinematics (large $ m_{jj} $ and $ \Delta\eta_{jj} $, reduced central activity) and s-channel configurations (small $ \Delta R_{jj} $ and reconstructed $ m_{jj} \sim m_V $), and the BDT score threshold is optimised independently at each $ m_a $.

Our results provide the projected discovery reaches at $ 2\sigma $ and $ 5\sigma $ statistical significances on the ALP–W coupling $ g_{aWW} $ and, in a model-independent form, on $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $, over the mass range $ m_a\in [170,4000]\;\ {\rm{GeV}} $. Specifically, for $ 260\leqslant m_a\leqslant 1500 \ {\rm{GeV}} $, the sensitivity remains approximately constant at approximately $ 0.61\,(0.76)\;\ {\rm{TeV}}^{-1} $ at the $ 2\sigma\,(5\sigma) $ level. We also report discovery sensitivities on the fiducial production cross sections $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ across $ m_a\in[170,4000]\;\ {\rm{GeV}} $, enabling straightforward reinterpretations in alternative theory models.

For orientation, we compare our $ p p \to j j a\,(\to W^+ W^-) $ projection with recent HL-LHC studies and Run-1/2 recasts: the tri-W strategy $ pp \to W^\pm a(\to W^+W^-) $ [84], $ jj\,a(\to Z\gamma) $ channel [82], and related recasts [77, 81]. In short, our reach improves on the tri-W strategy (notably at high $ m_a $) but is weaker than the clean $ Z\gamma $ resonance channel (see Fig. 5). The $ WWjj $ topology remains complementary: it directly probes $ aWW $, carries largely orthogonal systematics to $ Z\gamma $, and provides model-independent thresholds in $ \sigma(pp \to jj\,a)\times{\rm{Br}}(a \to W^+W^-) $; full details and a quantitative comparison are deferred to Sec. VI.

In summary, the $ e\mu+jj+ \not{E}_T $ topology at the HL-LHC offers competitive sensitivity to heavy photophobic ALPs, with robust performance in the multi-TeV regime owing to improved signal–background separation at large $ m_a $. The model-independent $ \sigma\times{\rm{Br}} $ projections are intended to facilitate recasts beyond the specific benchmark ALP model investigated in this study.

We assessed the discovery potential for heavy photophobic ALPs at the HL-LHC ($ \sqrt{s}=14\;\ {\rm{TeV}} $), where the ALP–diphoton coupling vanishes, i.e., $ g_{a\gamma\gamma}=0 $, and interactions with electroweak gauge bosons [77]. In this limit, the ALP is produced in association with two jets, $ pp\to jj\,a $, via either s-channel electroweak-boson exchange or VBF-like topologies; for $ m_a\gtrsim 2\,m_W $, the decay channel $ a\to W^+W^- $ dominates. We target the opposite-sign, different-flavour dilepton final state $ e^\pm\mu^\mp+jj+ \not{E}_T $, which strongly suppresses Drell–Yan contamination and multijet backgrounds, and we adopt the HL-LHC baseline of $ \sqrt{s}=14\;\ {\rm{TeV}} $ and $ {\cal{L}}=3\;\ {\rm{ab}}^{-1} $ [83].

Events are generated with MadGraph5_aMC@NLO, showered and hadronised with PYTHIA 8.3, and passed through a Delphes 3 fast detector simulation [91, 94, 95]. The analysis applies a physics-motivated preselection for the $ e\mu+jj+ \not{E}_T $ topology and then a multivariate classification based on a BDT implemented in TMVA [97]. The BDT input features capture both VBF-like jet–gap kinematics (large $ m_{jj} $ and $ \Delta\eta_{jj} $, reduced central activity) and s-channel configurations (small $ \Delta R_{jj} $ and reconstructed $ m_{jj} \sim m_V $), and the BDT score threshold is optimised independently at each $ m_a $.

Our results provide the projected discovery reaches at $ 2\sigma $ and $ 5\sigma $ statistical significances on the ALP–W coupling $ g_{aWW} $ and, in a model-independent form, on $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $, over the mass range $ m_a\in [170,4000]\;\ {\rm{GeV}} $. Specifically, for $ 260\leqslant m_a\leqslant 1500 \ {\rm{GeV}} $, the sensitivity remains approximately constant at approximately $ 0.61\,(0.76)\;\ {\rm{TeV}}^{-1} $ at the $ 2\sigma\,(5\sigma) $ level. We also report discovery sensitivities on the fiducial production cross sections $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ across $ m_a\in[170,4000]\;\ {\rm{GeV}} $, enabling straightforward reinterpretations in alternative theory models.

For orientation, we compare our $ p p \to j j a\,(\to W^+ W^-) $ projection with recent HL-LHC studies and Run-1/2 recasts: the tri-W strategy $ pp \to W^\pm a(\to W^+W^-) $ [84], $ jj\,a(\to Z\gamma) $ channel [82], and related recasts [77, 81]. In short, our reach improves on the tri-W strategy (notably at high $ m_a $) but is weaker than the clean $ Z\gamma $ resonance channel (see Fig. 5). The $ WWjj $ topology remains complementary: it directly probes $ aWW $, carries largely orthogonal systematics to $ Z\gamma $, and provides model-independent thresholds in $ \sigma(pp \to jj\,a)\times{\rm{Br}}(a \to W^+W^-) $; full details and a quantitative comparison are deferred to Sec. VI.

In summary, the $ e\mu+jj+ \not{E}_T $ topology at the HL-LHC offers competitive sensitivity to heavy photophobic ALPs, with robust performance in the multi-TeV regime owing to improved signal–background separation at large $ m_a $. The model-independent $ \sigma\times{\rm{Br}} $ projections are intended to facilitate recasts beyond the specific benchmark ALP model investigated in this study.

In the strict photophobic limit, one imposes $ g_{a\gamma\gamma}=0 $ at the matching scale, which corresponds to the coupling relation $ c_{\widetilde{B}}/c_{\widetilde{W}}=-t_\theta^2 $ with $ t_\theta\equiv\tan\theta $. To quantify a controlled departure, we introduce a dimensionless parameter δ in Eq. (10), which directly measures the mismatch of the cancellation in the photon direction.

Using the standard relations of Eqs. (3) and (6) between the gauge-eigenstate operators and the mass-eigenstate couplings, one finds that Eq. (10) implies

Therefore, for fixed $ g_{aWW} $, the diphoton partial width scales as $ \Gamma(a\to\gamma\gamma)\propto |g_{a\gamma\gamma}|^2\propto\delta^2 $, and the effect of $ \delta\neq 0 $ is to open/enhance $ a\to\gamma\gamma $ and correspondingly dilute the dominant electroweak modes.

Figure A1 shows the branching ratios of the ALP into $ W^+W^- $, $ ZZ $, $ Z\gamma $, and $ \gamma\gamma $ as functions of $ m_a $ for representative values $ \delta=0.1,\,0.3,\,0.5 $ (with other assumptions as in the main text). Above the $ WW $ threshold, $ a\to W^+W^- $ remains dominant for modest departures from the strict photophobic relation, while increasing δ enhances the diphoton mode and induces a corresponding reduction in $ {\rm{Br}}(a\to W^+W^-) $. In the heavy-mass region of interest for our HL-LHC study, $ \delta\lesssim0.3 $ keeps $ {\rm{Br}}(a\to\gamma\gamma) $ at the percent level or below, implying only mild changes to our $ jj\,a(\to W^+W^-) $ sensitivity; for larger deviations (e.g. $ \delta\sim0.5 $), the $ {\rm{Br}}(a\to\gamma\gamma) $ reaches ~ 5% level and diphoton searches can become increasingly competitive.

Figure A1.  (color online) Branching ratios of the ALP into $ W^+W^- $, $ ZZ $, $ \gamma Z $, and $ \gamma\gamma $ as functions of the ALP mass $ m_a $ for representative deviations from the strict photophobic relation, $ c_{\widetilde{B}}/c_{\widetilde{W}}=-t_\theta^2+\delta $ (Eq. (10), with $ t_\theta\equiv\tan\theta $). Solid, dashed, and dotted curves denote $ \delta=0.1 $, $ 0.3 $, and $ 0.5 $, respectively. Increasing δ enhances the diphoton mode and correspondingly reduces the $ WW/ZZ/Z\gamma $ branching fractions, while $ a\to W^+W^- $ remains dominant above threshold in the mass range shown.

Finally, we emphasize that the above discussion concerns tree-level departures from the strict photophobic relation. Even if $ \delta=0 $ is imposed at a high scale, radiative effects induce rare decay modes, such as $ a\to\gamma\gamma $ (from W loops) and $ a\to f\bar f,\,gg $ (from renormalization group running and/or loop effects). As shown explicitly in Ref. [82], in the mass region relevant to this study, these radiatively induced channels remain below the percent level and thus only negligibly affect the dominant $ WW/ZZ/Z\gamma $ pattern.

In the strict photophobic limit, one imposes $ g_{a\gamma\gamma}=0 $ at the matching scale, which corresponds to the coupling relation $ c_{\widetilde{B}}/c_{\widetilde{W}}=-t_\theta^2 $ with $ t_\theta\equiv\tan\theta $. To quantify a controlled departure, we introduce a dimensionless parameter δ in Eq. (10), which directly measures the mismatch of the cancellation in the photon direction.

Using the standard relations of Eqs. (3) and (6) between the gauge-eigenstate operators and the mass-eigenstate couplings, one finds that Eq. (10) implies

Therefore, for fixed $ g_{aWW} $, the diphoton partial width scales as $ \Gamma(a\to\gamma\gamma)\propto |g_{a\gamma\gamma}|^2\propto\delta^2 $, and the effect of $ \delta\neq 0 $ is to open/enhance $ a\to\gamma\gamma $ and correspondingly dilute the dominant electroweak modes.

Figure A1 shows the branching ratios of the ALP into $ W^+W^- $, $ ZZ $, $ Z\gamma $, and $ \gamma\gamma $ as functions of $ m_a $ for representative values $ \delta=0.1,\,0.3,\,0.5 $ (with other assumptions as in the main text). Above the $ WW $ threshold, $ a\to W^+W^- $ remains dominant for modest departures from the strict photophobic relation, while increasing δ enhances the diphoton mode and induces a corresponding reduction in $ {\rm{Br}}(a\to W^+W^-) $. In the heavy-mass region of interest for our HL-LHC study, $ \delta\lesssim0.3 $ keeps $ {\rm{Br}}(a\to\gamma\gamma) $ at the percent level or below, implying only mild changes to our $ jj\,a(\to W^+W^-) $ sensitivity; for larger deviations (e.g. $ \delta\sim0.5 $), the $ {\rm{Br}}(a\to\gamma\gamma) $ reaches ~ 5% level and diphoton searches can become increasingly competitive.

Figure A1.  (color online) Branching ratios of the ALP into $ W^+W^- $, $ ZZ $, $ \gamma Z $, and $ \gamma\gamma $ as functions of the ALP mass $ m_a $ for representative deviations from the strict photophobic relation, $ c_{\widetilde{B}}/c_{\widetilde{W}}=-t_\theta^2+\delta $ (Eq. (10), with $ t_\theta\equiv\tan\theta $). Solid, dashed, and dotted curves denote $ \delta=0.1 $, $ 0.3 $, and $ 0.5 $, respectively. Increasing δ enhances the diphoton mode and correspondingly reduces the $ WW/ZZ/Z\gamma $ branching fractions, while $ a\to W^+W^- $ remains dominant above threshold in the mass range shown.

Finally, we emphasize that the above discussion concerns tree-level departures from the strict photophobic relation. Even if $ \delta=0 $ is imposed at a high scale, radiative effects induce rare decay modes, such as $ a\to\gamma\gamma $ (from W loops) and $ a\to f\bar f,\,gg $ (from renormalization group running and/or loop effects). As shown explicitly in Ref. [82], in the mass region relevant to this study, these radiatively induced channels remain below the percent level and thus only negligibly affect the dominant $ WW/ZZ/Z\gamma $ pattern.

In the strict photophobic limit, one imposes $ g_{a\gamma\gamma}=0 $ at the matching scale, which corresponds to the coupling relation $ c_{\widetilde{B}}/c_{\widetilde{W}}=-t_\theta^2 $ with $ t_\theta\equiv\tan\theta $. To quantify a controlled departure, we introduce a dimensionless parameter δ in Eq. (10), which directly measures the mismatch of the cancellation in the photon direction.

Using the standard relations of Eqs. (3) and (6) between the gauge-eigenstate operators and the mass-eigenstate couplings, one finds that Eq. (10) implies

Therefore, for fixed $ g_{aWW} $, the diphoton partial width scales as $ \Gamma(a\to\gamma\gamma)\propto |g_{a\gamma\gamma}|^2\propto\delta^2 $, and the effect of $ \delta\neq 0 $ is to open/enhance $ a\to\gamma\gamma $ and correspondingly dilute the dominant electroweak modes.

Figure A1 shows the branching ratios of the ALP into $ W^+W^- $, $ ZZ $, $ Z\gamma $, and $ \gamma\gamma $ as functions of $ m_a $ for representative values $ \delta=0.1,\,0.3,\,0.5 $ (with other assumptions as in the main text). Above the $ WW $ threshold, $ a\to W^+W^- $ remains dominant for modest departures from the strict photophobic relation, while increasing δ enhances the diphoton mode and induces a corresponding reduction in $ {\rm{Br}}(a\to W^+W^-) $. In the heavy-mass region of interest for our HL-LHC study, $ \delta\lesssim0.3 $ keeps $ {\rm{Br}}(a\to\gamma\gamma) $ at the percent level or below, implying only mild changes to our $ jj\,a(\to W^+W^-) $ sensitivity; for larger deviations (e.g. $ \delta\sim0.5 $), the $ {\rm{Br}}(a\to\gamma\gamma) $ reaches ~ 5% level and diphoton searches can become increasingly competitive.

Figure A1.  (color online) Branching ratios of the ALP into $ W^+W^- $, $ ZZ $, $ \gamma Z $, and $ \gamma\gamma $ as functions of the ALP mass $ m_a $ for representative deviations from the strict photophobic relation, $ c_{\widetilde{B}}/c_{\widetilde{W}}=-t_\theta^2+\delta $ (Eq. (10), with $ t_\theta\equiv\tan\theta $). Solid, dashed, and dotted curves denote $ \delta=0.1 $, $ 0.3 $, and $ 0.5 $, respectively. Increasing δ enhances the diphoton mode and correspondingly reduces the $ WW/ZZ/Z\gamma $ branching fractions, while $ a\to W^+W^- $ remains dominant above threshold in the mass range shown.

Finally, we emphasize that the above discussion concerns tree-level departures from the strict photophobic relation. Even if $ \delta=0 $ is imposed at a high scale, radiative effects induce rare decay modes, such as $ a\to\gamma\gamma $ (from W loops) and $ a\to f\bar f,\,gg $ (from renormalization group running and/or loop effects). As shown explicitly in Ref. [82], in the mass region relevant to this study, these radiatively induced channels remain below the percent level and thus only negligibly affect the dominant $ WW/ZZ/Z\gamma $ pattern.

Partial-wave unitarity (conservative guide). Following Ref. [88], which derives partial-wave unitarity bounds in the gauge-invariant ALP EFT up to dimension-6, we use its dimension-5 truncation as a baseline guide. The strongest perturbative constraint arises from $ VV\to Va $ with one external ALP (dominant $ J=1 $ eigenvalue), primarily bounding the SU(2)$ _L $ coefficient $ c_{\widetilde{W}} $. Imposing $ |T_J|\le 1 $ in the high-energy limit yields

using $ g_{aWW}=4c_{\widetilde{W}}/f_a $. Because$ \sqrt{\hat{s}} $ is not directly observable event-by-event, we treat Eq. (A2) as a conservative guide for where the minimal dimension-5 truncation is expected to require UV unitarisation, rather than as a strict cut on our event selection.

Unitarity vs. EFT cutoff and UV sensitivity. The cutoff Λ controls the breakdown of the derivative expansion and is UV dependent. Parametrically one expects $ \Lambda \lesssim {4}/{|g_{aWW}|} $ for weakly-coupled UV completions, while NDA for strong dynamics suggests $ \Lambda \lesssim {16\pi}/{|g_{aWW}|} $ (cf. Eq. (11)). By contrast, Eq. (A2) constrains the perturbative description of specific high-energy $ 2\to2 $ channels in the truncated EFT and can be saturated already below the NDA cutoff estimate for the minimal dimension-5 truncation [88]. In that regime, additional UV ingredients (new states, higher-dimensional operators, or form factors) are expected to unitarise the amplitude, and the numerical bound can be relaxed once such effects are included (see, e.g., the modified limits when allowing dimension-6 contributions and cancellations in Ref. [88]). For collider applications, it is thus useful to treat the unitarity curve as a conservative validity guide and, if desired, to assess robustness by restricting to events below a chosen hard scale (e.g., an upper cut on a diboson invariant mass as in ALP-VBS studies [80]). Finally, because the high-energy rate in a UV-complete (unitarised) realization need not map one-to-one onto a single EFT coupling, we also provide the sensitivity in terms of $ \sigma(pp\to jja)\times {\rm{Br}}(a\to WW) $ for reinterpretation.

Total width estimate. Although $ m_a $ is not directly reconstructible due to the two neutrinos, the on-shell (factorised) signal modeling motivates verifying whether $ \Gamma_a/m_a $ remains small in the quoted sensitivity region. For dimension-5 gauge couplings, one expects $ \Gamma_a\propto g_{aWW}^2\,m_a^3 $ at large $ m_a $ (and hence $ \Gamma_a/m_a\propto g_{aWW}^2\,m_a^2 $). We compute the total width $ \Gamma_a $ in the strict photophobic benchmark with $\tt{MadGraph5\_aMC@NLO} $ and report $ \Gamma_a/m_a $ along the $ 2\sigma $ and $ 5\sigma $ reach curves in Table A1. We find $ \Gamma_a/m_a\lesssim$ 5% (8%) on the $ 2\sigma $ ($ 5\sigma $) curve for $ m_a\lesssim 2\;\ {\rm{TeV}} $, supporting the narrow-width factorization of $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $. For $ m_a\gtrsim 2.5\;\ {\rm{TeV}} $, the required couplings lead to a rapidly growing width ($ \Gamma_a/m_a\gtrsim 0.1 $ and up to $ {\cal{O}}(1) $ in the extreme corner), where finite-width and/or UV-completion effects are expected to be important. In channels where the resonance mass cannot be reconstructed (as in $ VV\to 2\ell2\nu $), off-shell contributions and signal–background interference can become relevant once the width is not small [100]; hence, the broad-width corner should be interpreted as an extrapolation. Accordingly, our main results are also presented as model-independent thresholds on $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $, which can be mapped onto specific UV completions once the relevant scales and Wilson coefficients are specified.

Partial-wave unitarity (conservative guide). Following Ref. [88], which derives partial-wave unitarity bounds in the gauge-invariant ALP EFT up to dimension-6, we use its dimension-5 truncation as a baseline guide. The strongest perturbative constraint arises from $ VV\to Va $ with one external ALP (dominant $ J=1 $ eigenvalue), primarily bounding the SU(2)$ _L $ coefficient $ c_{\widetilde{W}} $. Imposing $ |T_J|\le 1 $ in the high-energy limit yields

using $ g_{aWW}=4c_{\widetilde{W}}/f_a $. Because$ \sqrt{\hat{s}} $ is not directly observable event-by-event, we treat Eq. (A2) as a conservative guide for where the minimal dimension-5 truncation is expected to require UV unitarisation, rather than as a strict cut on our event selection.

Unitarity vs. EFT cutoff and UV sensitivity. The cutoff Λ controls the breakdown of the derivative expansion and is UV dependent. Parametrically one expects $ \Lambda \lesssim {4}/{|g_{aWW}|} $ for weakly-coupled UV completions, while NDA for strong dynamics suggests $ \Lambda \lesssim {16\pi}/{|g_{aWW}|} $ (cf. Eq. (11)). By contrast, Eq. (A2) constrains the perturbative description of specific high-energy $ 2\to2 $ channels in the truncated EFT and can be saturated already below the NDA cutoff estimate for the minimal dimension-5 truncation [88]. In that regime, additional UV ingredients (new states, higher-dimensional operators, or form factors) are expected to unitarise the amplitude, and the numerical bound can be relaxed once such effects are included (see, e.g., the modified limits when allowing dimension-6 contributions and cancellations in Ref. [88]). For collider applications, it is thus useful to treat the unitarity curve as a conservative validity guide and, if desired, to assess robustness by restricting to events below a chosen hard scale (e.g., an upper cut on a diboson invariant mass as in ALP-VBS studies [80]). Finally, because the high-energy rate in a UV-complete (unitarised) realization need not map one-to-one onto a single EFT coupling, we also provide the sensitivity in terms of $ \sigma(pp\to jja)\times {\rm{Br}}(a\to WW) $ for reinterpretation.

Total width estimate. Although $ m_a $ is not directly reconstructible due to the two neutrinos, the on-shell (factorised) signal modeling motivates verifying whether $ \Gamma_a/m_a $ remains small in the quoted sensitivity region. For dimension-5 gauge couplings, one expects $ \Gamma_a\propto g_{aWW}^2\,m_a^3 $ at large $ m_a $ (and hence $ \Gamma_a/m_a\propto g_{aWW}^2\,m_a^2 $). We compute the total width $ \Gamma_a $ in the strict photophobic benchmark with $\tt{MadGraph5\_aMC@NLO} $ and report $ \Gamma_a/m_a $ along the $ 2\sigma $ and $ 5\sigma $ reach curves in Table A1. We find $ \Gamma_a/m_a\lesssim$ 5% (8%) on the $ 2\sigma $ ($ 5\sigma $) curve for $ m_a\lesssim 2\;\ {\rm{TeV}} $, supporting the narrow-width factorization of $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $. For $ m_a\gtrsim 2.5\;\ {\rm{TeV}} $, the required couplings lead to a rapidly growing width ($ \Gamma_a/m_a\gtrsim 0.1 $ and up to $ {\cal{O}}(1) $ in the extreme corner), where finite-width and/or UV-completion effects are expected to be important. In channels where the resonance mass cannot be reconstructed (as in $ VV\to 2\ell2\nu $), off-shell contributions and signal–background interference can become relevant once the width is not small [100]; hence, the broad-width corner should be interpreted as an extrapolation. Accordingly, our main results are also presented as model-independent thresholds on $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $, which can be mapped onto specific UV completions once the relevant scales and Wilson coefficients are specified.

Partial-wave unitarity (conservative guide). Following Ref. [88], which derives partial-wave unitarity bounds in the gauge-invariant ALP EFT up to dimension-6, we use its dimension-5 truncation as a baseline guide. The strongest perturbative constraint arises from $ VV\to Va $ with one external ALP (dominant $ J=1 $ eigenvalue), primarily bounding the SU(2)$ _L $ coefficient $ c_{\widetilde{W}} $. Imposing $ |T_J|\le 1 $ in the high-energy limit yields

using $ g_{aWW}=4c_{\widetilde{W}}/f_a $. Because$ \sqrt{\hat{s}} $ is not directly observable event-by-event, we treat Eq. (A2) as a conservative guide for where the minimal dimension-5 truncation is expected to require UV unitarisation, rather than as a strict cut on our event selection.

Unitarity vs. EFT cutoff and UV sensitivity. The cutoff Λ controls the breakdown of the derivative expansion and is UV dependent. Parametrically one expects $ \Lambda \lesssim {4}/{|g_{aWW}|} $ for weakly-coupled UV completions, while NDA for strong dynamics suggests $ \Lambda \lesssim {16\pi}/{|g_{aWW}|} $ (cf. Eq. (11)). By contrast, Eq. (A2) constrains the perturbative description of specific high-energy $ 2\to2 $ channels in the truncated EFT and can be saturated already below the NDA cutoff estimate for the minimal dimension-5 truncation [88]. In that regime, additional UV ingredients (new states, higher-dimensional operators, or form factors) are expected to unitarise the amplitude, and the numerical bound can be relaxed once such effects are included (see, e.g., the modified limits when allowing dimension-6 contributions and cancellations in Ref. [88]). For collider applications, it is thus useful to treat the unitarity curve as a conservative validity guide and, if desired, to assess robustness by restricting to events below a chosen hard scale (e.g., an upper cut on a diboson invariant mass as in ALP-VBS studies [80]). Finally, because the high-energy rate in a UV-complete (unitarised) realization need not map one-to-one onto a single EFT coupling, we also provide the sensitivity in terms of $ \sigma(pp\to jja)\times {\rm{Br}}(a\to WW) $ for reinterpretation.

Total width estimate. Although $ m_a $ is not directly reconstructible due to the two neutrinos, the on-shell (factorised) signal modeling motivates verifying whether $ \Gamma_a/m_a $ remains small in the quoted sensitivity region. For dimension-5 gauge couplings, one expects $ \Gamma_a\propto g_{aWW}^2\,m_a^3 $ at large $ m_a $ (and hence $ \Gamma_a/m_a\propto g_{aWW}^2\,m_a^2 $). We compute the total width $ \Gamma_a $ in the strict photophobic benchmark with $\tt{MadGraph5\_aMC@NLO} $ and report $ \Gamma_a/m_a $ along the $ 2\sigma $ and $ 5\sigma $ reach curves in Table A1. We find $ \Gamma_a/m_a\lesssim$ 5% (8%) on the $ 2\sigma $ ($ 5\sigma $) curve for $ m_a\lesssim 2\;\ {\rm{TeV}} $, supporting the narrow-width factorization of $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $. For $ m_a\gtrsim 2.5\;\ {\rm{TeV}} $, the required couplings lead to a rapidly growing width ($ \Gamma_a/m_a\gtrsim 0.1 $ and up to $ {\cal{O}}(1) $ in the extreme corner), where finite-width and/or UV-completion effects are expected to be important. In channels where the resonance mass cannot be reconstructed (as in $ VV\to 2\ell2\nu $), off-shell contributions and signal–background interference can become relevant once the width is not small [100]; hence, the broad-width corner should be interpreted as an extrapolation. Accordingly, our main results are also presented as model-independent thresholds on $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $, which can be mapped onto specific UV completions once the relevant scales and Wilson coefficients are specified.

We estimate the impact of a residual background-normalization uncertainty by treating the post-selection background yield $ N_b $ as a Gaussian-constrained nuisance parameter with width $ \sigma_b \equiv \delta_b\,N_b $. For a single-bin counting experiment, the median (Asimov) discovery significance from the profile-likelihood-ratio test statistic is [98, 99]

which reduces for $ N_s\ll N_b $ and $ \sigma_b^2\ll N_b $ to the familiar approximation $ \sigma_{\rm stat+syst}\simeq N_s/\sqrt{N_b+\sigma_b^2} $ [99].

As an illustration, we repeat the analysis at the benchmark point $ m_a=750\;\ {\rm{GeV}} $ and, for each assumed $ \delta_b $, choose the BDT score threshold by maximizing $ \sigma_{\rm stat+syst} $ in Eq. (A3) (rather than the statistical-only significance). As a result, the optimal post-selection background yield $ N_b $ depends on $ \delta_b $. Table A2 reports the corresponding $ 2\sigma $ and $ 5\sigma $ sensitivities to $ g_{aWW} $ and to the fiducial rate $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ for $ \delta_b $ = 0, 2, 5, 8, 10%. In realistic experimental analyses, dominant backgrounds are typically constrained with dedicated control regions, which can reduce the effective $ \delta_b $; therefore, we use statistical-only reaches as a transparent baseline in this work and quote this benchmark study only to indicate the possible degradation from percent-level systematics.

We estimate the impact of a residual background-normalization uncertainty by treating the post-selection background yield $ N_b $ as a Gaussian-constrained nuisance parameter with width $ \sigma_b \equiv \delta_b\,N_b $. For a single-bin counting experiment, the median (Asimov) discovery significance from the profile-likelihood-ratio test statistic is [98, 99]

which reduces for $ N_s\ll N_b $ and $ \sigma_b^2\ll N_b $ to the familiar approximation $ \sigma_{\rm stat+syst}\simeq N_s/\sqrt{N_b+\sigma_b^2} $ [99].

As an illustration, we repeat the analysis at the benchmark point $ m_a=750\;\ {\rm{GeV}} $ and, for each assumed $ \delta_b $, choose the BDT score threshold by maximizing $ \sigma_{\rm stat+syst} $ in Eq. (A3) (rather than the statistical-only significance). As a result, the optimal post-selection background yield $ N_b $ depends on $ \delta_b $. Table A2 reports the corresponding $ 2\sigma $ and $ 5\sigma $ sensitivities to $ g_{aWW} $ and to the fiducial rate $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ for $ \delta_b $ = 0, 2, 5, 8, 10%. In realistic experimental analyses, dominant backgrounds are typically constrained with dedicated control regions, which can reduce the effective $ \delta_b $; therefore, we use statistical-only reaches as a transparent baseline in this work and quote this benchmark study only to indicate the possible degradation from percent-level systematics.

We estimate the impact of a residual background-normalization uncertainty by treating the post-selection background yield $ N_b $ as a Gaussian-constrained nuisance parameter with width $ \sigma_b \equiv \delta_b\,N_b $. For a single-bin counting experiment, the median (Asimov) discovery significance from the profile-likelihood-ratio test statistic is [98, 99]

which reduces for $ N_s\ll N_b $ and $ \sigma_b^2\ll N_b $ to the familiar approximation $ \sigma_{\rm stat+syst}\simeq N_s/\sqrt{N_b+\sigma_b^2} $ [99].

As an illustration, we repeat the analysis at the benchmark point $ m_a=750\;\ {\rm{GeV}} $ and, for each assumed $ \delta_b $, choose the BDT score threshold by maximizing $ \sigma_{\rm stat+syst} $ in Eq. (A3) (rather than the statistical-only significance). As a result, the optimal post-selection background yield $ N_b $ depends on $ \delta_b $. Table A2 reports the corresponding $ 2\sigma $ and $ 5\sigma $ sensitivities to $ g_{aWW} $ and to the fiducial rate $ \sigma(pp\to jj\,a)\times{\rm{Br}}(a\to W^+W^-) $ for $ \delta_b $ = 0, 2, 5, 8, 10%. In realistic experimental analyses, dominant backgrounds are typically constrained with dedicated control regions, which can reduce the effective $ \delta_b $; therefore, we use statistical-only reaches as a transparent baseline in this work and quote this benchmark study only to indicate the possible degradation from percent-level systematics.

We thank Bin Diao, Zilong Ding, Yiheng Xiong and Ye Lu for helpful discussions.

Note added. Following standard practice in high energy physics, authors are listed in strict alphabetical order by surname. This ordering should not be interpreted as indicating any ranking of contribution, seniority, or leadership.

We thank Bin Diao, Zilong Ding, Yiheng Xiong and Ye Lu for helpful discussions.

Note added. Following standard practice in high energy physics, authors are listed in strict alphabetical order by surname. This ordering should not be interpreted as indicating any ranking of contribution, seniority, or leadership.

We thank Bin Diao, Zilong Ding, Yiheng Xiong and Ye Lu for helpful discussions.

Note added. Following standard practice in high energy physics, authors are listed in strict alphabetical order by surname. This ordering should not be interpreted as indicating any ranking of contribution, seniority, or leadership.