# Extreme Event Prediction with Multi-agent Reinforcement Learning-based
Parametrization of Atmospheric and Oceanic Turbulence

Rambod Mojgani
Department of Mechanical Engineering
Rice University
<EMAIL_ADDRESS>
&Daniel Waelchli
Chair of Computational Science
ETH Zurich
<EMAIL_ADDRESS>
&Yifei Guan
Department of Mechanical Engineering
Rice University
<EMAIL_ADDRESS>
&Petros Koumoutsakos
School of Engineering and Applied Sciences
Harvard University
<EMAIL_ADDRESS>
&Pedram Hassanzadeh
Department of Mechanical Engineering and
Department of Earth, Environmental and Planetary Sciences
Rice University
<EMAIL_ADDRESS>

###### Abstract

Global climate models (GCMs) are the main tools for understanding and
predicting climate change. However, due to limited numerical resolutions,
these models suffer from major structural uncertainties; e.g., they cannot
resolve critical processes such as small-scale eddies in atmospheric and
oceanic turbulence. Thus, such small-scale processes have to be represented as
a function of the resolved scales via closures (parametrization). The accuracy
of these closures is particularly important for capturing climate extremes.
Traditionally, such closures are based on heuristics and simplifying
assumptions about the unresolved physics. Recently, supervised-learned
closures, trained offline on high-fidelity data, have been shown to outperform
the classical physics-based closures. However, this approach requires a
significant amount of high-fidelity training data and can also lead to
instabilities. Reinforcement learning is emerging as a potent alternative for
developing such closures as it requires only low-order statistics and leads to
stable closures. In Scientific Multi-Agent Reinforcement Learning (SMARL)
computational elements serve a dual role of discretization points and learning
agents. Here, we leverage SMARL and fundamentals of turbulence physics to
learn closures for canonical prototypes of atmospheric and oceanic turbulence.
The policy is trained using only the enstrophy spectrum, which is nearly
invariant and can be estimated from a few high-fidelity samples (these few
samples are far from enough for supervised/offline learning). We show that
these closures lead to stable low-resolution simulations that, at a fraction
of the cost, can reproduce the high-fidelity simulations’ statistics,
including the tails of the probability density functions (PDFs). These results
demonstrate the high potential of SMARL for closure modeling for GCMs,
especially in the regime of scarce data and indirect observations.

##  1 Introduction

Predicting extreme weather and climate change effects demands simulations that
account for complex interactions across nonlinear processes occurring over a
wide range of spatiotemporal scales. Turbulence as manifested in atmospheric
and oceanic flows, is a prominent example of such nonlinear and multi-scale
processes, and plays a critical role in transporting and mixing momentum and
heat in the climate system. While the governing equations of turbulent flows
are known, GCMs, which are the main tools for predicting climate variability,
cannot resolve all the relevant scales. For example, in the atmosphere alone,
these scales span from $10^{-4}$ (and smaller) to $10^{4}$ km [1, 2]. Despite
advances in our computing capabilities for climate modeling, this limitation
is expected to persist for decades.

The effect of unresolved small scales, often referred to as sub-grid scales
(SGSs), cannot be ignored in a nonlinear system, and their two-way
interactions with the resolved scales have to be accurately accounted for in
order for the GCMs to produce stable simulations and the right statistics
(climate) and extreme events. Current GCMs use semi-empirical and physics-
based representations of SGSs using closures [3]. The input to a closure
function is the resolved scales and the output is the SGSs fluxes of momentum,
heat, etc. The current closures for many Earth system processes, particularly
turbulent flows, fall short of accurately representing the two-way
interactions, due to oversimplifications and incomplete theoretical
understanding [4, 5]. For example, a major shortcoming is that the current
closures are too diffusive (dissipative) and also do not represent a real and
important phenomenon called backscattering (basically, anti-diffusion) [6],
preventing the GCMs from capturing the extreme events [7, 8]. Recently, there
has been growing interest in using machine learning (ML), particularly deep
neural networks (DNNs), to learn closures from data. There are two general
approaches to ML-based data-driven closure modeling:

Figure 1: Schematic of training an SGS closure using SMARL. The invariants of
gradient and Hessian of velocity ($\lambda^{\nabla\overline{u}}$ and
$\lambda^{\nabla\nabla\overline{u}}$) are the states, the actions are the
localized (and dynamic) coefficient $c_{s}(x,y,t)$ of classical Smagorinsky or
$c_{l}(x,y,t)$ of Leith closures of the SGS term $\Pi(x,y,t)$, and the policy,
$\pi$ , learns matching the enstrophy spectrum $\hat{\epsilon}$ to that of the
DNS.

Supervised (offline) learning of closures: In this approach, many snapshots of
high-fidelity data (e.g., from direct numerical simulation (DNS)) are
collected as the “truth”, then filtered and coarse-grained to extract the SGS
terms [9]. In turn, these data are used to train a DNN to match the SGS terms
$\Pi$ from the closures and from the truth (e.g., using a mean-square-error
loss [4]). Once trained, the DNN is coupled to the low-resolution solver to
perform large eddy simulation (LES). Studies using a variety of architectures
and test cases, from canonical turbulent systems to atmospheric and oceanic
flows, have shown the possibility of outperforming classical physics-based
closures such as Smagorinsky and Leith [e.g., 4, 8, 10, 11, 12, 13, 14].
However, for many critical climate processes, such high-quality datasets are
scarce, extracting the SGS terms is not straightforward, and offline-learned
closures can lead to unstable runs [15, 16]. While adding physics to the DNN,
transfer learning, and other techniques can address these issues to some
degree [14, 17, 18, 19, 20], overall, offline learning remains a promising but
challenging approach to data-driven SGS modeling for climate applications.

Online learning of closures: Online learning is emerging as a potent
alternative to supervised learning with the closures learned while they
operate on the LES. The goal is not to match detailed flow quantities (such as
the velocity field or the SGS terms) but instead match the low-order
statistics of the high-fidelity simulations or observations. In the context of
climate modeling, depending on the application, these statistics could be key
properties of climate variability [20] or spectra of turbulent flows [21, 3].
Online learning requires running the numerical model (e.g., a GCM) during the
DNN training, which can be challenging. In general, 3 approaches to online
learning exists: 1) using a differentiable LES solver/GCM [22, 23, 24], 2)
using Ensemble Kalman Inversion (EKI) [25, 26, 27, 20], and 3) using
reinforcement learning [28, 29, 30, 31]. While 1-2 have shown promising
results, they face major challenges. For example, current GCMs are not
differentiable and this approach requires major development in the climate
modeling infrastructure [16].

Multi-agent reinforcement learning (MARL), however, has accomplished
previously unattainable solutions in ML tasks [e.g., 32, 33, 34], as well as
success in improving the parametric and structural uncertainties of closures
for 3D homogeneous and wall-bounded turbulent flows [28, 29, 31, 30]. However,
its potential in climate-relevant applications, particularly in capturing
extreme events, has remained unexplored. Here,

  * •

We train a SMARL-based SGS closure using low-order statistics in climate-
relevant flows, i.e., 2D quasi-geostrophic turbulence with different forcing
or $\beta$ effects, producing multi-scale jets and vortices like those
observed in the Earth’s atmosphere and ocean,

  * •

We remark that we use as input states to the SMARL invariants of the flow and
learn the flow-dependent coefficient $c$ of two classic physics-based closures
(Leith and Smagorinsky) by matching the enstrophy spectrum the LES solver with
that of the DNS (obtained from only 10 true snapshots).

  * •

To test the performance of the data-driven closure, we compare the kinetic
energy spectrum and vorticity PDF of the DNS with $160$ to $10240\times$
spatio-temporally coarser LES-MARL. We particularly focus on comparing the
tails of these PDFs, as they represent extreme (weather) events in these
prototypes. As baselines, we use localized dynamics Smagorinsky and Leith
closures, which approximate $c$ as a function of the flow based on physical
arguments.

##  2 Scientific Multi-Agent Reinforcement Learning (SMARL)

In deep-SMARL, a DNN is trained to learn a policy that maps the states to
actions. States are fed into a DNN and actions of the agents maximize a
reward, see the schematic in Figure 1. Below, we describe the main elements of
the training, for which we have utilized Korali, a general-purpose, high-
performance framework for deep-SMARL [35].

### State:

The state vector consists of a combination of local and global variables. As
local states, instantaneous invariants $\lambda$ of filtered velocity
gradients [36] and velocity Hessians [28] (5 non-zero local variables) are
used. This choice embeds Galilean invariance into the closure. As global
states, enstrophy spectrum $\hat{\epsilon}$ is used. We have found the use of
these physically motivated invariants, rather than $(\bar{u},\bar{v})$ or
$\bar{\psi}$ or their derivatives, to be key in learning successful closures.

### Action:

The SGS values at each grid point are required to evolve the governing
equations of the environment in time. To retain some degree of
interpretability and reduce the computational complexity of training, two
classical physics-based closures are employed as the main structure of the SGS
closure: (i) Smagorinsky [37], which uses
$\nu_{e}=c_{s}\Delta^{2}|\bar{\mathcal{S}}|$, and (ii) Leith [38], which uses
$\nu_{e}=c_{l}\Delta^{3}{|\nabla\bar{\omega}|}$, where $\nu_{e}(x,y,t)$ is the
eddy viscosity, $\Delta$ is the filter size, $c_{s}(x,y,t)$ and $c_{l}(x,y,t)$
are the key coefficients, and $|\bar{\mathcal{S}}|$, and
$|\nabla\bar{\omega}|$ are the magnitudes of the filtered strain rate and
$\bar{\omega}$ gradient tensors. The coefficients $c$, which cannot be
obtained from first principles, are considered as actions (learned as a
function of the state). The actions are interpolated between the agents on to
the LES grid via a bilinear scheme, and used to calculate the SGS stress
tensor, ${\tau}^{\text{SGS}}=-2\nu_{e}\bar{\mathcal{S}}$, which is then used
to compute $\Pi$ that is needed in the low-resolution LES solver, Eq. (1).

### Reward:

The goal of our SGS closure is to match a target enstrophy spectrum, which can
be calculated from a short-time high-fidelity simulation or few observations.
We have computed the spectrum using 10 snapshots from a short DNS run, which
is known to be insufficient to learn a successful closer offline in these
flows [17]. The reward $r_{t}$ at each time step is defined as the cumulative
sum $r_{t}=\Sigma_{s=0}^{t}r^{\prime}_{s}$ of the inverse of the
$\mathcal{L}_{2}$ errors of the logarithms of the spectra, i.e.,
$r^{\prime}_{s}=1/\|\log(\hat{\varepsilon}_{\text{DNS}})-\log(\hat{\varepsilon}_{s,\text{RL}})\|^{2}_{2}$,
where $\hat{\varepsilon}_{\text{DNS}}$ is the time-averaged enstrophy spectrum
and $\hat{\varepsilon}_{s,\text{RL}}$ is the instantaneous enstrophy spectrum
at step $s$. Note that both local variables of the state and the actions are
defined at the location of the agents, and agents are uniformly distribution
in the domain.

### Environment:

The vorticity–streamfunction formulation of the 2D Navier-Stokes equation
(NSE) is solved using a pseudo-spectral method. In all cases,
$\text{Re}=20000$ and the DNS resolution is $1024$ collocation points in each
direction. The solver is coupled with Korali as the environment. Briefly, the
environment provides the dynamics of LES given the actions and the states,

$\displaystyle\frac{\partial\overline{\omega}}{\partial
t}=\mathscr{F}(\overline{\omega},\overline{\psi})+\Pi,$ (1)

where
$\Pi(\overline{\psi},\overline{\omega})=\nabla\times(\nabla\cdot\tau^{\text{SGS}})$,
$\mathscr{F}(.)$ represents the linear and nonlinear terms of the NSE (see Eq.
(2)), and $\overline{\psi}$ and $\overline{\omega}$ are the resolved
streamfunction and vorticity on the coarse grid.

##  3 Experiments

We have developed closures for 3 different forced 2D turbulent flows on the
$\beta$-plane (details in Table 1). These cases are commonly used to evaluate
the SGS closures of geophysical turbulence [e.g., 39, 17, 13] and exhibit
distinct behaviors and dynamics, as seen in snapshots of vorticity fields,
$\omega$, in Figure 2. Training is performed with the objective of achieving
an LES enstrophy spectrum close to the target (true, DNS) spectrum. We have
used the kinetic energy spectra as one unseen test metric. More importantly,
the PDFs of the resolved vorticity are compared. The tails of these PDFs
represent rare, extreme events, i.e., significantly large $\omega$ with a
small probability of occurrence. Note that these vortices resemble the weather
system’s high- and low-pressure anomalies, which can cause various extreme
weather events [40]. LES with $16$ to $1024\times$ coarser spatial resolution
and $10\times$ larger time steps coupled to the learned closures are then ran,
and their statistics are compared with those of the DNS and LES with classical
dynamic Smagorinsky and Leith closures. As summarized in Figure 2, the tails
of the vorticity PDFs clearly show the advantage of the SMARL-based closures,
suggesting that these closures have the right amount of diffusion and
backscattering (anti-diffusion). The classical closures are too diffusive (a
known problem), leading to much less frequent extreme events. The energy
spectra also show the better ability of LES with SMARL-based closures in
capturing the energy across the scales.

(a)  Case 1: $\kappa_{f}=4$, $\beta=0$, $N=32$ ($10240\times$ spatio-
temporally coarser)

(b)  Case 2: $\kappa_{f}=4$, $\beta=20$, $N=32$ ($10240\times$ spatio-
temporally coarser)

(c)  Case 3: $\kappa_{f}=25$, $\beta=0$, $N=256$ ($160\times$ spatio-
temporally coarser)

Figure 2:  Comparison of LES statistics with SMARL-based closures and the
classical closures. Rows (a)-(c) correspond to the cases in Table 1.

##  4 Conclusion and future work

We have trained a DNN-based SMARL to develop closures for climate-relevant
turbulent flows. We show that these closures enable LES with much fewer
degrees of freedom than DNS to produce statistics, including energy spectra,
PDF, and most importantly, tails of the PDFs, that closely match those of the
DNS. Particularly, in terms of capturing extreme events, LES with MARL-based
closures significantly outperform LES with classical physics-based closures.
The classical closures and even many offline-learned DNNs produce unstable LES
unless they are made overly diffusive (e.g., by eliminating backscattering),
which comes at the cost of under-representing extreme events [10, 8]. With a
small number of samples from DNS, which were not enough to even train an DNN
offline [17], SMARL develops closures capable of capturing the statistics of
such extreme events, suggesting that both diffusion and backscattering, i.e.,
the two-way interaction of the resolved scales and SGS, are accurately
represented (analysis is in progress to quantify the interscale
energy/enstrophy transfers).

Immediate next steps include further analysis and interpretability of the $c$
distributions learned using SMARL and examining the out-of-distribution
generalizability of these closures. While offline-learned closures are not
expected to extrapolate (e.g., to a different climate) unless methods like
transfer learning are used [41], the online-learned closures can be made
generalizable by proper scaling of the invariants and spectra [42], which
might be possible with enough theoretical understanding of the changing
system, e.g., the warming climate [43]. Future work will focus on applying
this framework to intermediate-complexity and comprehensive climate models to
learn SGS closures (specially for offline-online learning [20]) and
systematically calibrate GCMs [44].

## Acknowledgments and Disclosure of Funding

PH acknowledges support from an ONR Young Investigator Award (No.
N00014-20-1-2722), a grant from the NSF CSSI Program (no. OAC-2005123), and by
the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt
Futures program. PK gratefully acknowledges support from the Air Force Office
of Scientific Research (MURI grant no. FA9550-21-1-005). Computational
resources were provided by NSF XSEDE/ACCESS (Allocations ATM170020 and
PHY220125).

## References

  * Klein [2010] Rupert Klein.  Scale-dependent models for atmospheric flows.  _Annual Review of Fluid Mechanics_ , 42(1):249–274, 2010.
  * Vallis [2017] Geoffrey K Vallis.  _Atmospheric and Oceanic Fluid Dynamics_.  Cambridge University Press, 2017.
  * Schneider et al. [2017] Tapio Schneider, Shiwei Lan, Andrew Stuart, and Joao Teixeira.  Earth system modeling 2.0: A blueprint for models that learn from observations and targeted high-resolution simulations.  _Geophysical Research Letters_ , 44(24):12–396, 2017.
  * Rasp et al. [2018] Stephan Rasp, Michael S Pritchard, and Pierre Gentine.  Deep learning to represent subgrid processes in climate models.  _Proceedings of the National Academy of Sciences_ , 115(39):9684–9689, 2018.
  * Hewitt et al. [2020] Helene T Hewitt, Malcolm Roberts, Pierre Mathiot, Arne Biastoch, Ed Blockley, Eric P Chassignet, Baylor Fox-Kemper, Pat Hyder, David P Marshall, Ekaterina Popova, et al.  Resolving and parameterising the ocean mesoscale in earth system models.  _Current Climate Change Reports_ , pages 1–16, 2020.
  * Jansen et al. [2015] Malte F. Jansen, Isaac M. Held, Alistair Adcroft, and Robert Hallberg.  Energy budget-based backscatter in an eddy permitting primitive equation model.  _Ocean Modelling_ , 94:15–26, 2015.  ISSN 1463-5003.
  * O’Gorman and Dwyer [2018] Paul A O’Gorman and John G Dwyer.  Using machine learning to parameterize moist convection: Potential for modeling of climate, climate change, and extreme events.  _Journal of Advances in Modeling Earth Systems_ , 10(10):2548–2563, 2018.
  * Guan et al. [2022] Yifei Guan, Ashesh Chattopadhyay, Adam Subel, and Pedram Hassanzadeh.  Stable a posteriori LES of 2D turbulence using convolutional neural networks: Backscattering analysis and generalization to higher Re via transfer learning.  _Journal of Computational Physics_ , 458:111090, 2022.
  * Grooms et al. [2021] Ian Grooms, Nora Loose, Ryan Abernathey, JM Steinberg, Scott Daniel Bachman, Gustavo Marques, Arthur Paul Guillaumin, and Elizabeth Yankovsky.  Diffusion-based smoothers for spatial filtering of gridded geophysical data.  _Journal of Advances in Modeling Earth Systems_ , 13(9):e2021MS002552, 2021.
  * Maulik et al. [2019] Romit Maulik, Omer San, Adil Rasheed, and Prakash Vedula.  Subgrid modelling for two-dimensional turbulence using neural networks.  _Journal of Fluid Mechanics_ , 858:122–144, 2019.
  * Bolton and Zanna [2019] Thomas Bolton and Laure Zanna.  Applications of deep learning to ocean data inference and subgrid parameterization.  _Journal of Advances in Modeling Earth Systems_ , 11(1):376–399, 2019.
  * Yuval and O’Gorman [2020] Janni Yuval and Paul A O’Gorman.  Stable machine-learning parameterization of subgrid processes for climate modeling at a range of resolutions.  _Nature communications_ , 11(1):3295, 2020.
  * Srinivasan et al. [2023] Kaushik Srinivasan, Mickael D. Chekroun, and James C. McWilliams.  Turbulence closure with small, local neural networks: Forced two-dimensional and $\beta$-plane flows.  _preprint arXiv:2304.05029_ , 2023.
  * Subel et al. [2022] Adam Subel, Yifei Guan, Ashesh Chattopadhyay, and Pedram Hassanzadeh.  Explaining the physics of transfer learning a data-driven subgrid-scale closure to a different turbulent flow.  _arXiv preprint arXiv:2206.03198_ , 2022.
  * Sun et al. [2023] Y Qiang Sun, Pedram Hassanzadeh, M Joan Alexander, and Christopher G Kruse.  Quantifying 3d gravity wave drag in a library of tropical convection-permitting simulations for data-driven parameterizations.  _Journal of Advances in Modeling Earth Systems_ , 15(5):e2022MS003585, 2023.
  * Schneider et al. [2023] Tapio Schneider, Swadhin Behera, Giulio Boccaletti, Clara Deser, Kerry Emanuel, Raffaele Ferrari, L Ruby Leung, Ning Lin, Thomas Müller, Antonio Navarra, et al.  Harnessing AI and computing to advance climate modelling and prediction.  _Nature Climate Change_ , 13(9):887–889, 2023\.
  * Guan et al. [2023] Yifei Guan, Adam Subel, Ashesh Chattopadhyay, and Pedram Hassanzadeh.  Learning physics-constrained subgrid-scale closures in the small-data regime for stable and accurate LES.  _Physica D: Nonlinear Phenomena_ , 443:133568, 2023.  ISSN 0167-2789.
  * Beucler et al. [2021a] Tom Beucler, Michael Pritchard, Stephan Rasp, Jordan Ott, Pierre Baldi, and Pierre Gentine.  Enforcing analytic constraints in neural networks emulating physical systems.  _Physical Review Letters_ , 126(9):098302, 2021a.
  * Pedersen et al. [2023] Christian Pedersen, Laure Zanna, Joan Bruna, and Pavel Perezhogin.  Reliable coarse-grained turbulent simulations through combined offline learning and neural emulation.  _preprint arXiv:2307.13144_ , 2023.
  * Pahlavan et al. [2023] Hamid A Pahlavan, Pedram Hassanzadeh, and M Joan Alexander.  Explainable offline-online training of neural networks for parameterizations: A 1D gravity wave-QBO testbed in the small-data regime.  _preprint arXiv:2309.09024_ , 2023.
  * Schneider et al. [2021] Tapio Schneider, Andrew M Stuart, and Jin-Long Wu.  Learning stochastic closures using ensemble Kalman inversion.  _Transactions of Mathematics and Its Applications_ , 5(1):tnab003, 2021.
  * Frezat et al. [2022] Hugo Frezat, JL Sommer, Ronan Fablet, Guillaume Balarac, and Redouane Lguensat.  A posteriori learning for quasi-geostrophic turbulence parametrization.  _April_ , 2(3):6, 2022.
  * Shankar et al. [2023] Varun Shankar, Romit Maulik, and Venkatasubramanian Viswanathan.  Differentiable turbulence.  _preprint arXiv:2307.03683_ , 2023.
  * Gelbrecht et al. [2022] M. Gelbrecht, A. White, S. Bathiany, and N. Boers.  Differentiable programming for earth system modeling.  _EGUsphere_ , 2022:1–17, 2022.
  * Kovachki and Stuart [2019] Nikola B Kovachki and Andrew M Stuart.  Ensemble kalman inversion: a derivative-free technique for machine learning tasks.  _Inverse Problems_ , 35(9):095005, 2019.
  * Dunbar et al. [2021] Oliver RA Dunbar, Alfredo Garbuno-Inigo, Tapio Schneider, and Andrew M Stuart.  Calibration and uncertainty quantification of convective parameters in an idealized GCM.  _Journal of Advances in Modeling Earth Systems_ , 13(9):e2020MS002454, 2021.
  * Schneider et al. [2020] Tapio Schneider, Andrew M Stuart, and Jin-Long Wu.  Imposing sparsity within ensemble Kalman inversion.  _arXiv preprint arXiv:2007.06175_ , 2020.
  * Novati et al. [2021a] Guido Novati, Hugues Lascombes de Laroussilhe, and Petros Koumoutsakos.  Automating turbulence modelling by multi-agent reinforcement learning.  _Nature Machine Intelligence_ , 3(1):87–96, 2021a.
  * Kurz et al. [2023] Marius Kurz, Philipp Offenhäuser, and Andrea Beck.  Deep reinforcement learning for turbulence modeling in large eddy simulations.  _International Journal of Heat and Fluid Flow_ , 99:109094, 2023.  ISSN 0142-727X.
  * Bae and Koumoutsakos [2022] H Jane Bae and Petros Koumoutsakos.  Scientific multi-agent reinforcement learning for wall-models of turbulent flows.  _Nature Communications_ , 13(1):1443, 2022.
  * Kim et al. [2022] Junhyuk Kim, Hyojin Kim, Jiyeon Kim, and Changhoon Lee.  Deep reinforcement learning for large-eddy simulation modeling in wall-bounded turbulence.  _Physics of Fluids_ , 34(10), 2022.
  * Silver et al. [2017] David Silver, Thomas Hubert, Julian Schrittwieser, Ioannis Antonoglou, Matthew Lai, Arthur Guez, Marc Lanctot, Laurent Sifre, Dharshan Kumaran, Thore Graepel, Timothy P. Lillicrap, Karen Simonyan, and Demis Hassabis.  Mastering chess and Shogi by self-play with a general reinforcement learning algorithm.  _CoRR_ , abs/1712.01815, 2017.
  * Brown and Sandholm [2019] Noam Brown and Tuomas Sandholm.  Superhuman AI for multiplayer poker.  _Science_ , 365(6456):885–890, 2019.  ISSN 0036-8075.
  * OpenAI et al. [2019] OpenAI, :, Christopher Berner, Greg Brockman, Brooke Chan, Vicki Cheung, Przemysław Debiak, Christy Dennison, David Farhi, Quirin Fischer, Shariq Hashme, Chris Hesse, Rafal Józefowicz, Scott Gray, Catherine Olsson, Jakub Pachocki, Michael Petrov, Henrique P. d. O. Pinto, Jonathan Raiman, Tim Salimans, Jeremy Schlatter, Jonas Schneider, Szymon Sidor, Ilya Sutskever, Jie Tang, Filip Wolski, and Susan Zhang.  Dota 2 with Large Scale Deep Reinforcement Learning.  _arXiv e-prints_ , art. arXiv:1912.06680, December 2019.
  * Martin et al. [2022] Sergio M. Martin, Daniel Wälchli, Georgios Arampatzis, Athena E. Economides, Petr Karnakov, and Petros Koumoutsakos.  Korali: Efficient and scalable software framework for bayesian uncertainty quantification and stochastic optimization.  _Computer Methods in Applied Mechanics and Engineering_ , 389:114264, 2022.  ISSN 0045-7825.
  * Ling et al. [2016] Julia Ling, Andrew Kurzawski, and Jeremy Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  _Journal of Fluid Mechanics_ , 807:155–166, 2016.
  * Smagorinsky [1963] Joseph Smagorinsky.  General circulation experiments with the primitive equations: I. the basic experiment.  _Monthly Weather Review_ , 91(3):99 – 164, 1963\.
  * Leith [1996] C.E. Leith.  Stochastic models of chaotic systems.  _Physica D: Nonlinear Phenomena_ , 98(2):481–491, 1996.  ISSN 0167-2789.  Nonlinear Phenomena in Ocean Dynamics.
  * Frezat et al. [2021] Hugo Frezat, Guillaume Balarac, Julien Le Sommer, Ronan Fablet, and Redouane Lguensat.  Physical invariance in neural networks for subgrid-scale scalar flux modeling.  _Physical Review Fluids_ , 6(2):024607, 2021.
  * Woollings et al. [2018] Tim Woollings, David Barriopedro, John Methven, Seok-Woo Son, Olivia Martius, Ben Harvey, Jana Sillmann, Anthony R Lupo, and Sonia Seneviratne.  Blocking and its response to climate change.  _Current climate change reports_ , 4:287–300, 2018.
  * Subel et al. [2021] Adam Subel, Ashesh Chattopadhyay, Yifei Guan, and Pedram Hassanzadeh.  Data-driven subgrid-scale modeling of forced Burgers turbulence using deep learning with generalization to higher Reynolds numbers via transfer learning.  _Physics of Fluids_ , 33(3):031702, 2021.
  * Novati et al. [2021b] Guido Novati, Hugues Lascombes de Laroussilhe, and Petros Koumoutsakos.  Automating turbulence modelling by multi-agent reinforcement learning.  _Nature Machine Intelligence_ , pages 1–10, 2021b.
  * Beucler et al. [2021b] Tom Beucler, Michael Pritchard, Janni Yuval, Ankitesh Gupta, Liran Peng, Stephan Rasp, Fiaz Ahmed, Paul A O’Gorman, J David Neelin, Nicholas J Lutsko, et al.  Climate-invariant machine learning.  _arXiv preprint arXiv:2112.08440_ , 2021b.
  * Balaji et al. [2022] V Balaji, Fleur Couvreux, Julie Deshayes, Jacques Gautrais, Frédéric Hourdin, and Catherine Rio.  Are general circulation models obsolete?  _Proceedings of the National Academy of Sciences_ , 119(47):e2202075119, 2022.
  * Chandler and Kerswell [2013] Gary J Chandler and Rich R Kerswell.  Invariant recurrent solutions embedded in a turbulent two-dimensional Kolmogorov flow.  _Journal of Fluid Mechanics_ , 722:554–595, 2013.
  * Kochkov et al. [2021] Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer.  Machine learning–accelerated computational fluid dynamics.  _Proceedings of the National Academy of Sciences_ , 118(21), 2021.
  * Pope [2001] Stephen B Pope.  _Turbulent Flows_.  IOP Publishing, 2001.
  * Sagaut [2006] Pierre Sagaut.  _Large eddy simulation for incompressible flows: An introduction_.  Springer Science & Business Media, 2006.

##  Appendix A 2D turbulence

We consider the dimensionless governing equations in the vorticity ($\omega$)
and streamfunction ($\psi$) formulation in a doubly periodic square domain
with length $L=2\pi$, i.e.,

$\displaystyle\frac{\partial\omega}{\partial t}+\mathcal{N}(\omega,\psi)$
$\displaystyle=$ $\displaystyle\frac{1}{\text{Re}}\nabla^{2}\omega-
f-r\omega+\beta\frac{\partial\psi}{\partial x},$ (2a)
$\displaystyle\nabla^{2}\psi$ $\displaystyle=$ $\displaystyle-\omega,$ (2b)

where
$\displaystyle\mathcal{N}(\omega,\psi)=\left(\nicefrac{{\partial\psi}}{{\partial
y}}\right)\left(\nicefrac{{\partial\omega}}{{\partial
x}}\right)-\left(\nicefrac{{\partial\psi}}{{\partial
x}}\right)\left(\nicefrac{{\partial\omega}}{{\partial y}}\right),$ is the
nonlinear advection term, and
$f(x,y)=\kappa_{f}[\cos{(\kappa_{f}x)}+\cos{(\kappa_{f}y)}]$ is a
deterministic forcing [e.g., 45, 46].

To derive the equations for large eddy simulation (LES), we apply sharp
spectral filtering [47, 48], denoted by $\overline{(\cdot)}$, to Eq. (2) to
obtain

$\displaystyle\frac{\partial\overline{\omega}}{\partial
t}+\mathcal{N}(\overline{\omega},\overline{\psi})$ $\displaystyle=$
$\displaystyle\frac{1}{\text{Re}}\nabla^{2}\overline{\omega}-\overline{f}-r\overline{\omega}+\beta\frac{\partial\overline{\psi}}{\partial
x}+\underbrace{\mathcal{N}(\overline{\omega},\overline{\psi})-\overline{\mathcal{N}({\omega},{\psi})}}_{\Pi=\nabla\times(\nabla\cdot\tau^{\text{SGS}})},$
(3a) $\displaystyle\nabla^{2}\overline{\psi}$ $\displaystyle=$
$\displaystyle-\overline{\omega}.$ (3b)

The LES is solved on a coarse resolution with the sub-grid scale (SGS) term,
$\Pi$, being the unclosed term, requiring a model connecting it to the
resolved flow variables, i.e., closure.

For eddy viscosity models:

$\displaystyle{\tau}^{\text{SGS}}=-2\nu_{e}\bar{\mathcal{S}},$ (4)

where $\mathcal{S}_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial
x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)$ and
$\bar{\mathcal{S}}=\sqrt{2\mathcal{S}_{ij}\mathcal{S}_{ij}}$.

##  Appendix B Test Cases

The studied cases are summarized in Table 1.

Table 1: The test cases and hyper-parameters in training. $\sigma$ is standard deviation. Case | Re | $\beta$ | $\kappa_{f}$ | $r$ | $\sigma(\omega)$ | $\Delta t_{\text{RL}}/\Delta t_{\text{DNS}}$ | Training horizon | Updates policy every
---|---|---|---|---|---|---|---|---
1 | $20\times 10^{3}$ | $0$ | $4$ | $0.1$ | $5.51$ | $10$ | $1\times 10^{4}\Delta t_{\text{RL}}$ | $10\Delta t_{\text{RL}}$
2 | $20\times 10^{3}$ | $20$ | $4$ | $0.1$ | $10.75$ | $10$ | $2\times 10^{4}\Delta t_{\text{RL}}$ | $20\Delta t_{\text{RL}}$
3 | $20\times 10^{3}$ | $0$ | $25$ | $0.1$ | $13.01$ | $10$ | $1\times 10^{4}\Delta t_{\text{RL}}$ | $10\Delta t_{\text{RL}}$


∎

11institutetext: C. Ding 22institutetext: Badminton School 33institutetext: B.
Fuller 44institutetext: Sherborne School 55institutetext: E. Jones
66institutetext: University of Warwick 77institutetext: A. Martin
88institutetext: The Royal Latin School 99institutetext: W. Murray
1010institutetext: STFC Rutherford Appleton Laboratory & University of Warwick
1010email<EMAIL_ADDRESS>
# Electroweak-QCD interference in hadronic vector bosons decays at the LHC

Cathy Ding  Ben Fuller  Eleanor Jones  Aran Martin  William Murray

(Received: date / Accepted: date)

###### Abstract

The analysis of hadronic vector bosons decays at the LHC does not normally
consider interference with QCD $q\overline{q}$ production. This paper explores
the effect of this interference on the reconstructed peak positions and rates
for several production modes. In particular, boosted vectors bosons and vector
boson pairs are considered for the first time. Shifts of several GeV/c2 are
seen in the positions of the $W$ and $Z$ boson’s peaks, with a magnitude that
depends strongly on the kinematics involved. For boosted vector bosons in
regions currently studied, the effects are all found to be very small or
negligible. If experiments were to access lower transverse momenta, for
example in the experimental trigger systems, or in studies of low-recoil $VV$
decaying semileptonically, the effect of interference could be much larger.

##  1 Introduction

This paper presents a study of interference effects between electroweak
production of hadronic vector bosons and QCD background at the LHC. Such
effects are neglected in the simulation strategy generally adopted by the LHC
experimental collaborations, where vector bosons are simulated independently
from QCD. The goal of this paper is test the validity of that approach. Early
theoretical studies of interference ranft79  explored its impact on total $W$
and $Z$ (collectively referred to here as $V$) cross-sections in $pp$ and
$p\overline{p}$ collisions at energies below 1 TeV/c. More details were
explored theoretically at energies of 630 GeV/c and 1.8 TeV/c glover89  and
these studies concluded that “effects caused by the interference between
electroweak and QCD amplitudes must be taken into account when data are
compared with theoretical predictions”. A shift in the masses of the resonance
peaks of about 0.3 GeV/c2 was noted and it was also observed that this effect
was increased by experimental resolution. First evidence for hadronic vector
boson decay was shown in 1990 by UA2 in two-jet decays Alitti:210698 , using
4.66pb-1 of sp$\overline{\mathrm{p}}$s data. This result allowed for a 2.2 GeV
mass shift in the experimental peak position from interference, based on Ref.
glover89 . Pumplin PhysRevD.45.806  considered hadronic vector boson
observation at the Tevatron as a mass calibration channel and remarked that
interference causes a shift in mass of approximately 0.35 GeV/c2 downwards.
All of these studies were concerned with inclusive production at leading order
in QCD. No experimental evidence for this inclusive process has been published
by the LHC experiments.

There are several production modes of boosted hadronic vector bosons at the
LHC. Baur baur:2006sn  considered hadronic $V+X$ production, finding it to be
at the order of 1% of the jet rate. Interest in hadronic vector bosons at the
LHC was stimulated by the proposal to analyse boosted hadronic Higgs decays
Butterworth:2008iy . Inclusive $W\rightarrow q\bar{q}$ decays at high momentum
have been used as a calibration in hadronic $X\rightarrow
VV\mbox{$\rightarrow$}\mbox{$q\overline{q}$}\mbox{$q\overline{q}$}$ searches
Aad:2019fbh ; Sirunyan:2019jbg . Both CMS and ATLAS have also used boosted
$Z\rightarrow\mbox{$b\overline{b}$}$ decays as a method of validating
techniques to search for highly boosted Higgs bosons decaying to the same
final state PhysRevLett.120.071802 ; ATLAS-CONF-2018-052 . Recently, the ATLAS
collaboration performed a cross-section measurement in the $Z(\rightarrow
b\bar{b})\gamma$ channel using 36 fb-1 of data Aad:2019wdr  and this same
channel has also been used to validate $b$-tagging developments Aad:2680245 .
None of these studies considered interference.

Semileptonic $VV$ production provides experimentally accessible signals. For
example, the $VH,H\rightarrow\mbox{$b\overline{b}$}$ discovery papers by ATLAS
and CMS Aaboud:2018zhk ; Sirunyan:2018kst  obtained clear evidence for the
$VZ,Z\rightarrow\mbox{$b\overline{b}$}$ process, summing over
$V\mbox{$\rightarrow$}\nu\nu$, $l\nu$ and $l^{+}l^{-}$. These searches used
lower pT thresholds on the bosons than the inclusive ones discussed above,
which it will be shown gives rise to larger interference effects.

The multiple distinctive features of top quark pairs allow the extraction of a
clear hadronic $W$ boson signal, which has been used to validate jet
performance CMS-PAS-JME-16-003 ; ATLAS-CONF-2016-008  and calibrate the top
quark mass measurement Aaboud2019 ; Sirunyan2018 . A shift in the $W$ mass
peak here could have a significant impact. This is considered briefly in the
concluding section of this paper.

At this stage, it may be helpful to highlight an important point: if two
amplitudes $a$ and $b$ contribute to a process, the cross-section is given by
$\sigma=aa^{*}+bb^{*}+ab^{*}+ba^{*}$. If $a$ is the process under study and
$b$ represents a background, then the fractional contribution of the
interference term is at most $2|b|/|a|$. Thus, (fractional) interference
effects are expected to scale inversely with the square-root of the signal-to-
background ratio.

Simulations of several processes are studied at leading order using Sherpa
Gleisberg:2008ta  and the procedure is discussed in section 2. At higher
order, there are additional effects that are important in the prompt
production of heavy quark mesons at hadron colliders. This is especially
evident in $\psi(2s)$ production, where at high pT even Next-to-Leading-Order
QCD (NLO) calculations under-predict the rate by a factor of two
Cacciari:2012ny . The OpenLoops Cascioli:2011va  package was explored as an
NLO generator, but no relevant states are currently interfaced to Sherpa. The
results for the inclusive vector bosons, for comparison with previous
literature, are given in section 3, and, in addition, the analysis methodology
is introduced there. This is then applied to inclusive boosted $W$ and $Z$
bosons in section 4, $V+\gamma$ in section 5 and $VV$ in section 6. Section 7
contains a discussion of the results and prospects for experimental
verification.

##  2 Methodology

The results presented in this paper are obtained using the Sherpa 2.2.8
package with Photos Schonherr:2008av  and normally the NNPDF3.0 Ball:2014uwa
Parton Distribution Function (PDF) set, with MMHT2014nlo68cl Harland-Lang2015
used for certain studies. The Comix Gleisberg:2008fv  matrix element generator
is used throughout to generate leading order cross-sections. As such,
$s$-channel gluon-mediated processes are pure colour octets and do not
interfere with vector bosons. Furthermore, the parton shower, hadronization
and underlying event are mostly disabled to save computing time, so the
resulting output is a small set of partons only. For a small number of
examples, the parton shower is enabled and the DelphesMC deFavereau2014
simulation is employed, using the ‘ATLAS’ detector card to approximate the
influence of a detector on the results. These cases require a jet treatment,
discussed in section 3.1.1.

Sherpa uses the concepts of a ‘process’ and a ‘selector’. The process defines
the initial and final state particles, with additional options such as
defining intermediate states or the order of the electroweak coupling at the
amplitude level. The aim of this paper is to study the resonant vector boson
peaks, and their interference with the continuum background, but there is no
simple way to include all non-resonant electroweak processes in the
background. Instead, the three samples generated are a total, a background
(which is referred to as QCD throughout), with the minimum electroweak order
to reach the specified final state, and the electroweak signal, where two
additional orders in the electroweak coupling are required corresponding to
the creation and decay of an internal electroweak boson. So, for example,
$\gamma\mbox{$u\overline{d}$}$ is generated at order one and order three in
the electromagnetic coupling, though most states involve order zero and order
two. The selector imposes kinematic selections on the outgoing particles. In
this paper, the selectors applied are:

  * •

a mass window on the relevant quark-antiquark pair of 50 GeV$/c^{2}<m<130$
GeV$/c^{2}$. This wide window enables studies involving the experimental
resolution.

  * •

all final state particles must have a pT greater than 25 GeV/c. This mimics
experimental selection and stabilises the Sherpa integration.

  * •

all charged objects must have $|\eta|<2.5$. This mimics central detector
acceptance and also serves to stabilise the Sherpa integration in some cases.

  * •

when studying boosted bosons, the quark-antiquark pair is required to pass a
pT selection. In the case where there are three coloured objects in the final
state, this requirement is instead placed on the leading one, which, at this
matrix element level, has identical pT to the recoiling pair.

Simulation is done at 13 TeV centre-of-mass energy for proton-proton
collisions. A full ‘Run.dat’ file defining one process can be seen in appendix
A. Within this, it can be seen that the output is in hepMC format. DelphesMC
is then used to reformat this into a file that can be processed by root 6.16
Brun:1997pa .

The analysis extracts the interference component as the QCD background and
electroweak signal subtracted from the total and comparisons are made between
the pure electroweak signal and the electroweak plus the interference term
($EW+I$). The invariant mass of the $q\overline{q}$ pair is examined. For
processes where there are only two quarks in the final state, there is only
one combination to form the invariant mass. However, for processes with a
$q\overline{q}$ pair plus an extra quark in the final state, there are two
possible combinations and both are considered.

##  3 Vector boson production at rest

In this section, the simplest case of $W$ and $Z$ production with no recoil,
$p\overline{p}\rightarrow V$, is considered. This mode shows clear
interference effects and so is also used to check the impact of the selection,
the experimental resolution and additional gluon radiation.

At leading order, $Z$ or $W$ boson production can only proceed through quark-
antiquark annihilation. The experimental background has large contributions
from processes involving gluons in the initial or final state. However,
interference only occurs between processes with identical initial and final
states, so these background processes with gluons will not contribute. To
study interference, it is sufficient to focus on
$q\overline{q}$$\rightarrow$$q\overline{q}$. The total jet-jet cross-section
is two orders of magnitude larger than the $u\overline{u}$ cross-section so
this represents a considerable calculation simplification. Figure 1 shows the
Feynman diagrams for the electroweak and QCD processes.

Figure 1: The leading order $q\overline{q}$$\rightarrow$ $q\overline{q}$
diagrams, with $s$\- and $t$-channel $Z$ and gluon processes. The $Z$ could be
replaced by a $W$ or photon.

Only the $s$-channel electroweak process produces a resonant mass signature in
the $q\overline{q}$ final state and it is a colour-singlet. However, the
similar $s$-channel QCD diagram is a colour octet and, as a result, it does
not interfere. This does not hold in general for the $t$-channel gluon
exchange process. Therefore, interference requires identical initial and final
state particles. However, while the experimentalist has some possibility to
distinguish final state quarks, all $q\overline{q}$ initial states must be
included. This means, for example, that the $t$-channel QCD diagrams for a
$b\overline{b}$ final state are suppressed compared to the $u\overline{u}$
final state as a result of the relatively small $b$ quark PDF. Consequently,
it is to be expected that the apparent interference is reduced for heavy
quarks.

###  3.1 Inclusive $Z$ production

To study $Z$ production, the five final state quark-pairs are simulated
separately. The differential cross-section for
$q\overline{q}$$\rightarrow$$u\overline{u}$ is shown in figure 2, with the QCD
and electroweak channels compared with the combined process.

Figure 2: The process $q\overline{q}\rightarrow\mbox{$u\overline{u}$}$. Left
shows the pure electroweak (pale/green), pure QCD (dark/red) and total (black)
cross-sections. Right has the same electroweak component, but the hollow black
points show the total minus the QCD component, and the dark (blue) is the
interference term. The fit is discussed in the text.

There are sizeable interference effects in the mass spectrum, with an
enhancement below the peak and destructive interference above. These effects
are largest about $\pm\Gamma/2$ from the pole mass, but they also induce a
small shift in the peak position. They are modelled using the sum of a
relativistic Breit-Wigner function and an arbitrary constant amplitude:

$\mathcal{A}=\frac{2m_{0}\sqrt{\kappa(m/m_{0})^{3}}}{(m^{2}-m_{0}^{2})+im_{0}\Gamma}+\sqrt{b}+\sqrt{c}i$
(1)

In this formula, $\kappa$ is the normalisation of the Breit-Wigner, $m$ is the
invariant mass of the system, $m_{0}$ and $\Gamma$ are the particle mass and
width, respectively. $b$ and $c$ are the coefficients of the real and
imaginary amplitude components of an interfering background, respectively.
This formula neglects variation in the interfering amplitude with mass. More
importantly, in the identification of the $b$ and $c$ terms with the
interfering part of the QCD contribution, it ignores $t$-channel electroweak
processes. The motivation for the $m^{3}$ dependence of the Breit-Wigner rate
is heuristic.

The cross-section is not just $\mathcal{A}\mathcal{A}^{*}$, but has an
additional term, $d(m)$, which is a non-interfering background component
modelled with a second-order polynomial dependence on $m$:

$\displaystyle\frac{\partial\sigma}{\partial m}$ $\displaystyle=$
$\displaystyle\frac{2m_{0}\sqrt{\kappa(m/m_{0})^{3}}}{(m^{2}-m_{0}^{2})^{2}+m_{0}^{2}\Gamma^{2}}\left(2m_{0}\sqrt{\kappa(m/m_{0})^{3}}+\sqrt{b}(m^{2}-m_{0}^{2})+\sqrt{c}m_{0}\Gamma\right)~{}+b+c+d(m)~{}\mathrm{pb}~{}[\mbox{GeV/c${}^{2}$}]^{-1}.$
(2)

These nine parameters are simultaneously fit to the three distributions in
figure 2 (left), with terms dropped for the partial distributions where they
do not contribute. The extracted values are $m_{Z}=91.125\pm 0.001$ GeV/c2
(c/f 91.19 GeV/c2 used in simulation) and $\Gamma_{z}=2.687\pm 0.003$ GeV/c2
(c/f 2.50 GeV/c2 used in simulation). The increase in $\Gamma_{z}$ with
respect to the simulation is due to the fit attempting to describe the
$t$-channel exchange in the electroweak component.

The fit also finds $b=202\pm 2~{}\mathrm{pb}~{}[\mbox{GeV/c${}^{2}$}]^{-1}$
and $c=2.2\pm 0.2~{}\mathrm{pb}~{}[\mbox{GeV/c${}^{2}$}]^{-1}$; a sizeable
real term and a significant imaginary one, which is not expected and is
interpreted as being due to the $t$-channel electroweak exchange. The size of
the interfering QCD component $b$ is approximately one eighth of the total QCD
contribution to $q\overline{q}$$\rightarrow$$u\overline{u}$ at $m_{Z}$. The
components of the fit are shown in figure 2 (right) and it can be seen that
the approximate features are captured, but the limitations are clearly
visible. The fit is informative, but due to the missing electroweak
$t$-channel, it is not quantitatively reliable and simpler methods are
generally used for the rest of this paper.

A parameter, $\eta_{\mbox{scale}}$, characterising the scaling of the signal
strength is defined as the integral of the $EW+I$ component in a window of
$\pm$10 GeV/c2 around the nominal mass, divided by the same integral for the
electroweak term alone. This is the same as the definition of $\eta$ presented
by the ‘New Physics’ Working Group Brooijmans:2018xbu , but the integration
limits have been reduced from the $\pm 10\Gamma$ window to one of $\pm$10
GeV/c2 which is comparable to the typical experimental resolution. This
revised definition has less sensitivity to the $t$-channel component. The
other parameter extracted is the shift of the peak position. This is
established by fitting a Gaussian to the observed data, in a mass range
defined by the half-maxima of the distribution. The shifts in the peak
positions due to interference are evaluated for all the $Z$ decay modes and
given in table 1.

| $d\overline{d}$ | $u\overline{u}$ | $s\overline{s}$ | $c\overline{c}$ | $b\overline{b}$
---|---|---|---|---|---
signal-to-background | 0.60 | 0.43 | 1.08 | 1.18 | 2.6
Shift, GeV/c2 |  $-0.384\pm$0.014 |  $-0.405\pm$0.009 |  $-0.202\pm$0.010 |  $-0.135\pm$0.011 |  $-0.066\pm$0.007
$\eta_{\mbox{scale}}$ | 0.86$\pm$0.01 | 0.77$\pm$0.01 | 0.95$\pm$0.01 | 0.91$\pm$0.01 | 1.009$\pm$0.003
Table 1: Interference effects on the
$\mbox{$q\overline{q}$}\mbox{$\rightarrow$}W$ peak by final state quark
flavour. Signal-to-background is defined at the resonance peak without
experimental resolution and using only the initial and final states quoted,
with negligible statistical error.

The peak shift, averaged over the final state quark species, is -0.23$\pm$0.01
GeV/c2, similar to the $0.3-0.35$ GeV/c2 shifts reported in the previous work
discussed in the introduction, but smaller with the higher beam energy, as
expected. The magnitude of the effects falls as the quark mass rises, since
the PDF density reduces, and so the non-interfering $s$-channel QCD process
becomes more dominant.

These shifts are for the theoretical peak in the cross-section. Estimating the
expected change in the measured peak position requires allowing for detector
resolution. This was discussed by the ‘New Physics’ Working Group
PhysRevD.86.073016 , and Baur and Glover glover89 , where the latter adopted a
10% Gaussian smearing to approximate detector resolution effects. The same 10%
is used here, noting that ATLAS found a mass resolution between 8% and 16% for
boosted dijets with pT between 300 and 1000 GeV/c Aad:2019fbh  and CMS shows
between 7.5% and 10% for $W$ and $Z$ from the decay of 1200 - 5000 GeV/c2
bosons Sirunyan:2019jbg .

Figure 3: Left: The process $q\overline{q}\rightarrow u\overline{u}$ after
smearing by an assumed 10% resolution. The electroweak component is
pale/green, but the hollow black points show the $EW+I$. Right is the same
distribution for $b\overline{b}\rightarrow b\overline{b}$, showing a
significantly larger effect.

Figure 3 (left) shows that after smearing, the $EW+I$ component is negative
both above and below the peak, unlike what is seen in the original
distribution, as the smearing makes the $t$-channel contribution fractionally
larger. The $t$-channel imposes a decrease in differential cross-section
across the mass spectrum, and effectively a baseline shift under the peaks.
This gives rise to the reported reduction in $\eta_{\mbox{scale}}$, 0.69. If
the background had been extracted from a fit to sideband data, this reduction
in cross-section would not be seen. The fit to the electroweak alone yields
good agreement with the $Z$ mass, while the $EW+I$ distribution is down-
shifted by $3.35\pm 0.14$ GeV/c2, as a result of the interference.

This can be compared with the shift in the
$\mbox{$q\overline{q}$}\mbox{$\rightarrow$}Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$
peak position before smearing of $-0.405\pm 0.009$ GeV/c2 which shows the
experimental resolution increases the shift by an order of magnitude. The
larger number can be compared with the systematic error obtained by the
experiments on jet energy scale, which in the best measured regions are around
1% Aaboud:2017jcu ; Khachatryan:2016kdb  or even a little less. It can be seen
that the use of vector boson resonances as a standard candle for mass scale is
potentially problematic.

The analysis is repeated for only the subprocess
$\mbox{$b\overline{b}$}\mbox{$\rightarrow$}Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$
where large effects are expected as the initial and final state quarks are
identical. This finds a peak shift of $-0.81\pm 0.02$ GeV/c2, larger than any
shift reported in table 1. This demonstrates that the small interference in
the $b\overline{b}$ final state is because its initial state which is
dominated by the non-interfering lighter-quarks. Figure 3 (right) shows the
mass distribution for $b\overline{b}$$\rightarrow$$b\overline{b}$ after
smearing by 10%, where the change in the peak position is $-6.9\pm 0.4$
GeV/c2. Measuring this distribution would require knowledge of the initial
state quarks, which is not experimentally accessible.

####  3.1.1 Variation of procedure

This section discusses the stability of the results presented above as the
physics assumptions involved are varied.

The sensitivity of the peak position to the PDF uncertainties has been
explored using the 50 PDF eigenvector variations of MMHT2014nlo68cl.
111Attempting to use the PDF variations of NNPDF30NNLO caused Sherpa
termination for memory overflow on a 128 GB RAM computer. The quadratic sum of
the resulting shifts in the
$\mbox{$q\overline{q}$}\mbox{$\rightarrow$}Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$
peak position is 0.006 GeV/c2, corresponding to a 1.5% uncertainty on the
shift. PDF uncertainties are therefore considered negligible.

Variation of the pT threshold used to select quarks in the process generation
has also been studied. This was normally set to 25 GeV/c, but selections of 20
and 30 GeV/c have been investigated. The change in the pure electroweak cross-
section is $\pm$10%, but for the QCD it is a factor of two at 90 GeV/c2. The
effect is much larger for QCD because the $t$-channel production tends to have
large angular separation. Thus, lower pT thresholds enlarge the QCD sample in
the kinematic region where the colour structure can interfere. The result is
that the peak shift has quite a strong dependence on the pT selection, as seen
in table 2.

Quark pT | Inclusive | $|\Delta\eta|<1$ | $|\Delta\eta|>2$
---|---|---|---
threshold | cross-section | peak shift | cross-section | peak shift | cross-section | peak shift
20 GeV/c | 1882 pb |  $-0.465\pm 0.017$GeV/c2 | 954 pb |  $-0.280\pm 0.017$GeV/c2 | 283 pb |  $-0.933\pm 0.045$GeV/c2
25 GeV/c | 1741 pb |  $-0.405\pm 0.009$GeV/c2 | 945 pb |  $-0.309\pm 0.008$GeV/c2 | 156 pb |  $-0.737\pm 0.056$GeV/c2
30 GeV/c | 1564 pb |  $-0.351\pm 0.009$GeV/c2 | 945 pb |  $-0.302\pm 0.010$GeV/c2 | 0.2 pb | n.a.
Table 2: The electroweak cross-section accepted within a window of 89 to 93
GeV/c2 in $m_{\mbox{$u\overline{u}$}}$, selecting the $Z$ peak, as well as the
fitted shift in that peak due to interference. These are given for various
quark pT thresholds and diquark $\Delta\eta$ selections. A pT threshold of 30
GeV/c combined with a $|\Delta\eta|>2$ essentially excludes the $Z$ boson
peak, hence there is no entry for that region.

The dependence of the mass shift on the $|\Delta\eta|$ between the quarks is
also explored in table 2. It can be seen that for events where the quarks have
similar $\eta$, the pT threshold has little impact, while at large
$|\Delta\eta|$, where the pT can be lower, the threshold is much more
important. Thus, the results have a significant dependence on the kinematic
selection used.

The shift is also a function of the assumed resolution and a smearing width
varying from 0% to 13% is tested in figure 4. This is compared with the shift
extracted from convolving equation 2 with a Gaussian of variable width. The
match is good, apart from a small offset, until the width reaches 10%, where
the impact of the $t$-channel interference, not included in the equation,
becomes important.

Figure 4: The shift in the peak positions of Gaussian distributions fitted to
$u\overline{u}$ (left)and $b\overline{b}$ (right) peaks, plotted against the
fitted width. The squares show the impact of different smearing resolutions,
and the star symbols the effect of gluon radiation or Delphes simulation as
discussed in the text. The continuous lines show the impact of a numerical
smearing applied to the fitted interference formula.

The figure also shows the result of three progressively more complete ways of
introducing resolution: allowing a parton shower, allowing a parton shower and
an additional hard parton in the matrix element, and a full hadronization plus
the use of Delphes detector simulation, with the ATLAS detector simulation
employed. The resulting particles are all analysed using Delphes, with the
anti-kt Cacciari:2008gp  jet algorithm employed to order the resulting
particles. The events are required to have precisely two jets above the 25
GeV/c pT threshold and, in the case of the Delphes simulation, the generator
threshold is lowered to 15 GeV/c, with the reconstructed jets still required
to exceed 25 GeV/c2, to allow for threshold effects. The two jets are required
to be coplanar to 0.1 radians, which reduces events where a lot of energy is
missed from the jet. Finally, in the two cases where Delphes was not used, a
10% detector smearing is applied to the mass.

The resulting peak widths are 12.0 GeV/c2 for the parton shower, 12.1 GeV/c2
with additional hard radiation and 10.2$\pm$0.2 GeV/c2 for the Delphes
simulation. This happens to match almost exactly the width resulting from the
10% smearing applied to the quarks. The shift in the peak position is
3.91$\pm$0.28 GeV/c2 for the parton shower, 4.10$\pm$0.27 GeV/c2 for the
parton shower plus possible hard jet, and 2.42$\pm$0.42 GeV/c2 for the Delphes
events. These numbers can be visualised in figure 4 (left), which shows a
generally consistent pattern, with the Delphes results somewhat higher. This
has been traced to reduced acceptance for events from Delphes simulation at
large $|\Delta\eta|$, presumably because such events have a pT close to the
acceptance threshold and can easily fall below it as a result. This is the
region where the interference effects are larger, see table 2. 222These
Delphes studies in particular consumed many thousands of CPU hours and over a
hundred terabytes of disk.

The right hand side of figure 4 shows a comparison of no smearing, 10%
smearing and parton shower plus Delphes for the
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ case. A very similar pattern to
the $u\overline{u}$ is seen, with the effects scaled down by a factor of six.

###  3.2 Inclusive $W$ production

Inclusive $W$ production is analysed in a manner very similar to the $Z$, but
because the $W$ boson is charged it does not have any $s$-channel QCD
background. However, this background did not create interference in the $Z$
case and so the actual effects observed are similar. Sherpa is run in its
default form, with a diagonal CKM matrix, so the possible signatures are
$u\overline{d}$, $d\overline{u}$, $c\overline{s}$ and $s\overline{c}$. The
peak shifts and scale changes are given in table 3.

| $u\overline{d}$ | $d\overline{u}$ | $c\overline{s}$ | $s\overline{c}$
---|---|---|---|---
signal-to-background | 1.97 | 2.05 | 5.2 | 4.6
Shift, GeV/c2 | -0.287$\pm$0.005 | -0.242$\pm$0.006 | -0096$\pm$0.005 | -0.099$\pm$0.005
$\eta_{\mbox{scale}}$ | 1.018$\pm$0.003 | 1.027$\pm$0.004 | 0.995$\pm$0.002 | 1.039$\pm$0.003
Table 3: Interference effects on the
$\mbox{$q\overline{q}$}\mbox{$\rightarrow$}W$ peak by final state quark
flavour. Signal-to-background is defined at the resonance peak without
experimental resolution and using only the initial and final states quoted,
with negligible statistical error.

The shifts are found to be similar to, but smaller than, those in the $Z$
case, and the reduction of the effects for second-generation quarks is clear,
as was also the case for the $Z$ boson.

##  4 Boosted vector bosons

This section describes the dominant boosted production modes, $Vq$ and $Vg$.
These modes have different interference properties, but experimentally, they
are almost impossible to distinguish. Recoil against photons or other vector
bosons is considered later.

Reference Aad:2019fbh  used a pT selection of 600 GeV/c, but reference
PhysRevLett.120.071802  clearly shows that an inclusive hadronic $W$ peak
could be measured with a pT selection of 450 GeV/c. ATLAS and CMS used
unprescaled jet trigger thresholds of approximately 400 GeV/c in run 2 at the
LHC. In order to imitate the experiments, a pT threshold of 400 GeV/c is used
as the default here. The extraction of a hadronic $Z$ boson signal
experimentally often relies on $b$-quark identification to reduce QCD
backgrounds suppress the larger $W$ peak, so this is given particular
attention here.

###  4.1 Inclusive Boosted $W$

The processes requested for $W\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$ are
$\mbox{$q\overline{q}$}\mbox{$\rightarrow$}\mbox{$u\overline{d}$}g$,
$qg\mbox{$\rightarrow$}\mbox{$u\overline{d}$}q$ or
$\overline{q}g\mbox{$\rightarrow$}\mbox{$u\overline{d}$}\overline{q}$. It is
required that there is a $u\overline{d}$ quark pair with a mass of 50 to 150
GeV/c2, and as none of the diagrams allow for two identical quarks this is
unambiguous. In addition, at least one of the partons must have a pT above the
400 GeV/c threshold. In general, this means that the highest pT parton is
recoiling against the $u\overline{d}$ pair.

Figure 5: Left: example boosted $W$ boson plus quark production diagram.
Centre and Right: example QCD diagrams leading to the same quark state.

Sample Feynman diagrams are shown in figure 5. It is an amusing feature of the
QCD process in the centre that the colour label of the internal $u$ quark is
transferred first to the gluon and then to the final $d$ quark. Thus there is
a $u$ and $\overline{d}$ pair that is a colour singlet in both the electroweak
and this QCD diagram. This does not, however, apply to the right-hand diagram
as the incoming gluon changes the colour structure.

The cross-sections of the electroweak process ($W\rightarrow u\overline{d}$)
are 0.11 pb for the $Wg$ process and 0.23 pb for the $Wq$ process in a 2
GeV/c2 window around the $W$ peak. The signal-to-background ratios at the peak
of the resonance are very large, 115 and 681 respectively. This high signal-
to-background ratio implies that interference effects must be small. Indeed,
this is observed: the values of $\eta_{\mbox{scale}}$ are 1.004$\pm$0.002 for
the $Wg$ process and 1.000$\pm$0.002 for the $Wq$ process, and the
corresponding mass shifts are $-0.005\pm$0.004 GeV/c2 and 0.007$\pm$0.004
GeV/c2 respectively.

For the $Wg$ process, the interfering QCD diagrams have the same quark
flavours in the initial and final state, much like the inclusive production.
Therefore, interference effects will be smaller for
$W\mbox{$\rightarrow$}\mbox{$c\overline{s}$}$ than for
$W\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$, in the same way as described
for production at rest. Conversely, for the $Wq$ process, only one of the
final state quarks comes from the initial state and as a result, the
suppression will not be as strong. The effects are so small because the
$W$-like quark pair in QCD must come from different vertices.

The analysis has been repeated with a threshold of 200 GeV/c, where the
interfering QCD component is less suppressed. The signal-to-background ratios
fall by a factor of four and the resulting peak shift for
$Wg\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$ is 0.011$\pm$0.005 GeV/c2, and
for $Wq$ it is $-0.001\pm 0$.004 GeV/c2. Thus, even for a threshold far below
that currently employed for an unprescaled jet trigger, interference can be
safely neglected.

###  4.2 Inclusive Boosted $Z$

The inclusive boosted $Z$ production has many features in common with the $W$,
with possible recoil on a gluon, quark or antiquark considered. The
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ channel is studied here and it
is also compared with $Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$, which is
expected to show larger interference. The simulation is very similar to the
$W$ above, except that the $q\overline{q}$ pair of interest is either
$b\overline{b}$ or $u\overline{u}$. Unlike the $W$ case, two identical quarks
can be produced, and as the Sherpa selector acts on all matching combinations,
it is changed to require 50 GeV/c2 or more. If this had not been done, the
joint requirements that one of the partons has a pT above 400 GeV/c and the
upper mass cut on both $q\overline{q}$ pairs would together have eliminated
like-quark combinations.

Figure 6: Feynman diagrams showing the production of a $q\overline{q}$ pair
recoiling against either a quark or a gluon. The QCD examples drawn are
selected such that they can exhibit interference.

Sample Feynman diagrams are shown in figure 6. The rightmost diagram creates a
colour singlet $b\overline{b}$ final state without requiring a $b$ quark in
the initial state; there are other diagrams that produce a colour octet
$b\overline{b}$ from gluon splitting. These are important differences from the
$W$ case.

The signal-to-background ratios and peak shifts are given in table 4, where
results for the $W$ are also included. The signal-to-background ratio is one
to two orders of magnitude lower than for boosted $W$ production, so larger
shifts might be seen, but most of this comes from gluon splitting. The largest
shift is 0.022 GeV/c2, which after applying a 10% smearing produces a 0.18
GeV/c2 effect. This is negligible compared with the experimental systematic
errors.

Recoil |  pT, | Boson | signal / | $\eta_{\mbox{scale}}$ | Peak shift | Boson | signal / | $\eta_{\mbox{scale}}$ | Peak Shift
---|---|---|---|---|---|---|---|---|---
particle | GeV/c | mode | background |  | GeV/c2 | mode | background |  | GeV/c2
gluon | 400 | $Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$ | 9.3 | 1.004$\pm$0.002 |  $+0.004\pm$0.005 | $W\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$ | 116 | 1.000$\pm$0.002 | $-0.005\pm$0.004
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ | 21. | 0.996$\pm$0.002 |  $-0.009\pm$0.005
200 | $Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$ | 3.9 | 1.022$\pm$0.032 |  $+0.022\pm$0.005 | $W\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$ | 29 | 0.996$\pm$0.002 | $+0.011\pm$0.004
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ | 21. | 0.996$\pm$0.002 |  $+0.009\pm$0.004
quark | 400 | $Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$ | 1.05 | 1.016$\pm$0.009 |  $-0.004\pm$0.007 | $W\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$ | 680 | 1.018$\pm$0.002 | $0.007\pm$0.004
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ | 1.2 | 1.023$\pm$0.007 |  $-0.002\pm$0.007
200 | $Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ | 0.91 | 0.958$\pm$0.008 |  $+0.009\pm$0.004 | $W\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$ | 160 | 1.004$\pm$0.002 |  $-0.001\pm$0.004
Table 4: Interference effects on the $Z$ and $W$ boson peaks with an initial
state of either $qg$ plus $\overline{q}g$ or $q\overline{q}$ and separated by
final state quark flavour. Signal-to-background is defined at the resonance
peak without experimental resolution and using only the initial and final
states quoted, with negligible statistical error.

There is however one more thing which can potentially be measured in the
$\mbox{$b\overline{b}$}q$ state - the flavour of the outgoing quark. If a
third $b$ jet is identified, then diagrams like the one in figure 6 (right)
have enhanced possibilities for interference, as there are two $b\overline{b}$
combinations. Figure 7 shows the electroweak and $EW+I$ contributions when the
minimum jet pT requirement is reduced to 50 GeV/c.

Figure 7: The process $bg\rightarrow b\overline{b}b$. Left shows events where
$|\Delta\eta|<1$ between the $Z$ candidate quarks, while right requires
$|\Delta\eta|>2$. Each event generates two $b\overline{b}$ combinations and
two entries in these plots if they both fall within the mass range.

The $b\overline{b}$ pairs with similar pseudorapidity, $|\Delta\eta|<1$, show
a significant shift of $-0.17$ GeV/c2 at parton level, which will correspond
to a shift of more than a GeV/c2 after detector effects are considered, which
is similar to energy scale systematic errors. However, those with a large
$|\Delta\eta|$ separation show a dramatic effect, with a shift in the peak
position of around $+0.5$ GeV/c2. This would produce a shift of several GeV in
a measured spectrum. This effect might be observable at the LHC. For example,
the ATLAS collaboration in 2016 ATL-DAQ-PUB-2017-001  operated an unprescaled
trigger requiring two b-jets and a third jet with thresholds of 65, 65 and 110
GeV/c2. If an analysis could be performed in the trigger system, or prescaled,
it is conceivable that those thresholds could be halved.

##  5 $V$ \+ $\gamma$ production

Studying hadronic vector boson decays recoiling against a photon has
significant experimental advantages. Photons couple to charge, so the gluon
background is suppressed. ATLAS analysed
$\gamma\mbox{$b\overline{b}$}$Aad:2019wdr , with a trigger on a single photon,
an offline pT selection of 175 GeV/c, and a $Z$ candidate pT threshold of 200
GeV/c. For the current leading order quark-level analysis, the pT of the $V$
and the photon are exactly equal, and 200 GeV/c is used as a benchmark for
these processes. Sample Feynman diagrams are shown in figure 8, and results
are given in this section for $u\overline{u}$, $b\overline{b}$,
$u\overline{d}$ and $c\overline{s}$ quark combinations. As usual, the initial
state of $gg$, which does not interfere with the electroweak process at
leading order, is not included.

Figure 8: Example Feynman diagrams for
$\mbox{$q\overline{q}$}\mbox{$\rightarrow$}\mbox{$q\overline{q}$}\gamma$. Left
is an electroweak process and centre and right are referred to here as QCD.

As always, the colour structure of $s$-channel gluon production of the quark
pair means it does not contribute to interference, although its rate is
sizeable. The $t$-channel gluon exchange diagram can interfere, but it
preserves the quark species, so while it is important for
$\mbox{$u\overline{u}$}\gamma$, the small PDF suppress its contribution to the
$b\overline{b}$ state.

Figure 9: The cross-section for the process $q\overline{q}\rightarrow
u\overline{u}\gamma$, with the $u\overline{u}$ pair required to have pT$>$200
GeV/c. Left shows the pure electroweak (pale/green), pure QCD (dark/red) and
total (black) cross-sections. Right has the same electroweak component, but
the hollow black points show the total minus the QCD component, and the
dark/blue is the interference term.

The case of
$\mbox{$q\overline{q}$}\mbox{$\rightarrow$}Z\rightarrow\mbox{$u\overline{u}$}$
is shown in figure 9. The signal-to-background ratio, on peak, is nearly 3:1,
but the interfering fraction of the background, fitted as discussed in section
3, is only $0.7\pm 0.1$%. This is consistent with much of the background
arising from colour octet gluon splitting. The relative phase of the
background is reversed in comparison to the inclusive case, giving a positive
shift to the mass peak. The fitted shifts in the peak position are given in
table 5.

Recoil |  pT, | Boson | signal / | $\eta_{\mbox{scale}}$ | Peak shift | Boson | signal / | $\eta_{\mbox{scale}}$ | Peak Shift
---|---|---|---|---|---|---|---|---|---
particle | GeV/c | mode | background |  | GeV/c2 | mode | background |  | GeV/c2
photon | 400 | $Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$ | 5.2 |  $0.999\pm$0.002 |  $+0.016\pm$0.005 |  |  |  |
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ | 10.9 |  $0.997\pm$0.002 |  $-0.007\pm$0.004 |  |  |  |
200 | $Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$ | 2.7 |  $0.973\pm$0.003 |  $+0.045\pm$0.005 | $W\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$ | 35 |  $1.026\pm$0.002 | $+0.007\pm$0.004
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ | 9.7 |  $1.002\pm$0.002 |  $-0.008\pm$0.004 | $W\mbox{$\rightarrow$}\mbox{$c\overline{s}$}$ | 666 |  $0.999\pm$0.002 | $-0.001\pm$0.004
100 | $Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$ | 1.15 |  $0.918\pm$0.005 |  $+0.065\pm$0.007 |  |  |  |
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ | 8.1 |  $1.010\pm$0.002 |  $+0.001\pm$0.004 |  |  |  |
50 | $Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$ | 0.46 |  $0.769\pm$0.011 |  $-0.073\pm$0.014 | $W\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$ | 4.3 | $0.953\pm$0.003 | $-0.020\pm$0.005
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ | 5.5 |  $0.987\pm$0.002 |  $-0.004\pm$0.004
Table 5: Interference effects on the $Z$ and $W$ boson peaks in
$\mbox{$q\overline{q}$}\mbox{$\rightarrow$}V\gamma$, separated by final state
quark flavour. Signal-to-background is defined at the resonance peak without
experimental resolution and using only the initial and final states quoted,
with negligible statistical error.

The shift is 0.045 GeV/c2, which, when smeared by 10%, moves the peak position
by 0.26$\pm$0.06 GeV/c2. This is significantly below the roughly 1%
experimental systematic errors on the mass scale. The shift in
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ is only 0.008 GeV/c2, far below
any observable effects.

A range of pT selections is explored, from 400 GeV/c down to 50 GeV/c. The
magnitude of the shifts grow as the pT threshold is reduced. At 50 GeV/c the
$t$-channel electroweak exchange becomes important as evidenced by
$\eta_{\mbox{scale}}$ departing from 1. Figure 10 shows the differential
cross-sections for the sub-processes
$\mbox{$u\overline{u}$}\mbox{$\rightarrow$}\mbox{$u\overline{u}$}\gamma$ and
$\mbox{$d\overline{d}$}\mbox{$\rightarrow$}\mbox{$u\overline{u}$}\gamma$.
There is a qualitative difference in the two processes which highlights the
role of the $t$-channel $W$ exchange that is only present in the
$\mbox{$d\overline{d}$}\mbox{$\rightarrow$}\mbox{$u\overline{u}$}\gamma$
process; however, both channels show clear interference. With the 50 GeV/c pT
threshold requirement, the $t$-channel exchange components contribute mostly
for masses above the $Z$ mass, and so too does the interference, which
suppresses the high-mass cross-section by approximately 20%. It must be
stressed that experimentally these two processes are indistinguishable.

Figure 10: The differential cross-section for the process
$q\overline{q}\rightarrow u\overline{u}\gamma$, with the $u\overline{u}$ pair
required to have pT$>$50 GeV/c. Left is
$\mbox{$u\overline{u}$}\mbox{$\rightarrow$}\mbox{$u\overline{u}$}\gamma$ and
right is
$\mbox{$d\overline{d}$}\mbox{$\rightarrow$}\mbox{$u\overline{u}$}\gamma$. The
differential cross-section of pure electroweak is pale/green, pure QCD is
dark/red and the total is black. The electroweak $t$-channel component is much
more prominent when $t$-channel $W$ exchange is possible, but in each case the
destructive interference is clear.

Overall, the interference effects in the $u\overline{u}$ channel are small,
but in the experimentally accessed $b\overline{b}$ they are close to zero.
Again, the $b$ quarks PDF suppresses the interfering background.

The $W+\gamma$ channel is probably experimentally accessible, although there
is no published study of the peak. However, as seen in table 5, the signal-to-
background ratio is large, and the $s$-channel interference is vanishingly
small. In the same way as for the $Z+\gamma$ channel, at 50 GeV/c, a
$t$-channel electroweak component appears and complicates the picture.

##  6 Vector boson associated production

An interesting alternative source of hadronic vector bosons is diboson
production, with one of the bosons decaying leptonically. $WW$ and $ZZ$ decays
are explored here, requiring that the leptonic decay involves muons. Such
events are triggered efficiently by the LHC detectors, right down to zero
recoil. All leptons, charged or neutral, are required to have a pT above 25
GeV/c and the angular selections on the muons is $|\eta|<2.5$.

The background is similar to the previously-discussed modes with either
$q\overline{q}$ $s$-channel interaction via a colour-octet gluon state, or a
$t$-channel diagram in which, for $WW$, the final state $d\overline{u}$
inherits the flavours of the initial $u\overline{u}$ system. The two diagrams
are shown on the right-hand side of figure 11 and, within the limitations of a
diagonal CKM-matrix, only $u\overline{u}$ and $d\overline{d}$ initial states
can exhibit the $t$-channel process and thus interfere. The signal has many
diagrams, involving one $W^{+}$ along with $W^{-}$, $Z$ or $\gamma$
contributions, but the on-shell process sought has either a triple-boson
vertex or two separate $W$ boson emissions, as shown in figure 11 (left).

Figure 11: The leading order $q\overline{q}$$\rightarrow$
$W^{+}$$u\overline{d}$ diagrams. The first two diagrams show electroweak
processes giving rise to on-shell $WW$ pair production, and the second two are
the basic diagrams involving a gluon, referred to here as QCD .

The phase of the two electroweak diagrams shown is opposite, and their
relative importance depends upon the kinematics, as does the QCD contribution.
The inclusive shift in the peak position in $W^{+}$$u\overline{d}$ is
$0.111\pm 0.004$ GeV/c2, which is significantly smaller than that seen in the
prompt production. However, it is dependent on the event kinematics. This is
highlighted in figure 12, which shows (left) the relative phase space density
and (right) the shift in the mean mass. The shift is calculated here by
finding the mean mass of $d\overline{u}$ pairs within $\pm$2 GeV/c2 of 80
GeV/c2, and plotting the difference between the electroweak and the difference
between total and QCD as follows:

$\Delta m=\frac{\langle m_{total}\rangle\sigma_{total}-\langle
m_{QCD}\rangle\sigma_{QCD}}{\sigma_{total}-\sigma_{QCD}}-\langle
m_{EW}\rangle.$ (3)

This approach, calculating the deviation from the mean in a window,
understates any deviation from zero, but the bias has been found to be small.

Figure 12: The process $WW\rightarrow\mu^{+}\nu\mbox{$d\overline{u}$}$. Left
is the relative density of the electroweak signal as a function of the
$|\Delta\eta|$ and pT of the quark pair. Right is the shift in the mean mass
due to interference, as defined in the text.

The most probable kinematics has roughly the shift of $-0.1$ GeV/c2 found for
the mean, but when there is a pseudorapidity difference of 1.5 or more between
the quarks the shifts range from $-0.25$ GeV/c2 at rest to $+0.1$ GeV/c2 at
recoil momenta approaching 100 GeV/c. Such effects, corresponding to two or
three GeV/c2 after allowing for detector resolution, are potentially
observable.

pT, | Boson | signal / | $\eta_{\mbox{scale}}$ | Peak shift | Boson | signal / | $\eta_{\mbox{scale}}$ | Peak Shift
---|---|---|---|---|---|---|---|---
GeV/c | mode | background |  | GeV/c2 | mode | background |  | GeV/c2
100 | $ZZ\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$ | 6.9 | 0.981$\pm$0.002 |  $-0.005\pm$0.005 | $WW\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$ | 48 | 1.003$\pm$0.001 | $+0.001\pm$0.003
$ZZ\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ | 15.2 | 0.981$\pm$0.002 |  $+0.007\pm$0.005
0 | $ZZ\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$ | 4.0 | 0.885$\pm$0.002 | $-0.117\pm$0.005 | $WW\mbox{$\rightarrow$}\mbox{$u\overline{d}$}$ | 20.5 | 0.996$\pm$0.002 | $-0.080\pm$0.005
$WW\mbox{$\rightarrow$}\mbox{$d\overline{u}$}$ | 19.1 | 0.994$\pm$0.002 | $-0.085\pm$0.002
$ZZ\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ | 13.1 | 1.006$\pm$0.002 | $-0.010\pm$0.004 | $WW\mbox{$\rightarrow$}\mbox{$c\overline{s}$}$ | 86 | 1.011$\pm$0.002 | $-0.027\pm$0.004
$WW\mbox{$\rightarrow$}\mbox{$s\overline{c}$}$ | 92 | 0.997$\pm$0.002 | $-0.019\pm$0.004
Table 6: Interference effects on the hadronic $Z$ and $W$ boson peaks in
$\mbox{$q\overline{q}$}\mbox{$\rightarrow$}VV$ with one $V$ decaying
muonically and the other separated by final state quark flavour. Signal-to-
background is defined at the hadronic resonance peak without experimental
resolution and using only the initial and final states quoted, with negligible
statistical error.

The mean shifts observed in a set of diboson states are summarised in table 6.
$WW$ and $ZZ$ states are explored; $WZ$ are not expected to show qualitatively
different behaviour.

The experimental study of such a state would not be able, in general, to
distinguish jets from quarks or gluons. The
$W^{+}\mbox{$\rightarrow$}\mu^{+}\nu jj$ dijet mass spectrum is analysed
inclusively for the hadronic $W$ and $Z$ decays, and including all quark or
gluon initial state combinations. The 10% mass resolution is applied, and a
mass region from 70.775 GeV/c2 to 92.525 GeV/c2 is chosen, defined by the
masses where the binned electroweak signal differential cross-section has
halved compared with its peak. In that region, the pure electroweak cross-
section is 1.3 pb, while that for QCD production of the additional jets is 20
pb. This corresponds to production of 180,000 events in 140 fb-1 of LHC data,
or 700,000 if the electron mode and their charge conjugates are also
considered, which is may be sufficient to allow experimental study.

##  7 Discussion

The impact of electroweak-QCD interference on a variety of hadronic vector
boson production modes has been explored. The effect is mostly to move the
peak position without changing the integral. There can be destructive
interference effects between QCD and $t$-channel electroweak diagrams, for
example in figure 10, but they are not the main focus of this study.

The mass shifts are summarised in figure 13, where the left-hand plot shows
the production at rest, plotted against the signal-to-background ratio of the
quark final state in question. The different decay modes of the $W$ and $Z$
decay at rest are linked to highlight the correlation between large signal-to-
background ratio and small interference Brooijmans:2018xbu . Thus for the $Z$
boson the interference is largest for $u\overline{u}$ decay and smallest for
$b\overline{b}$. The $W$ boson does not have a gluon-splitting background, and
thus has higher signal-to-background than the $Z$, but the background is more
likely to have a colour singlet form, so the interference effects are larger
at a given signal-to-background. The peak shifts fall inside the range spanned
by the $Z$ decay modes.

Figure 13: Summary of the observed shifts in peak positions. Left shows the
peak shift as a function of the signal-to-background ratio in the final state,
with inclusive production and recoil against a $W$ or $Z$. For both the $W$
and $Z$ produced singly, a dashed line links the different decay modes. Right
shows hadronic vector boson recoiling against $q$, gluon, $\gamma$, $W$ or
$Z$, as a function of the pT of the vector boson. For the
$Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$, a dashed line links the
different pT selections.

Experiments generally uses $b$-tagging to identify $Z$ bosons. This suppresses
the $t$-channel QCD processes and reduces the interference effects rather
accidentally. If there were a pure and efficient way to identify $u$-quarks
available, experiments would use that to identify $Z$ candidates too, and
would then experience larger interference effects.

It may be possible to study inclusive
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ using LHC data. The ATLAS and
CMS collaborations record data at about 1 KHz, from an LHC bunch-crossing rate
of around 30 MHz. Thus one event in every 30,000 is recorded as pile-up in the
detectors, or approximately 4 pb-1. If a 10% resolution is assumed the
$Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$ cross-section in $80-100$ GeV/c2
at leading order is 2500 pb, and the QCD background in $b\overline{b}$ is only
8000 pb. If non-$b$ background can be suppressed using $b$-tagging, this
should be clearly visible.

In the light of the dependence on the signal-to-background ratio in 13 (left),
an explorative simulation of $t\overline{t}$ has been performed. This required
the process
$gg\rightarrow\mu^{-}\overline{\nu_{\mu}}\overline{b}\mbox{$d\overline{u}$}b$,
with selections on the two $W$ candidates masses between $50$ GeV/c2 and $120$
GeV/c2, and the top candidates masses between 160 GeV/c2 and 190 GeV/c2. The
cross-section with two electroweak vertices, one $W$ boson plus QCD jets, is
0.26$\pm$0.01 fb, while that at fourth order in electroweak, which includes
two on-shell top quarks, is 3.01$\pm 0.05$ pb. This signal-to-background
ratio, 104, strongly suggests that interference will have a completely
negligible impact.

The mass shifts calculated in various boosted states are shown against the
vector boson minimum pT requirement in figure 13 (right). As with the bosons
at rest, they are largest for $Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$,
and smallest for $Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$, with the $W$
bosons being intermediate. When vector bosons are boosted to 200 GeV/c or
more, there is generally a positive shift, but it is too small to be
experimentally observable. For transverse momenta of 50 to 100 GeV/c the
patterns are complicated, with the interference swapping sign, partly because
the electroweak process is no longer dominated by $s$-channel resonances. The
2D space of $|\Delta\eta$ between the quarks and the vector boson pT is
presented in figure 12 for the $WW$ case, where the varying contributions of
TGC production and emission of two independent $W$ bosons, which have opposite
phase, leads to opposite sign shifts in different regions.

The shifts in the peak positions are increased by roughly an order of
magnitude when smearing the results to mimic resolution effects. For example,
an intrinsic shift of $-0.409\pm 0.005$ GeV/c2 in
$Z\mbox{$\rightarrow$}\mbox{$u\overline{u}$}$ is enlarged to $-3.4\pm 0.2$
GeV/c2. This is reproduced by an approximate $s$-channel-only interference
formula. In this context, the parton shower simulation acts similarly to
detector resolution. This is not surprising as it is independently simulated
and knows nothing of the interference in the matrix elements. Simulation using
Delphes suggests that slightly smaller effects would occur as the jets with
large $|\Delta\eta|$ are less likely to be accepted. This could be offset by
working with a lower jet pT threshold than the 25 GeV/c generally used here.

For the original question that considers whether the boosted hadronic bosons
are good standard candles for the calibration of the detectors, the answer is
yes. The $W$ decays, especially boosted, have excellent intrinsic signal-to-
background ratio, while the $Z\mbox{$\rightarrow$}\mbox{$b\overline{b}$}$
decay mode is protected by the colour structure of QCD and the small $b$ PDF.
The impact of interference on kinematics close to the published measurements
has always been found to be small. The one caveat is that NLO effects might
give rise to colour-singlet $q\overline{q}$ pairs with a significant cross-
section, which would reopen the question for the $Z$ boson. There is no
analogous QCD process for $u\overline{d}$ production, so the $W$ boson results
should be regarded as more robust. Higgs to $b\overline{b}$ decay has a cross-
section significantly lower than $Z$, but its 4 MeV width improves the signal-
to-background ratio by 3 orders of magnitude and interference can have no
practical impact there.

In future trigger-level analyses, for example in the $W\gamma$ state with a 50
GeV/c pT threshold, the $s$\- and $t$-channel interference effects might
perhaps be possible to observe. More challenging to record is the
$b\mbox{$b\overline{b}$}$ final state, where large effects should be
observable if significant luminosity can be recorded with a pT threshold of 50
GeV/c. This $Z$ mode might be enhanced by NLO QCD corrections.

However, the $WW\mbox{$\rightarrow$}l\nu\mbox{$q\overline{q}$}$ process seems
to have the greatest promise for experimental study. The
$W\mbox{$\rightarrow$}\mbox{$d\overline{u}$}$ mass peak, after detector
resolution, is expected to change position by two or more GeV/c2 across the
kinematic plane. Smaller effects will occur for $c\overline{s}$, but it might
be possible to distinguish these final states using charm tagging. When only
irreducible background is considered, and pileup neglected, the signal-to-
background ratio exceeds 5%, so a more detailed study seems warranted. This
particular example includes a leptonically decaying boson, so the experiments
have recorded this data, and the experimental collaborations are encouraged to
see if the predicted effects can be observed.

###### Acknowledgements.

Great thanks are due to J. Quevillon who clarified many of these concepts for
us and suggested studies and cross-checks, L. Xia, for the helpful observation
that at NLO colour octet and singlet merge and to S. Schramm, for the
observation that pileup events contribute a significant trigger-free event
rate. We also wish to thank our journal referees, whose insistence that we
explore the impact of the jet pT threshold and the use of Delphes simulation
highlighted important issues.

## References

  * (1) G. Ranft, J. Ranft, Physics Letters B 87(1-2), 122 (1979).  DOI 10.1016/0370-2693(79)90034-0.  URL https://doi.org/10.1016/0370-2693(79)90034-0
  * (2) U. Baur, E. Glover, Physics Letters B 232(4), 519 (1989)
  * (3) J. Alitti, et al., Z. Phys. C 49(CERN-PPE-90-105), 17 (1990).  DOI 10.1007/BF01570793.  URL https://cds.cern.ch/record/210698
  * (4) J. Pumplin, Phys. Rev. D 45, 806 (1992).  DOI 10.1103/PhysRevD.45.806.  URL https://link.aps.org/doi/10.1103/PhysRevD.45.806
  * (5) U. Baur, Phys. Rev. D75, 013005 (2007).  DOI 10.1103/PhysRevD.75.013005.  URL https://link.apr.org/doi/10.1103/PhysRevD.75.013005
  * (6) J.M. Butterworth, A.R. Davison, M. Rubin, G.P. Salam, Phys. Rev. Lett. 100, 242001 (2008).  DOI 10.1103/PhysRevLett.100.242001.  URL https://link.aps.org/doi/10.1103/PhysRevLett.100.242001
  * (7) The ATLAS collaboration, Aad, G. et al.  Search for diboson resonances in hadronic final states in 139 fb-1 of $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector (2019).  URL https://arxiv.org/abs/1906.08589
  * (8) The CMS collaboration, Sirunyan, A. et al.  A multi-dimensional search for new heavy resonances decaying to boosted WW, WZ, or ZZ boson pairs in the dijet final state at 13 TeV (2019).  URL https://arxiv.org/abs/1906.05977
  * (9) The CMS collaboration, Sirunyan, A. et al., Phys. Rev. Lett. 120, 071802 (2018).  DOI 10.1103/PhysRevLett.120.071802.  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.071802
  * (10) The ATLAS collaboration, Search for boosted resonances decaying to two b-quarks and produced in association with a jet at $\sqrt{s}=13$ TeV with the ATLAS detector.  Tech. Rep. ATLAS-CONF-2018-052, CERN, Geneva (2018).  URL https://cds.cern.ch/record/2649081
  * (11) The ATLAS collaboration, Aad, G. et al.  Measurement of the jet mass in high transverse momentum $Z(\rightarrow b\overline{b})\gamma$ production at $\sqrt{s}=13$ TeV using the ATLAS detector (2019).  URL https://arxiv.org/abs/1907.07093
  * (12) The ATLAS collaboration, Aad, G. et al., Eur. Phys. J. C 79(arXiv:1906.11005. 10), 836. 54 p (2019).  DOI 10.1140/epjc/s10052-019-7335-x.  URL https://cds.cern.ch/record/2680245.  54 pages in total, author list starting page 38, 20 figures, 3 tables. All figures including auxiliary figures are available at http://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/PERF-2017-04
  * (13) M. Aaboud, et al., Phys. Lett. B786, 59 (2018).  DOI 10.1016/j.physletb.2018.09.013.  URL https://link.aps.org/doi/10.1016/j.physletb.2018.09.013
  * (14) The CMS collaboration, Sirunyan, A. et al., Phys. Rev. Lett. 121(12), 121801 (2018).  DOI 10.1103/PhysRevLett.121.121801.  URL https://link.aps.org/doi/10.1103/PhysRevLett.121.121801
  * (15) The CMS collaboration, Jet algorithms performance in 13 TeV data.  Tech. Rep. CMS-PAS-JME-16-003, CERN, Geneva (2017).  URL https://cds.cern.ch/record/2256875
  * (16) The ATLAS collaboration, Measurement of large radius jet mass reconstruction performance at $\sqrt{s}=8$ TeV using the ATLAS detector.  Tech. Rep. ATLAS-CONF-2016-008, CERN, Geneva (2016).  URL https://cds.cern.ch/record/2139642
  * (17) The ATLAS collaboration, Aad, G. et al., The European Physical Journal C 79(4), 290 (2019).  DOI 10.1140/epjc/s10052-019-6757-9.  URL https://doi.org/10.1140/epjc/s10052-019-6757-9
  * (18) The CMS collaboration, Sirunyan, A. et al., The European Physical Journal C 78(11), 891 (2018).  DOI 10.1140/epjc/s10052-018-6332-9.  URL https://doi.org/10.1140/epjc/s10052-018-6332-9
  * (19) T. Gleisberg, et al., JHEP 02, 007 (2009).  DOI 10.1088/1126-6708/2009/02/007.  URL http://www.slac.stanford.edu/spires/find/hep/www?eprint+0811.4622
  * (20) M. Cacciari, S. Frixione, N. Houdeau, M.L. Mangano, P. Nason, G. Ridolfi, JHEP 10, 137 (2012).  DOI 10.1007/JHEP10(2012)137.  URL https://link.aps.org/doi/10.1007/JHEP10(2012)137
  * (21) F. Cascioli, P. Maierhöfer, S. Pozzorini, Phys. Rev. Lett. 108, 111601 (2012).  DOI 10.1103/PhysRevLett.108.111601.  URL https://link.aps.org/doi/10.1103/PhysRevLett.108.111601
  * (22) M. Schönherr, F. Krauss, JHEP 12, 018 (2008).  DOI 10.1088/1126-6708/2008/12/018.  URL https://link.apr.org/doi/10.1088/1126-6708/2008/12/018
  * (23) R.D. Ball, et al., JHEP 04, 040 (2015).  DOI 10.1007/JHEP04(2015)040.  URL https://link.apr.org/doi/10.1007/JHEP04(2015)040
  * (24) L.A. Harland-Lang, A.D. Martin, P. Motylinski, R.S. Thorne, The European Physical Journal C 75(5), 204 (2015).  DOI 10.1140/epjc/s10052-015-3397-6.  URL https://doi.org/10.1140/epjc/s10052-015-3397-6
  * (25) T. Gleisberg, S. Höche, JHEP 12, 039 (2008).  DOI 10.1088/1126-6708/2008/12/039.  URL https://link.apr.org/doi/10.1088/1126-6708/2008/12/039
  * (26) The DELPHES 3 collaboration, J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, M. Selvaggi, Journal of High Energy Physics 2014(2), 57 (2014).  DOI 10.1007/JHEP02(2014)057.  URL https://doi.org/10.1007/JHEP02(2014)057
  * (27) R. Brun, F. Rademakers, Nucl. Instrum. Meth. A389, 81 (1997).  DOI 10.1016/S0168-9002(97)00048-X.  URL https://link.aps.org/doi/10.1016/S0168-9002(97)00048-X
  * (28) G. Brooijmans, et al., in _Les Houches 2017: Physics at TeV Colliders New Physics Working Group Report, Section 19, p145_ (2018).  URL http://lss.fnal.gov/archive/2017/conf/fermilab-conf-17-664-ppd.pdf
  * (29) S.P. Martin, Phys. Rev. D 86, 073016 (2012).  DOI 10.1103/PhysRevD.86.073016.  URL https://link.aps.org/doi/10.1103/PhysRevD.86.073016
  * (30) M. Aaboud, et al., Phys. Rev. D96(7), 072002 (2017).  DOI 10.1103/PhysRevD.96.072002.  URL https://link.aps.org/doi/10.1103/PhysRevD.96.072002
  * (31) The CMS collaboration, Khachatryan, V. et al., JINST 12(02), P02014 (2017).  DOI 10.1088/1748-0221/12/02/P02014.  URL https://link.aps.org/doi/10.1088/1748-0221/12/02/P02014
  * (32) M. Cacciari, G.P. Salam, G. Soyez, JHEP 04, 063 (2008).  DOI 10.1088/1126-6708/2008/04/063.  URL https://link.aps.org/doi/10.1088/1126-6708/2008/04/063
  * (33) The ATLAS collaboration, Trigger Menu in 2016.  Tech. Rep. ATL-DAQ-PUB-2017-001, CERN, Geneva (2017).  URL https://cds.cern.ch/record/2242069

##  Appendix A Example Run.dat

An example ‘run.dat’, for $\mu^{+}\nu\mbox{$d\overline{u}$}$ production with a
pT selection set to 0 GeV/c in the second last line.

    (run){
      # general settings
      EVENT_OUTPUT=HepMC_GenEvent[sb]
      EVENTS 25;

      # me generator setup
      ME_SIGNAL_GENERATOR Comix;
      SCALES VAR{Abs2(p[0]+p[1])};

      # Five lines added to switch off decays.
      SHOWER_GENERATOR=None
      FRAGMENTATION=Off
      MI_HANDLER=None
      ME_QED=Off
      BEAM_REMNANTS=0

      # LHC beam setup:
      BEAM_1 2212; BEAM_ENERGY_1 6500;
      BEAM_2 2212; BEAM_ENERGY_2 6500;
    }(run)

    (processes){
    Process 94 94 -> 1 -2 -13 14  ;
      Order (*,4);
      Print_Graphs graphs;
    End process;
    }(processes)

    (selector){
    ΨMass 1 -2  50 150
    ΨMass -13 14  70 90
    ΨPT 93 25 E_CMS
    ΨPT 90 25  E_CMS
    ΨPseudoRapidity 93 -2.5 2.5
    Ψ"Calc(PPerp(p[0]+p[1])>0)" 13,-14 1,1
    }(selector)