aainstitutetext: Rudolf Peierls Centre for Theoretical Physics, Clarendon
Laboratory, Parks Road, University of Oxford, Oxford OX1 3PU,
UKbbinstitutetext: CERN, Theoretical Physics Department, CH-1211 Geneva 23,
Switzerland

# Next-to-leading-logarithmic PanScales showers for Deep Inelastic Scattering
and Vector Boson Fusion

Melissa van Beekveld b Silvia Ferrario Ravasio
<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Received: date / Accepted: )

###### Abstract

We introduce the first family of parton showers that achieve next-to-leading
logarithmic (NLL) accuracy for processes involving a $t$-channel exchange of a
colour-singlet, and embed them in the PanScales framework. These showers are
applicable to processes such as deep inelastic scattering (DIS), vector boson
fusion (VBF), and vector boson scattering (VBS). We extensively test and
verify the NLL accuracy of the new showers at both fixed order and all orders
across a wide range of observables. We also introduce a generalisation of the
Cambridge-Aachen jet algorithm and formulate new DIS observables that exhibit
a simple resummation structure. The NLL showers are compared to a standard
transverse-momentum ordered dipole shower, serving as a proxy for the current
state-of-the-art leading-logarithmic showers available in public codes.
Depending on the observable, we find discrepancies at NLL of the order of
$15\%$. We also present some exploratory phenomenological results for Higgs
production in VBF. This work enables, for the first time, to resum
simultaneously global and non-global observables for the VBF process at NLL
accuracy.

###### Keywords:

QCD, Parton Shower, Resummation, LHC, HERA, DIS, VBF, VBS

††preprint: CERN-TH-2023-079, OUTP-23-05P

##  1 Introduction

Particle-physics collider experiments provide us with a unique opportunity to
test our knowledge of the fundamental interactions of elementary particles.
Accurately predicting signatures originating from Standard Model (SM) physics
is crucial to fully harness the potential of the data, and have sensitivity to
possible signals originating from beyond-the-Standard-Model scenarios.
Deciphering the nature of the Higgs boson and its interaction with other SM
particles is indeed one of the main objectives of the physics programme of the
Large Hadron Collider (LHC). General purpose Monte Carlo (GPMC) event
generators are fundamental tools in this context. They play a crucial role in
our understanding of the phenomenology of colliders, thanks to their ability
to describe much of the data from the LHC and its predecessors.

Parton showers lie at the core of GPMCs, and describe the energy degradation
of highly-energetic partons that are produced in the hard scattering process,
through radiation of soft and/or collinear partons. They enable us to simulate
arbitrarily complex collider events, characterised by a large multitude of
particles, whose modelling involves physics across a broad range of scales.
Despite their fundamental role in collider phenomenology, only in the recent
years more attention has been dedicated to understanding and improving the
formal accuracy of parton showers. Parton showers resum logarithmic terms in
the perturbative series. These terms arise from soft and collinear
divergences, and can become large when exploring physics across a wide range
of scales. Consequentially, terms of the form
$\alpha_{s}^{n}\sum_{i=1}^{n+1}c_{i}L^{i}$, where $L$ is a large logarithm of
a ratio of two disparate scales, will occur at each perturbative order $n$,
with $c_{i}$ a coefficient that depends on the observable. For example, with
the strong coupling $\alpha_{s}(m_{Z})=0.118$, a hard scale of $m_{Z}=91.18$
GeV, and a hadronic scale of $\Lambda=1$ GeV, the next-to-leading logarithmic
(NLL) corrections of the form $\alpha_{s}\ln(m_{Z}/\Lambda)$ become of the
order of $0.5$. This constitutes an $\mathcal{O}(1)$ correction, which needs
to be resummed to all orders in the perturbative coupling to obtain an
accurate prediction.

The Herwig7 Bahr:2008pv ; Bellm:2019zci  angular-ordered shower, which is
based on the coherent branching formalism Gieseke:2003rz , correctly resums
NLL terms (provided a careful interpretation of the ordering variable is
performed Bewick:2019rbu ; Bewick:2021nhc ) for global observables, but not
for non-global ones Banfi:2006gy . The NLL terms of the latter cannot be
captured by the coherent branching formalism, but require a dipole approach
Gustafson:1987rq , as they are sensitive to the full angular distribution of
soft emissions. Another practical advantage of a dipole shower, is that
matching with fixed-order matrix elements is much simpler. This has led to
numerous techniques to perform multi-jet merging at next-to-leading order
(NLO) Catani:2001cc ; Krauss:2002up ; Lavesson:2008ah ; Hoeche:2009rj ;
Giele:2011cb ; Platzer:2012bs ; Lonnblad:2012ix ; Frederix:2012ps ;
Lonnblad:2012ng ; Bellm:2017ktr ; Brooks:2020mab  , and matching to fixed-
order accuracies, i.e. (next-to-)next-to-next-to-leading order ((N)NNLO)
Hamilton:2012rf ; Alioli:2013hqa ; Hoche:2014dla ; Monni:2019whf ;
Campbell:2021svd ; Prestel:2021vww .111LO multi-jet merging techniques that
can be applied to an angular-ordered showers are discussed in Refs.
Mangano:2001xp ; Hamilton:2009ne ; Martinez:2021chk . Despite the large body
of work on improving the fixed-order accuracy of showers, until recently,
relatively little has been done on improving their logarithmic accuracy.
Indeed, all standard dipole showers currently embedded in public GPMC
Schumann:2007mg ; Platzer:2009jq ; Hoche:2015sya ; Cabouat:2017rzi ;
Brooks:2020upa ; Sjostrand:2006za ; Sjostrand:2014zea ; Gleisberg:2008ta ;
Sherpa:2019gpd ; Bierlich:2022pfr  are only leading logarithmic (LL), i.e.
they resum terms that are proportional to $\alpha_{s}^{n}L^{n+1}$.

The focus of this paper is to design the first dipole showers for processes
characterised by the $t$-channel exchange of a colour-singlet that reach NLL
(here defined to be single-logarithmic, $\alpha_{s}^{n}L^{n}$) accuracy for
both global and non-global observables. This work follows on earlier
developments of parton showers with a controlled logarithmic accuracy by the
PanScales collaboration, like those for dijet production in $e^{+}e^{-}$
collisions Dasgupta:2020fwr . Key components in the design of showers is a
careful construction of the recoil distribution after a new emission has been
generated, and its interplay with the ordering variable of the shower
Dasgupta:2018nvj . Several other groups have also been investigating NLL-
accurate showers in the context of $e^{+}e^{-}$ collisions. Ref.
Forshaw:2020wrq  introduced an algorithm that is shown to analytically
reproduce NLL accuracy for the thrust distribution and subject multiplicity.
The $\Lambda$-ordered Deductor shower is shown to be NLL accurate for the
thrust distribution in Ref. Nagy:2020dvz . Ref. Herren:2022jej  presents
Alaric, a shower that is proven to be NLL accurate for a wide range of global
observables. Designing NLL showers for hadronic colliders bring additional
complications with respect to that for $e^{+}e^{-}$ collisions, as the
treatment of the recoil from the initial-state is more subtle. In Refs.
vanBeekveld:2022zhl ; vanBeekveld:2022ukn , we presented NLL-accurate dipole
showers for the production of a colour singlet at the LHC. At present, no
other showers exists with demonstrated NLL accuracy (for both global and non-
global observables) for such processes.222A discussion on the treatment of the
transverse-momentum recoil for initial-state radiation in an angular-ordered
shower can be found in Ref. Bewick:2021nhc , as well as in the Herwig++ manual
Bahr:2008pv . This topic has also been addressed in Ref. Nagy:2009vg  in the
context of the transverse-momentum distribution of the $Z$ boson in Drell Yan
production, albeit without a claim on the logarithmic accuracy that is
achieved.  The present work extends the set of NLL accurate showers to
processes characterised by a $t$-channel exchange of a colour singlet, such
that now all processes with at most two partons at the first contributing
order can be described at single-logarithmic accuracy for generic observables.

(a)

(b)

(c)

Figure 1: Feynman diagrams contributing to (a) Higgs production in VBF, (b)
VBS and (c) DIS. Colour-connected quark lines are highlighted in red and blue.
Leptons are represented with black solid lines.

Examples of processes with a $t$-channel exchange of a colour-singlet, are
vector-boson fusion (VBF, Fig. 1(a)) and vector-boson scattering (VBS, Fig.
1(b)). Measurements of VBF and VBS processes can provide valuable insight into
the electroweak (EW) and Higgs sectors BuarqueFranzosi:2021wrv . Higgs
production in VBF is the second most abundant production mode for the Higgs
boson at the LHC. It has a clear experimental signature, given by the presence
of two back-to-back hard jets in the forward and backward regions of the
detectors, and hence it is the preferred channel to for example measure the
Higgs coupling to muons ATLAS:2020fzp ; CMS:2020xwi , tauons CMS:2022kdi ;
ATLAS:2022yrq , and to study Higgs to invisible decays ATLAS:2022yvh ;
CMS:2022qva . Run III data will enable us to perform precise determinations of
the size of gauge-bosons self-interactions. Modifications to VBS processes are
predicted in models of physics beyond the SM, through changes to the Higgs
boson couplings to gauge bosons and the resonant production of new particles.
Differences in the parton-shower modelling turn out to be one of the leading
uncertainty in these kind of processes Ballestrero:2018anz ; Jager:2020hkz ;
Hoche:2021mkv ; CMS:2020xwi ; ATLAS:2022tnm ; ATLAS:2022ooq ; ATLAS:2022yrq ;
ATLAS:2020rej ; CMS:2022kdi ; CMS:2021kom ; CMS:2022uhn ; Buckley:2021gfw .
Having access to several NLL-accurate showers, such as the ones we present in
this article, will help to assess more realistic shower uncertainties, and
possibly reduce them.

Higher-order QCD contributions to VBF and VBS processes are computed using a
factorised approach Han:1992hr ; Bolzoni:2010xr ; Bolzoni:2011cu , i.e.
neglecting non-factorisable corrections stemming from the exchange of partons
between the two hadronic sectors, which are coloured in red and blue in Figs.
1(a) and 1(b). Non-factorisable corrections appear only from NNLO and are
typically colour-suppressed with respect to the factorisable ones Liu:2019tuy
; Dreyer:2020urf . In the factorised approximation, radiative corrections to
VBF and VBS are closely related to those for deep inelastic scattering (DIS,
Fig. 1(c)).333 VBF can be treated using the structure function approach, i.e.
treating the two quark lines as two seperate copies of DIS. Instead, for VBS,
interferences between diagrams with a $t$-channel exchange of a boson and the
tagged bosons being emitted from the quark legs directly spoil this picture.
From the QCD point of view, one can still use the factorised approximation and
neglect gluon exchanges between the two quark lines, also see Ref.
BuarqueFranzosi:2021wrv  and references therein. However, for the VBS
topology, one needs to decide whether a boson needs to be boosted together
with all partons belonging to one hadronic sector, to preserve its angle with
respect to the original-final state quark, which is sensible when that boson
was emitted directly from the quark line. Otherwise, the shower may reshuffle
the post-branching momenta in such a way that the boson becomes collinear to
the jet, or that a collinear boson is suddenly emitted at a wide angle: in
both situations, the matrix element used to generate the process would then no
longer describe the post-branching situation with the new four-momenta.  For
this reason, the first part of this article focuses on the formulation of the
PanScales showers for DIS processes. We then extend our showers to handle two
copies of DIS, so that they can also be used to describe VBF and VBS events in
the factorised approach. After including mass effects, which we leave for
future work, this framework will also enable us to handle $t$-channel single-
top production. These showers will have full colour accuracy at LL, and
leading colour accuracy at NLL, i.e. single-logarithmic in our definition.
This is the first time this logarithmic accuracy is reached, either
analytically or numerically, for the VBF and VBS processes.

In addition to its relevance to hadron-collider processes involving a
colourless $t$-channel exchange, an improved theoretical framework for deep
inelastic scattering (DIS) can be directly applied to interpret data gathered
from electron-proton ($ep$) colliders like the Hadron Electron Ring
Accelerator (HERA). These colliders provide an ideal environment for
investigating the internal structure of hadrons and conducting accurate
studies of quantum chromodynamics (QCD). The new-generation of lepton-hadron
colliders, such as the Electron-Ion Collider (EIC), whose construction is
planned to start next year, would also benefit from more accurate MC
generators.

The paper is structured as follows. In Sec. 2 we review the standard kinematic
variables used in DIS. In Sec. 3 we present the common building-blocks of a
generic dipole-shower for DIS. In Sec. 3.1, we summarise the main features of
a standard dipole shower, while in Sec. 3.2 we present the new PanScales
showers for DIS, whose extension to VBF/VBS is discussed in Sec. 3.3. We then
perform fixed-order tests of such showers in Sec. 4. The all-orders
validations of the PanScales showers are instead performed in Sec. 5. We test
numerous observables, designed to probe soft and/or collinear emissions that
an NLL shower should accurately describe: the DGLAP evolution of the parton
distribution functions (Sec. 5.1), average particle multiplicity (Sec. 5.2),
several continuously-global observables (Sec. 5.3), and the amount of
radiation in a rapidity slice (Sec. 5.4), which we consider as a proxy for a
generic non-global observable. In Sec. 6 we show some phenomenological results
for Higgs production in VBF, and in Sec. 7 we present our conclusions.

##  2 DIS definition and kinematics

We consider the DIS process $\gamma^{*}p\to X$, where $p$ is the incoming
proton, $\gamma^{*}$ is a space-like photon (or more generically a colourless
boson) and $X$ is a generic hadronic final-state. We denote with
$q_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}^{\mu}$
the photon momentum, and we define

$Q^{2}=-q_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}^{2}>0\,.$
(1)

The (massless) proton momentum is denoted by $P^{\mu}$. We introduce two
light-like reference vectors $n_{1}$ and $n_{2}$, such that $n_{1}$ is
parallel to the incoming proton i.e.

$n_{1}^{\mu}=x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}P^{\mu}\,,$
(2)

and

$n_{2}^{\mu}=q_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}^{\mu}+n_{1}^{\mu}\,.$
(3)

Requiring that $n_{2}^{2}=0$ leads to

$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}=\frac{Q^{2}}{2\,q_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}\cdot
P}\,.$ (4)

Furthermore, we see that $n_{1}\cdot n_{2}=Q^{2}/2$. At the partonic level,
the lowest-order contribution reads $\gamma^{*}\,q\to q$, where the incoming
(outgoing) quark momentum is precisely $n_{1}^{\mu}$ ($n_{2}^{\mu}$). More
generally, if the sum of the momenta of the final-state partons has an
invariant mass
$M_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}$,
one may always parameterise their collective momentum
$p_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}$
as

$p_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}=n_{2}^{\mu}+\frac{M_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}{Q^{2}}n_{1}^{\mu}\,.$
(5)

Momentum conservation then implies that the incoming parton has momentum

$p_{1}^{\mu}=\left(1+\frac{M_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}{Q^{2}}\right)n_{1}^{\mu}=x_{1}P^{\mu}\,,\quad\mbox{with
}\quad
x_{1}=x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}\left(1+\frac{M^{2}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}}{Q^{2}}\right).$
(6)

Indeed, with these definitions we have

$p_{1}^{\mu}+q_{{\rm
dis}}^{\mu}=p_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}\,.$
(7)

We work in the Breit frame Webber:1993bm , which is defined as the frame where
$n_{1}$ and $n_{2}$ are back-to-back and aligned along the $z$-axis, i.e.

$\displaystyle n_{1}^{\mu}=\frac{Q}{2}\left(0,0,-1;1\right),$
$\displaystyle\qquad n_{2}^{\mu}=\frac{Q}{2}\left(0,0,+1;1\right),$ (8a)
$\displaystyle
q_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}^{\mu}$
$\displaystyle=Q\left(0,0,+1;0\right),$ (8b)

with $Q\equiv\sqrt{Q^{2}}$. The final-state momenta can then always be
decomposed in terms of Sudakov variables, i.e.

$\displaystyle k_{i}^{\mu}=\alpha_{i}n_{1}^{\mu}+\beta_{i}n_{2}^{\mu}+k_{\perp
i}^{\mu}\,,$ (9)

where $k_{\perp i}^{\mu}$ is a generic space-like vector orthogonal to
$n_{1,2}^{\mu}$. Partons are said to reside in the current hemisphere when
$\beta_{i}>\alpha_{i}$. The remnant hemisphere instead contains partons with
$\alpha_{i}>\beta_{i}$. When constructing the shower, it will turn out to be
convenient to introduce the reference vector

$Q^{\mu}=n_{1}^{\mu}+n_{2}^{\mu}\,,$ (10)

which has norm
$Q^{2}=-q_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}^{2}$.
We note that the energy component in the Breit frame of a given momentum
$p^{\mu}$ is obtained through

$E_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
Breit}=\frac{p\cdot Q}{Q}\,.$ (11)

We will use the reference vector $Q^{\mu}$ in the formulation of the PanScales
showers to measure angular distances in the Breit frame instead of the dipole
frame.

##  3 Dipole showers for DIS and VBF/VBS

The fundamental building block for a dipole shower is a $2\to 3$ branching
kernel Gustafson:1987rq . In these types of showers each emitter
${\tilde{\imath}}$ is colour-connected (understood in the limit of a large
number of colours $N_{\text{{c}}}$) to a spectator ${\tilde{\jmath}}$, such
that the branching is ${\tilde{\imath}}{\tilde{\jmath}}\to ijk$, with $k$ the
radiated parton. Dipole showers for DIS (and VBF or VBS) need to handle
emissions that come from dipoles that contain two final-state partons (FF
dipoles), or one parton in the initial state and one in the final state (IF
dipoles).

First, we consider an FF dipole with pre-splitting momenta $\tilde{p}_{i}$,
$\tilde{p}_{j}$, and post-splitting momenta $p_{i}$, $p_{j}$, $p_{k}$. The
momentum of the radiated parton $k$ may be decomposed as

$p_{k}^{\mu}=z_{i}\tilde{p}_{i}^{\mu}+z_{j}\tilde{p}_{j}^{\mu}+k_{\perp}^{\mu}\,.$
(12)

where $k_{\perp}^{\mu}$ is a space-like vector orthogonal to
$\tilde{p}_{i,j}^{\mu}$. In this notation with $\tilde{p}_{i}$ identified as
emitter, $z_{i}$ corresponds to the collinear momentum fraction carried away
by the emission $p_{k}$, defined relative to the pre-branching momentum. The
emission probability that describes correctly radiation in the soft and
collinear limit can be written as

${\rm d}\mathcal{P}_{{\tilde{\imath}}{\tilde{\jmath}}\to
ijk}=\frac{\alpha_{s}(\mu_{\text{{r}}}^{2})}{2\pi}\left(1+\frac{\alpha_{s}(\mu_{\text{{r}}}^{2})K}{2\pi}\right)\frac{{\rm
d}v^{2}}{v^{2}}{\rm d}\bar{\eta}\frac{{\rm
d}\phi}{2\pi}\left[g(\bar{\eta})z_{i}P^{\rm{FS}}_{ik}(z_{i})+g(-\bar{\eta})z_{j}P^{\rm{FS}}_{jk}(z_{j})\right],$
(13)

where $\mu_{\text{{r}}}^{2}$ is the renormalisation scale. The phase space of
the emission is parameterised by the shower variables $v$ and $\bar{\eta}$,
where $v$ can be identified with the shower ordering scale, and $\bar{\eta}$
is a rapidity-like auxiliary variable. The exact relation between
$z_{i},(z_{j}),k_{\perp}$ and the shower variables $v,\bar{\eta}$ is shower
dependent and will be detailed in the following sections. The azimuthal angle
of the radiation defined with respect to the plane spanned by the two pre-
branching dipole momenta is denoted by $\phi$. To achieve NLL accuracy one
needs to evaluate the running of the coupling at two loops, with $\mu_{{\rm
R}}$ equal to the transverse-momentum of the emission in the soft-collinear
limit. In addition, the soft-collinear gluon emission probability must include
an $\alpha_{s}K/(2\pi)$ correction term with
$K=(67/18-\pi^{2}/6)C_{A}-5n_{\\!\>\\!f}/9$ Catani:1990rr . We take $C_{A}=3$,
$T_{R}=1/2$ and we work with $n_{\\!\>\\!f}=5$ light flavours (unless
otherwise stated). The definition of the DGLAP final-state splitting functions
$P^{\rm{FS}}(z)$ that we use is given in appendix A of Ref.
vanBeekveld:2022zhl , with their soft limit obtained as $z\to 0$. Following
Ref. Dasgupta:2020fwr , we use two definitions of $g$:

$g(\bar{\eta})=g^{\text{dip.}}(\bar{\eta})\equiv\begin{cases}0&\text{ if
}\bar{\eta}<-1\,,\\\
\displaystyle\frac{15}{16}\left(\frac{\bar{\eta}^{5}}{5}-\frac{2\bar{\eta}^{3}}{3}+\bar{\eta}+\frac{8}{15}\right)&\text{
if }-1<\bar{\eta}<1\,,\\\ 1&\text{ if }\bar{\eta}>1\,,\end{cases}$ (14)

or

$g(\bar{\eta})=g^{\text{ant.}}(\bar{\eta})\equiv\frac{e^{\bar{\eta}}}{e^{\bar{\eta}}+e^{-\bar{\eta}}}=\frac{e^{2\bar{\eta}}}{e^{2\bar{\eta}}+1}\,.$
(15)

We use $g^{\text{ant.}}$ for antenna showers, where one makes no distinction
between an emitter/spectator to distribute the transverse-momentum recoil,
while we use $g^{\text{dip.}}$ when only one of the parent partons takes the
transverse recoil.

Instead, for IF dipoles where we define the parton ${\tilde{\imath}}$ as
belonging to the initial-state, we write the momentum of the radiated parton
as

$p_{k}^{\mu}=\frac{z_{i}}{1-z_{i}}\tilde{p}^{\mu}_{i}+z_{j}\tilde{p}_{j}^{\mu}+k^{\mu}_{\perp}\,.$
(16)

The emission probability then takes the form

${\rm d}\mathcal{P}_{{\tilde{\imath}}{\tilde{\jmath}}\to
ijk}=\frac{\alpha_{s}(\mu_{\text{{r}}}^{2})}{2\pi}\left(1+\frac{\alpha_{s}(\mu_{\text{{r}}}^{2})K}{2\pi}\right)\frac{{\rm
d}v^{2}}{v^{2}}{\rm d}\bar{\eta}\frac{{\rm d}\phi}{2\pi}\\\
\times\left[\frac{x_{i}f_{i}(x_{i},\mu_{\text{{f}}}^{2})}{\tilde{x}_{i}f_{{\tilde{\imath}}}(\tilde{x}_{i},\mu_{\text{{f}}}^{2})}g(\bar{\eta})z_{i}P^{\rm{IS}}_{ik}(z_{i})+g(-\bar{\eta})z_{j}P^{\rm{FS}}_{jk}(z_{j})\right],$
(17)

where $f_{i}$ is the PDF of the incoming parton $i$, and
$x_{i}=\frac{\tilde{x}_{i}}{1-z_{i}}$. The choice of the factorisation scale
$\mu_{\text{{f}}}$ depends on the ordering variable $v$, and it must be chosen
such that for a hard-collinear initial-state branching, $\mu_{\text{{f}}}$ is
equal to the transverse momentum of the emission. The definition of the DGLAP
initial-state splitting functions is also given in Appendix A of Ref.
vanBeekveld:2022zhl .

The DIS invariants
$q^{2}_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}$
and
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}$,
introduced in the previous section, determine the structure of the event. At
LO these two invariants are the only two quantities needed to describe the
interaction between the incoming lepton and nucleon. Although unitarity of the
shower is preserved, the real-virtual cancelation on more differential LO
observables is spoiled when modifying the DIS invariants through the
generation of an emission. Preserving the LO DIS invariants is therefore
crucial to ensure that the shower does not alter the description of any
inclusive observable. Furthermore, this also simplifies the inclusion of
higher-order corrections in future works, especially in the context of VBF if
the two hadronic sectors evolve completely independently. For these reasons,
the conservation of the DIS invariants
$q^{2}_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}$
and
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}$
is a fundamental property of all the showers we consider/develop here.

###  3.1 A standard transverse-momentum-ordered dipole shower

In this section, we briefly summarise the fundamental features of a standard
dipole shower, which we use as a proxy to illustrate the behaviour of current
standard leading-logarithmic showers that are used for DIS and VBF/VBS
phenomenological studies. We refer to this shower as “Dipole-$k_{t}$”, and its
kinematic maps (with local momentum conservation in the IF and FF dipoles) are
presented in Refs. vanBeekveld:2022zhl ; vanBeekveld:2022ukn .444The ordering
variable and kinematic maps coincide with those of the Dire-v1 shower
Hoche:2015sya . Here we review the basic characteristics of the Dipole-$k_{t}$
shower.

The ordering variable $v$ of the Dipole-$k_{t}$ shower algorithm is
transverse-momentum-like, i.e. it corresponds to the transverse momentum of
the emitted particle in the limit that this emission is soft-collinear. The
rapidity-like variable $\bar{\eta}$, used to partition the dipole in two
halves, can be related to the collinear momentum fraction $z$. This relation
reads

$\bar{\eta}=\begin{cases}\frac{1}{2}\ln\frac{z^{2}\,\tilde{s}_{ij}}{v^{2}}\,&\text{(final
state),}\\\
\frac{1}{2}\ln\frac{z^{2}\tilde{s}_{ij}}{(1-z)^{2}v^{2}}&\text{(initial
state)},\end{cases}$ (18)

with $\tilde{s}_{ij}=2\tilde{p}_{i}\cdot\tilde{p}_{j}$ the dipole invariant
mass. From this it is clear that $\bar{\eta}$ partitions the dipole in its
rest frame. In the Dipole-$k_{t}$ shower, momentum conservation is fully
local. This implies that when an emission occurs from a FF dipole, the
transverse momentum recoil is entirely taken up by the emitter, which
corresponds to the original dipole leg closer in angle to the emission.
Conversely, in an IF dipole, the transverse momentum recoil is always absorbed
by the final-state leg. The factorisation and renormalisation scales are set
equal to $v$. Further details on the kinematic maps used in the Dipole-$k_{t}$
shower are given in Appendix B.1 of Ref. vanBeekveld:2022zhl .555Conversely to
Refs. vanBeekveld:2022zhl ; vanBeekveld:2022ukn , here we use only the local
variant of Dipole-$k_{t}$ and we do not consider the global variant. The
latter algorithm does not preserve the DIS invariants (specifically
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}$
and $y_{\rm
dis}=\frac{Q^{2}}{x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}s}$)
and for this reason it has never been used in phenomenological applications to
DIS or VBF/VBS.

###  3.2 PanScales showers for DIS

The PanScales showers need a reference momentum $Q^{\mu}$ to define a common
frame where to measure angular distances for all the emissions, which here we
take to be the Breit frame. We set this equal to $Q^{\mu}$ as introduced in
eq. (10), which in the Breit frame reads $Q^{\mu}=(0,0,0;Q)$. The ordering
variable $v$ is defined in such a way that for a soft-collinear emission of
transverse momentum $k_{\perp}$ and rapidity $\eta$, we have

$v\approx
k_{\perp}e^{-\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}|\eta|}\,.$ (19)

We also introduce the shower variable
$\bar{\eta}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}Q}$,
that corresponds to the rapidity of a soft-collinear emission in the Breit
frame. Given the shower variables $v$ and
$\bar{\eta}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}Q}$,
we can define a transverse-momentum auxiliary variable

$\displaystyle\kappa_{\perp}\equiv\rho
ve^{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}|\bar{\eta}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}Q}|}\,,$
(20)

where we have used

$\displaystyle\tilde{s}_{i}=2\tilde{p}_{i}\cdot
Q\,,\quad\tilde{s}_{j}=2\tilde{p}_{j}\cdot
Q\,,\quad\tilde{s}_{ij}=2\tilde{p}_{i}\cdot\tilde{p}_{i}\,,\quad\rho=\left(\frac{\tilde{s}_{i}\tilde{s}_{j}}{\tilde{s}_{ij}Q^{2}}\right)^{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}/2}\,.$ (21)

To achieve NLL accuracy, the renormalisation scale at which the coupling
constant is evaluated is set to $\kappa_{\perp}$, while for the factorisation
scale we choose

$\mu_{\text{{f}}}=Q\left(\frac{v}{Q}\right)^{\frac{1}{1+\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}}}\,.$ (22)

We also introduce the variables

$\displaystyle\alpha_{k}\equiv\sqrt{\frac{\tilde{s}_{i}}{\tilde{s}_{j}\tilde{s}_{ij}}}\kappa_{\perp}e^{\bar{\eta}_{Q}}\,,\qquad\beta_{k}\equiv\sqrt{\frac{\tilde{s}_{j}}{\tilde{s}_{i}\tilde{s}_{ij}}}\kappa_{\perp}e^{-\bar{\eta}_{Q}},$
(23)

which enable us to write the light-cone momentum fraction at which we need to
evaluate the DGLAP splitting probabilities appearing in eqs. (13), (17). If
the original-dipole legs ${\tilde{\imath}}$ and ${\tilde{\jmath}}$ are final-
state partons, we define

$z_{i}=\alpha_{k}\,,\qquad z_{j}=\beta_{k}\,,$ (24)

while if ${\tilde{\imath}}$ is an incoming parton, the definition of $z_{i}$
is modified and we instead use

$z_{i}=\frac{\alpha_{k}}{1-\alpha_{k}}\,.$ (25)

####  3.2.1 PanGlobal

This section details the kinematic mapping of an antenna shower with global
transverse-momentum recoil. We refer to this shower as PanGlobal. The
treatment of the longitudinal recoil is identical to that of the PanGlobal
variant for hadron collisions of Refs. vanBeekveld:2022zhl ;
vanBeekveld:2022ukn . However, the DIS variant differs in the boost that is
performed to achieve momentum conservation in the perpendicular component.
This choice is motivated by the fact that in our case we want to preserve
$q_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}^{\mu}$,
i.e. the momentum of the $t$-channel exchanged boson, while for colour-singlet
production the most natural variables to preserve are the invariant mass and
the rapidity of the colour singlet system. One begins by introducing some
intermediate post-branching dipole momenta, which read

$\displaystyle\bar{p}_{i}^{\mu}$ $\displaystyle=(1\pm
a_{k})\tilde{p}_{i}^{\mu}\,,$ (26a) $\displaystyle\bar{p}_{j}^{\mu}$
$\displaystyle=(1\pm b_{k})\tilde{p}_{j}^{\mu}\,,$ (26b)
$\displaystyle\bar{p}_{k}^{\mu}$
$\displaystyle=a_{k}\tilde{p}_{i}^{\mu}+b_{k}\tilde{p}_{k}^{\mu}+k_{\perp}^{\mu}\,,$
(26c)

with $a_{k}=\alpha_{k}$, $b_{k}=\beta_{k}$ and

$k_{\perp}^{\mu}=\sqrt{a_{k}b_{k}\tilde{s}_{ij}}\left(\hat{k}_{\perp
1}^{\mu}\sin\phi+\hat{k}_{\perp 2}^{\mu}\cos\phi\right),$ (27)

with $\hat{k}_{\perp 1,2}$ two vectors with norm $-1$ orthogonal to
$\tilde{p}_{i,j}$ such that $\hat{k}_{\perp 1}\cdot\hat{k}_{\perp 2}=0$. The
signs in eq. (26) depend on whether the ${\tilde{\imath}}$ or
${\tilde{\jmath}}$ is incoming ($+$) or outgoing ($-$).

The momentum mapping of eq. (26) clearly does not conserve momentum. In both
cases the final state has received a perpendicular component. This is not
accounted for in the initial-state momentum sum, $\bar{p}_{1}^{\mu}+q_{\rm
dis}^{\mu}$ with $\bar{p}_{1}^{\mu}$ the momentum of the initial-state parton
after radiation. The imbalance for both the IF and FF kinematic maps is given
by

$\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}-\bar{p}_{1}^{\mu}=q_{\rm
dis}^{\mu}+k_{\perp}\,.$ (28)

To see this, note that for FF dipoles we get the imbalance

$\displaystyle\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}$
$\displaystyle\equiv\tilde{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}-\tilde{p}_{i}^{\mu}+\bar{p}_{i}^{\mu}-\tilde{p}_{j}^{\mu}+\bar{p}_{j}^{\mu}+\bar{p}_{k}^{\mu}$
(29) $\displaystyle=q_{\rm dis}^{\mu}+\bar{p}_{1}^{\mu}+k_{\perp}^{\mu}\,,$

with
$\tilde{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}$
the pre-branching sum of final-state momenta, and
$\bar{p}_{1}^{\mu}\equiv\tilde{p}_{1}^{\mu}$ the momentum of the initial-state
parton. Note that we have used eq. (7) to rewrite
$\tilde{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}$.
For IF dipoles we instead get

$\displaystyle\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}$
$\displaystyle=\tilde{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}-\tilde{p}_{j}^{\mu}+\bar{p}_{j}^{\mu}+\bar{p}_{k}^{\mu}\,,$
(30) $\displaystyle=q_{\rm dis}^{\mu}+\bar{p}_{1}^{\mu}+k_{\perp}^{\mu}\,,$

where $\bar{p}_{1}^{\mu}\equiv(1+a_{k})\tilde{p}_{1}^{\mu}$, see eq. (26).

To restore momentum conservation, we first address the requirement that the
system composed of the exchanged boson and the incoming parton have a squared
invariant mass consistent with that of the sum over all final-state momenta,
i.e.

$(p_{1}+q_{\rm
dis})^{2}=\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}\,.$
(31)

To this aim, we reset the momentum of the incoming quark $p_{1}^{\mu}$.
Requiring that $p_{1}^{\mu}$ only has a component in the direction of
$n_{1}^{\mu}$ gives us

$\displaystyle p_{1}^{\mu}$
$\displaystyle=\left(1+\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}{Q^{2}}\right)n_{1}^{\mu}=\left(\frac{Q^{2}+\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}{Q^{2}+\tilde{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}\right)\tilde{p}_{1}^{\mu}\,.$
(32)

In the second step we boost the partonic final-state momenta from
$\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}$
to

$p_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}=n_{2}^{\mu}+\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}{Q^{2}}n_{1}^{\mu}\,.$
(33)

There is considerable freedom in constructing the boost. We design our boost
such that the transverse momentum recoil is smoothly redistributed to the
partons in the current hemisphere, i.e. partons that are collinear to the
direction of $n_{2}^{\mu}$. Remnant-hemisphere partons are instead largely
unaffected in their perpendicular components, and only receive a rescaling in
the component along the direction of $n_{1}^{\mu}$. The boost that achieves
this is derived in Appendix A and reads

$\displaystyle\Lambda^{\mu\nu}=$ $\displaystyle
g^{\mu\nu}+\frac{2n_{1}^{\mu}}{Q^{2}}\left[(\beta-1)n_{2}^{\nu}+\frac{p_{t}^{2}}{\beta
Q^{2}}n_{1}^{\nu}+p_{\perp}^{\nu}\right]+\frac{2n_{2}^{\mu}n_{1}^{\nu}}{Q^{2}}\frac{1-\beta}{\beta}-\frac{2p_{\perp}^{\mu}n_{1}^{\nu}}{\beta
Q^{2}}\,.$ (34)

The components $\beta$ and $p_{\perp}$ are defined through a Sudakov
decomposition of the final-state momentum sum

$\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}=\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}+p_{t}^{2}}{\beta
Q^{2}}n_{1}+\beta n_{2}+p_{\perp}\,,$ (35)

with $p_{\perp}$ a space-like vector orthogonal to $n_{1}^{\mu}$,
$n_{2}^{\mu}$ with norm $-p_{t}^{2}$. Notice that for soft-collinear
emissions, $1-\beta\sim p_{t}\sim\kappa_{t}$.

####  3.2.2 PanLocal

In this section, we describe how to implement a dipole shower with local
recoil for DIS, which we refer to as PanLocal.666We omit the description of an
antenna PanLocal version, as future applications such as NLO matching or
multi-jet merging would be substantially more cumbersome, as found in Ref.
Hamilton:2023dwb . For all three dipole types (IF, FI and FF), we parameterise
the momentum of the radiated parton as

$\bar{p}_{k}^{\mu}=a_{k}\tilde{p}_{i}^{\mu}+b_{k}\tilde{p}_{k}^{\mu}+k_{\perp}^{\mu}\,.$
(36)

The emitter can either be an initial-state parton (for IF dipoles), or a
final-state one (for FI and FF dipoles). In the former case we use

$a_{k}=\alpha_{k},\qquad
b_{k}=\beta_{k}(1+\alpha_{k})^{\frac{2}{1+\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}}}.$ (37)

This choice was introduced in Ref. vanBeekveld:2022zhl  to restore transverse-
momentum ordering for very hard-collinear emissions, and serves to avoid
unphysical correlations between emissions in opposite hemispheres. Like in
Refs. vanBeekveld:2022zhl ; vanBeekveld:2022ukn , the new momenta of the
emitter and the spectator become

$\displaystyle\bar{p}_{i}^{\mu}=$
$\displaystyle(1+a_{k})\tilde{p}_{i}^{\mu}+\frac{a_{k}b_{k}}{1+a_{k}}\tilde{p}_{j}^{\mu}+k_{\perp}^{\mu}\,,$
(38) $\displaystyle\bar{p}_{j}^{\mu}=$
$\displaystyle\left(1-\frac{b_{k}}{1+a_{k}}\right)\tilde{p}_{j}^{\mu}\,.$ (39)

Although the momentum is locally conserved (i.e.
$\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}-\bar{p}_{1}^{\mu}=q_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}^{\mu}$),
we now end up in a situation where the incoming parton is no longer aligned
with the beam direction $n_{1}^{\mu}$. The Lorentz transformation that we
apply to realign the incoming parton with the beam differs from the one
applied for colour-singlet production, since in this case we want to preserve
the DIS invariants. After this transformation, the momentum of the incoming
parton becomes

$p_{i}^{\mu}=p_{1}^{\mu}=\left(1+\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}{Q^{2}}\right)n_{1}^{\mu}\,,$
(40)

and the sum of final-state partons

$p_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}=n_{2}^{\mu}+\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}{Q^{2}}n_{1}^{\mu}\,.$
(41)

To achieve this, we first rotate all momenta except the photon momentum
$q_{\rm dis}^{\mu}$ with a rotation matrix $R(\theta)$. This matrix is defined
to align the post-branching incoming parton momentum $\bar{p}_{1}$ along the
direction of $n_{1}^{\mu}$. This operation introduces a momentum imbalance
since we do not change the photon momentum $q_{\rm dis}^{\mu}$. To restore
momentum conservation we proceed in the same way as for PanGlobal, i.e. we
Sudakov decompose the sum of (rotated) final-state momenta as

$\displaystyle(R(\theta)\cdot\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X})^{\mu}=\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}+p_{t}^{2}}{\beta
Q^{2}}n_{1}^{\mu}+\beta n_{2}^{\mu}+p_{\perp}^{\mu}\,,$ (42)

and apply the boost of eq. (34) to all partons including the initial-state
one.

In the case where the emitter is instead a final-state parton we set

$a_{k}=\alpha_{k},\qquad b_{k}=\beta_{k}\,.$ (43)

The momenta of the emitter and spectator in the case of an FF/FI dipole reads

$\displaystyle\bar{p}_{i}^{\mu}=$
$\displaystyle(1-a_{k})\tilde{p}_{i}^{\mu}+\frac{a_{k}b_{k}}{1-a_{k}}\tilde{p}_{j}^{\mu}-k_{\perp}^{\mu}\,,$
(44) $\displaystyle\bar{p}_{j}^{\mu}=$
$\displaystyle\left(1\pm\frac{b_{k}}{1-a_{k}}\right)\tilde{p}_{j}^{\mu}\,,$
(45)

where the $+$($-$) sign is required for FI(FF) dipoles. Note that no
boost/rotation needs to be performed in these cases since the initial-state
parton does not acquire a transverse momentum component. Hence,
$\bar{p}_{i,j,k}^{\mu}=p_{i,j,k}^{\mu}$ for FI and FF dipoles.

###  3.3 Extension to VBF and VBS

From the point of view of QCD radiative corrections, VBF (Fig. 1(a)) and VBS
(Fig. 1(b)) can be seen as a double copy of DIS (Fig. 1(c)). Non-factorisable
corrections arising from the exchange of partons between the two hadronic
sectors are colour suppressed and only contribute from order $\alpha_{s}^{2}$.
They were computed for the first time in Ref. Liu:2019tuy  for Higgs
production in VBF, using the eikonal approximation, and found to be typically
ten times smaller than the NNLO factorisable corrections Cruz-Martinez:2018rod
; Dreyer:2018rfu ; Asteriadis:2021gpd , as confirmed by the phenomenological
study of Ref. Dreyer:2020urf , which also addresses double Higgs production.

For this reason, to shower a VBF or VBS process, we treat the two hadronic
sectors as two separate and independent DIS processes. This is done by
labelling the two sectors $a$ and $b$, and using for each of them a separate
reference vector $Q^{\mu}_{a}$ and $Q^{\mu}_{b}$ that reads

$\displaystyle Q^{\mu}_{a}=n_{1,a}^{\mu}+n_{2,a}^{\mu}\,,$ (46)

where $n_{1,a}^{\mu},n_{2,a}^{\mu}$ are the light-like reference vectors
introduced in Sec. 2 such that
$q^{\mu}_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS,a}=n_{2,a}^{\mu}-n_{1,a}^{\mu}$.
Here
$q^{\mu}_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS,a}$
describes the four-momentum of the $t$-channel colour-singlet boson that is
exchanged between the initial- and final-state partons. A similar definition
is employed for sector $b$, with different light-like reference vectors
$n_{1,b}^{\mu}$ and $n_{2,b}^{\mu}$. Since the showering of the two hadronic
sectors is completely factorised, we have the freedom to choose different
values for the shower starting scale for the two different sectors, as well as
the factorisation and renormalisation scales employed in the underlying fixed-
order calculation.

Once we have established SL accuracy for the showers in the DIS process,
extending this accuracy to VBF topologies is conceptually straightforward
(especially in the case of VBF), with the exception of non-factorisable
corrections, which are not included.777We remind the reader that interferences
between the diagrams where the tagged bosons in the VBS process are connected
to the $t$-channel propagator (Fig. 1(b)), and those where the tagged bosons
in the VBS process are emitted directly from the quark legs spoil the
structure-function picture, so higher-order corrections to the VBS process are
not simply those for DIS squared.  The terms that we are neglecting are
colour-suppressed NLL corrections. Describing correctly the factorisable
contributions is however sufficient to obtain full-colour accuracy at LL, and
NLL accuracy in the large $N_{\text{{c}}}$ limit. Additionally, non-
factorisable contributions are typically suppressed after applying VBF cuts
Dreyer:2020urf , and also unknown in their complete form (only approximated
results have been obtained at fixed order in Ref. Liu:2019tuy ).

##  4 Fixed-order tests

One fundamental property that NLL-accurate parton showers must satisfy is that
emissions widely separated in rapidity or transverse momentum must be
independent Dasgupta:2018nvj . This requirement follows from the factorisation
properties of the underlying QCD matrix element in the limit of soft and/or
collinear emissions. It can be directly translated into a “fixed-order”
criterion: “a soft emission can alter the momentum of a previously radiated
parton only if they are very close in rapidity and emitted with a commensurate
transverse momentum”.

To better understand this requirement, consider an emission with momentum
$\tilde{k}_{1}^{\mu}$ from a dipole. After this, the shower will try to
generate a second emission, with momentum $k_{2}^{\mu}$. The redistribution of
momenta (either directly in the kinematic map or via a momentum-conserving
boost) may alter the kinematics of the first emission, i.e.
$\tilde{k}_{1}^{\mu}\to k_{1}^{\mu}$. The result is that the event now has two
emissions with momenta $k_{1}^{\mu}$ and $k_{2}^{\mu}$, but the first emission
was generated with a matrix element corresponding to $\tilde{k}_{1}^{\mu}$.
The matrix element corresponding to $k_{1}^{\mu}$ and $\tilde{k}_{1}^{\mu}$
are not the same when these momenta are vastly different, hence this would
destroy the factorisation property it must obey when $k_{2}^{\mu}$ is either
collinear to any of the partons in the pre-branching event or very soft. This
process gets repeated for every emission that follows, leading to a wrong
logarithmic exponentiation starting at NLL (at leading colour accuracy). To
verify the condition above, in Sec. 4.1 we investigate the behaviour of the
showers introduced in the previous section in the presence of two gluons
emitted at commensurate values of the ordering variable $v$.888Note that this
is a necessary, but not sufficient requirement to reach NLL accuracy for
global observables. For the latter, it is important that any pair of emissions
are described accurately in the shower, while in the fixed-order test of Sec.
4.1 only the action of the second emission on the first emission’s momentum is
probed.

Another fixed-order criterion that we test is that subleading-colour
corrections are correctly implemented in case of strongly-ordered emissions.
Indeed, dipole showers are implemented in the large-$N_{\text{{c}}}$ limit,
but subleading (i.e. $\propto 1/N_{\text{{c}}}^{2}$) corrections in the LL
contribution have the same numerical size of NLL terms, thus from a
phenomenological point of view they have the same relevance and must be
included. For this reason, in Sec. 4.2 we considered two strongly-ordered
emissions and we assess whether subleading-colour corrections have been
correctly incorporated in the shower at fixed order through a comparison to
the exact analytic matrix element. Note that this algorithm is expected to not
only yield the correct result at LL full-colour, but also at NLL for many
observables, such as average particle multiplicity (Sec. 5.2) and global
observables (Sec. 5.3), as well as for the DGLAP evolution (Sec. 5.1).

The final ingredient necessary to reach NLL is the implementation of spin
correlations to correctly reproduce the azimuthal distributions of the
radiation. To this aim, in appendix D we show that the algorithm introduced in
Ref. vanBeekveld:2022zhl  for PanScales showers for hadron collisions can be
applied also to DIS-type processes. However, since spin correlations do not
impact the NLL resummation of the observables we use to validate the NLL
accuracy of our showers at all orders (detailed in Sec. 5), we do not include
them in our all-order tests.

###  4.1 Lund-plane contours

In this section, we investigate if emissions of commensurate hardness are
independent when widely separated in angle, as required to achieve NLL
accuracy. We consider the DIS process $q\,\gamma^{*}\to q$, with
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}=0.01$,
and we emit a gluon $g_{1}$ from the $qq$ dipole at a fixed shower evolution
scale of $\ln v_{1}/Q=-20$, azimuthal angle $\phi_{1}=0$, and a few fixed
values of $\bar{\eta}_{Q,1}$. We then examine the impact of a second gluon
$g_{2}$, emitted at a scale $\ln v_{2}/Q=\ln v_{1}/Q-2$, on the kinematics of
$g_{1}$, varying $\bar{\eta}_{Q,2}$. Its azimuthal angle is fixed at
$\phi_{2}=0$. We parameterise the available phase space of the emissions in
terms of the Lund variables $\ln k_{t,i}/Q$ and $\eta_{i}$ Andersson:1988gp .
Some care has to be taken in defining these variables, in particular for
emissions off the final-state quark. If only one gluon emission is present, we
calculate two angular distances, one that parameterises the distance from the
beam
($d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}kB}$),
and one to the final-state quark
($d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}kF}$)

$d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}kB}=1-\cos\theta_{k}\,,\qquad
d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}kF}=1-\cos\theta_{kq}\,.$
(47)

where $\theta_{k}$ is the angle between $k$ and the beam, while $\theta_{kq}$
is the angle between $k$ and the final-state quark $q$, both defined in the
Breit frame. We then define the angle

$\cos\theta\equiv\begin{cases}+\cos\theta_{kq}\qquad\mbox{ if
}d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}kF}<d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}kB}\,,\\\
-\cos\theta_{k}\qquad\mbox{ if
}d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}kF}>d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}kB}\,,\end{cases}$
(48)

so that the Lund variables read

$\eta=\frac{1}{2}\ln\frac{1+\cos\theta}{1-\cos\theta}\,,\qquad\ln\frac{k_{t}}{Q}=\frac{1}{2}\ln\frac{E_{k}^{2}(1-\cos^{2}\theta)}{Q^{2}}\,.$
(49)

The sign in eq. (48) is chosen such that partons in the current hemisphere
have positive rapidity, while partons in the remnant hemisphere have negative
rapidity. Now let us move to the case where two gluons are present. In this
case, we first calculate the angular distances of eq. (47) for both gluons.
Focusing on emissions on the primary Lund plane, there are two possible
scenarios: the smallest angular distance can either be to the beam, or the
final-state quark. Note that the third scenario, where the second gluon is
closest to the first emitted gluon, is instead described by a secondary Lund
plane, which is not considered here. If the smallest angular distance is
$d_{k_{i}B}$, we simply use the definition of the Lund variables of eq. (49)
for both gluon emissions. If instead the smallest angular distance is between
any of the two gluons and the final-state quark, we merge them into a new
‘final-state’ momentum $p_{F}^{\mu}$. We then recalculate
$d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}kF}$
using that momentum, and calculate the new Lund variables for the last
recombined emission.

Notice that this definition of the Lund variables shares some similarities
with the variables associated with the generalisation of the Cambridge/Aachen
algorithm for DIS, detailed in Appendix B. The main differences are first that
for the current case, we are only interested in describing the ‘primary’ Lund
plane, ignoring the case where the two emissions would first cluster together,
as we want to focus on the case in which the emissions are widely separated in
angle. Secondly, to define the transverse momentum of the splitting we always
use the energy of the radiated gluon, while in the flavour-blind algorithm of
Appendix B we take the energy of the softest parton.

Figure 2: Double-emission contours for Dipole-$k_{t}$ (top) and
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$) (bottom) for a first soft-collinear emission in the remnant hemisphere
(left) and in the current hemisphere (right) for a DIS process with $x_{\rm
DIS}=0.01$ and an incoming proton with negative rapidity in the Breit frame.
The phase-space contours are shown as a function of $\ln k_{t}/Q$ and $\eta$.
A red dot indicates the kinematics of the first emission which is fixed at
$\eta_{1}=-7$ (left) or $\eta_{1}=7$ (right) and $\ln v_{1}/Q=-20$ (in the
plot labels, values of $v_{i}$ are always expressed in units of $Q$). The
contour of the first emission, obtained by fixing $\ln v_{1}/Q$ but varying
$\eta_{1}$, is shown with the mostly horizontal red dotted line. That of the
second emission at $\ln v_{2}/Q=-22$ is drawn as a blue solid line. The colour
shading of the lines indicates the branching probability. In the bottom
panels, we show the logarithm of the ratio between the transverse momentum of
the first emission after ($k_{t,1}$) and before the second emission took place
($\tilde{k}_{t,1}$). This is expected to be zero (dashed grey line) except
when the two emissions are close in rapidity (grey vertical band, where the
rapidity of the first emission is indicated with a vertical red dotted line).

On the top two panels of Fig. 2 we illustrate the two-emission contours for
Dipole-$k_{t}$, a transverse-momentum ordered shower. We notice that the first
emission erroneously takes the transverse-momentum recoil if

$\eta_{2}<\frac{1}{2}\left(\eta_{1}-\ln\frac{v_{1}}{Q}\right).$ (50)

This is because this shower uses a fully local map, such that the first gluon
always absorbs the recoil whenever the second emission comes from the new
initial-final $(q_{I}g_{1})$ dipole in the region $\eta_{2}<\eta_{1}$, and
because the midpoint of the final-final dipole is assigned in the dipole frame
in the region where $\eta_{2}>\eta_{1}$.

This behaviour can be corrected by either choosing a different evolution
variable, or by conserving the transverse momentum globally, as done by the
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$) shower. As discussed in Sec. 3.2.1, the global boost that takes care
of redistributing the transverse momentum imbalance has the property of
affecting mostly those partons at very large (and positive) rapidity, while
leaving the ones living in the remnant hemisphere unchanged. Indeed, the
bottom two panels of Fig. 2 show that subsequent emissions widely separated in
rapidity leave the Lund variables associated with the first emission
unaffected. The results obtained after taking a different ordering variable
are shown in Fig. 3 for
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$) and PanLocal. Here the ordering variable is set to $v\approx
k_{t}e^{-\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}|\eta|}$, with
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$. As already discussed in Ref. Dasgupta:2020fwr  local transverse-
momentum conservation requires one to choose
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}>0$.

Figure 3:  Same as Fig. 2, but for
PanLocal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$, left) and
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$, right) for $\eta_{1}=-8.7$.

Finally in Fig. 4 we show results when choosing a different value for
$\eta_{1}$ resulting in a hard-collinear first emission. They are analogous to
the case where the first emission is soft-collinear, but note that care has to
be taken for the PanLocal shower. At variance with the
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$) case, for very collinear initial-state radiation the ordering
variable for PanLocal behaves like a transverse momentum. As discussed in Ref.
vanBeekveld:2022zhl  and in Sec. 3.2.2 around eq. (37), this modification in
the local shower is necessary to prevent a very hard-collinear emission from
significantly impacting the momentum of a softer previous emission.

We conclude that in the DIS PanScales showers, conversely to standard dipole
showers, emissions widely separated in angle leave the kinematics of the
previous emissions untouched. The novelty of our showers is that the
transverse-momentum recoil due to initial-state-radiation is effectively
absorbed by hard partons in the current hemisphere, while soft partons or
partons collinear to the initial state remain unaffected, as required from
colour coherence.

Figure 4:  Same as Fig. 3, but considering a hard-collinear first emission
with $\eta_{1}=12$.

###  4.2 Subleading-colour corrections

Dipole showers are developed using the large-$N_{\text{{c}}}$ approximation of
QCD, and designed to correctly describe not only collinear emissions, but also
soft wide-angle gluon emissions in the large-$N_{\text{{c}}}$ limit. In
standard dipole showers, subleading-colour corrections are included by
replacing $C_{A}/2=N_{\text{{c}}}/2$ with
$C_{F}=(N_{\text{{c}}}^{2}-1)/(2N_{\text{{c}}})$ when a quark leg is
identified as emitter (we refer to this as the colour-factor-from-emitter
scheme, CFFE). It has been known for quite some time Gustafson:1987rq  that
this choice is inconsistent with colour coherence and more recently
Dasgupta:2018nvj  it was observed that it leads to wrong (subleading-colour)
contributions already at LL, due to the incorrect assignment of the emitter.
Since $1/N_{\text{{c}}}^{2}\sim\alpha_{s}$, these LL mistakes have the same
size of NLL terms, so it stands to reason that these subleading colour terms
should be included in an NLL-accurate shower.

Two schemes that do result in full-colour accuracy for LL terms were
introduced in Ref. Hamilton:2020rcu  for showers applicable to $e^{+}e^{-}$
collisions, which also have been generalised for colour-singlet production at
hadron colliders in Ref. vanBeekveld:2022zhl : the so-called segment and NODS
schemes. These schemes furthermore result in correct full-colour accuracy at
NLL in the case of global observables, and next-to-double-logarithmic (NDL)
accuracy for jet multiplicity, the latter of which is sensitive to the
integrated rate of double-soft energy-ordered emissions at commensurate
angles. The segment colour scheme divides each dipole into an arbitrary number
of distinct segments, identified through the generation variable
$\bar{\eta}_{Q}$. These segments either have a $C_{F}$ or $C_{A}/2$ colour
factor, and are assigned respecting colour coherence. The nested-ordered-
double-soft (NODS) scheme was designed to not only get the correct integrated
rate of two soft and energy-ordered emissions that occur at commensurate
angles, but also describe them correctly at the differential level. This is
achieved by applying a matrix element correction that describes a pair of
energy-ordered commensurate-angle emissions, well separated in rapidity from
all other emissions. For more details on these schemes, see Refs.
Hamilton:2020rcu ; vanBeekveld:2022zhl . These algorithms are
straightforwardly extended to the DIS case.

Figure 5: Density for the emission of an ultra-soft gluon from a
$q\gamma^{*}\rightarrow qg_{1}$ configuration (left), or from a
$g\gamma^{*}\rightarrow q\bar{q}$ system (right). In both cases, the incoming
parton has a negative rapidity in the Breit frame, the final-state quark a
positive one, and the first-emitted final-state parton ($g_{1}$ in the left
pane, $\bar{q}$ in the right one) has a rapidity $\eta_{1}=-5$. The emission
density is illustrated as a function of the Lund variables $\eta$ and $\psi$
(defined in the main text). For each configuration, the left panel corresponds
to the primary Lund plane, while the right panel represents the secondary Lund
plane (populated by emissions from the parton located at $\eta_{1}=-5$, see
Appendix B). From top to bottom, the three rows show the results for the
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$) shower result with the segment, NODS and CFFE colour scheme
implemented For each row, the upper panel illustrates the ratio between the
parton shower differential cross section and the leading colour (LC) result,
$\;\mathrm{d}\sigma_{\rm PS}/\;\mathrm{d}\sigma_{\rm LC}$, obtained setting
$C_{A}=2C_{F}=3$, while the lower panels show the deviation from the full
colour (FC) differential matrix element, $\;\mathrm{d}\sigma_{\rm
PS}/\;\mathrm{d}\sigma_{\rm FC}$.

To validate these schemes for DIS, we perform tests of the differential matrix
element produced by the shower after two strongly-ordered emissions. The
kinematics of the first emission (which can either be a quark or a gluon) is
fixed at $\eta_{1}=-5$, $\phi_{1}=0$ and $\ln v_{1}/Q=-10$. The second
emission is then emitted at a fixed value for $\ln v_{2}/Q=-60$, and sampled
over $\bar{\eta}_{Q,2}$ and $\phi_{2}$. The resulting shower predictions are
then compared with the analytic result at full colour (FC).

In Fig. 5 we show the results obtained with
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$) shower employing the segment, NODS and CFFE colour schemes. The
emission density on the primary and secondary Lund plane are shown in terms of
the rapidity $\eta$ of the second emission, and its azimuthal angle $\psi$.
These variables are defined according to our generalised DIS Cambridge-Aachen
algorithm (Appendix B). Specifically, for emissions in the primary Lund plane,
we define $\psi=\phi_{2}-\phi_{1}+\pi$ (or $\psi=\phi_{2}-\phi_{1}-\pi$, if
$\phi_{2}-\phi_{1}$ is larger than $\pi$). For emissions in the secondary Lund
plane we instead use

$\psi=\arctan\frac{\phi_{2}-\phi_{1}}{\eta_{2}-\eta_{1}}+\pi\,.$ (51)

Let us first focus on the ratio between the parton shower and the analytic LC
result, obtained using $C_{A}=2C_{F}=3$. These are illustrated in the upper
set of panels (coloured green-blue). For the majority of phase space, the
value of this ratio is either $1$ or $8/9$, when the effective colour factor
is either $C_{A}/2$ or $C_{F}$. The lower set of panels (coloured red-white-
blue) shows the ratio between the parton shower and the analytic FC result. As
expected, the segment scheme assigns the correct colour factor everywhere
except for the region where the second gluon is close in angle to the first
emission. Note however that these deviations integrate to $0$ after averaging
over the angular phase space, therefore the all-order validations in Sec. 5
would not be sensitive to this effect. The NODS colour scheme corrects for
this effect and produces the analytic differential matrix element also in the
region where a pair of energy-ordered emissions occur at commensurate angles.
The last set of plots shows the deviation from the CFFE scheme, which is
currently implemented in standard dipole showers. This results in a wrong
colour assignment over a region in phase space that is logarithmically
extended. Hence this will result in a wrong LL term at FC, as already pointed
out in Ref. Dasgupta:2018nvj .

Note that although here we have only shown results for the
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$) shower, analogous results can be obtained for
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$) and PanLocal. Similar considerations also apply to the
Dipole-$k_{t}$ shower. However, the fraction of phase space that has the wrong
subleading colour corrections in the CFFE scheme is larger, as the emitter is
chosen by partitioning the dipole frame, instead of in the Breit frame.

##  5 All-order validation

In this section, we investigate the all-order behaviour of the PanScales
showers for DIS over a broad range of observables, targeting distinct classes
of next-to-leading logarithmic effects, i.e. those of the form
$\alpha_{s}^{n}L^{n}$. In Sec. 5.1 we test the capability of the showers to
reproduce the DGLAP evolution of the parton distribution functions at single-
logarithmic (SL) accuracy, which probes nested emissions in the hard-collinear
region. In Sec. 5.2 we focus on the average particle multiplicity at next-to-
double-logarithmic accuracy (NDL), targeting nested emissions in both the soft
and collinear regions. In Sec. 5.3 we compare the shower cumulative cross
section against the NLL predictions for several continuously-global
observables, which is the only test sensitive to both double- and single-
logarithmic terms in the Sudakov exponentiation (including running coupling
effects), created through soft and collinear emissions. Finally, in Sec. 5.4
we study a non-global observable at SL accuracy, which probes the description
of soft large-angle emissions in the shower.

For each observable, we compare the shower’s predictions to the known
resummation. Subleading colour corrections are accounted for using the NODS
scheme. Although it is possible to include spin correlations (see Appendix D),
this is not done for the tests described below as none of the observables is
sensitive to them at our targeted accuracy.

###  5.1 DGLAP evolution

Here we test the capability of our showers to reproduce the DGLAP
Gribov:1972ri ; Dokshitzer:1977sg ; Altarelli:1977zs  evolution of parton
distribution functions (PDFs). More specifically, we fix the starting point of
the shower to be the tree-level DIS process $q\,\gamma^{*}\to q$, where the
incoming quark $q$ is a down quark and has an energy fraction
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}$,
the photon has a space-like momentum equal to
$q_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}$
with
$Q^{2}=-q_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}^{2}$,
and we run our showers from the scale $Q$ until a given transverse-momentum
cutoff $p_{t,\rm{cut}}=Qe^{L}$. This cutoff is applied to the scale used as
the argument of the running coupling, and it corresponds to the ordering-scale
$v$ for Dipole-$k_{t}$, and to $\kappa_{\perp}$ defined in eq. (20) for the
PanScales showers. The flavour and energy fraction of the original quark
extracted from the proton can change throughout the shower evolution. When the
showering terminates the parton extracted from the proton will have a flavour
denoted by $i$, and an energy fraction $x\geq
x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}$.
The dominant term of the distribution over $i,x$ is SL ($(\alpha_{s}L)^{n}$
with $L=\ln p_{t,\rm{cut}}/Q$). The expected distribution will take the form

$\frac{1}{\sigma}\frac{d\sigma_{i}}{dx}=\frac{1}{f_{q}(x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS},Q^{2})}\int_{x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}}^{1}\frac{{\rm
d}z}{z}D_{qi}\left(z,\alpha_{s}L\right)f_{i}\left(\frac{x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}}{x},p_{t,{\rm
cut}}^{2}\right)\delta\left(\frac{x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}}{x}-z\right),$
(52)

where $f_{j}(y,\mu^{2})$ is the density of partons of flavour $j$, carrying
momentum fraction $y$ at a factorisation scale $\mu$. The (single-logarithmic)
DGLAP evolution operator $D_{ij}(z,\alpha_{s}L)$ is defined such that the
parton density functions satisfy

$\displaystyle f_{q}(x,\mu^{2})=\sum_{j}\int_{x}^{1}\frac{{\rm
d}z}{z}\,z\,D_{qj}(z,\alpha_{s}L)f_{j}\left(\frac{x}{z},p_{t,{\rm
cut}}^{2}\right)\,.$ (53)

Figure 6: The ratio of the DGLAP evolution produced by the
PanLocal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$) (upper panel) and
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$) (lower panel) showers versus the DGLAP evolution as calculated with
HOPPET. The results are shown as a function of the momentum fraction $x$
carried by the parton $i$ extracted from the proton. We fix the underlying
Born DIS process to have a $d$ quark in the initial-state with
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}=0.001$,
and fix $\alpha_{s}=5\cdot 10^{-6}$ with $\lambda=\alpha_{s}L=-0.5$. The three
columns show different extracted flavours $i$ that can lead to a $d$ quark
with energy fraction
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}$
entering the hard scattering process: we have a $d$ quark on the left panel, a
gluon in the middle panel, and any other flavour in the right panel.

In Fig. 6 we illustrate the results for the
PanLocal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$) and
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$) showers. The other showers,
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$) and Dipole-$k_{t}$, lead to analogous results. The DGLAP reference
prediction is obtained with HOPPET Salam:2008qg  at single logarithmic
accuracy, probing only those initial conditions that lead in the evolution to
$\delta_{{\tilde{\imath}}i}\delta(x-x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS})$
at the hard scale $\mu_{\text{{f}}}=Q$, with
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}=0.001$
and ${\tilde{\imath}}$ is a $d$ quark. The shower runs are performed by fixing
the underlying Born process to contain a $d$ quark in the initial-state with
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}=0.001$,
and the running of the strong coupling is performed at one loop. We show
results obtained by setting $\alpha_{s}(Q)=5\cdot 10^{-6}$ and a large value
of $L=-10^{5}$ ($\lambda=-0.5$), which sets any terms beyond SL accuracy to
zero. To further speed up the calculation we discard radiation with a momentum
fraction below some finite but small threshold $e^{-11}$, and only keep
radiation with an absolute rapidity larger than $18$ (i.e. very collinear to
the initial- or final-state quark). We have verified that these cuts do not
impact the results. Further details on the treatment of the PDFs may be found
in Appendix A of Ref. vanBeekveld:2022ukn . We see that agreement with the
HOPPET predictions is obtained for all showers to within the statistical
accuracy (below $0.1\%$ for the majority of the $x$ range).

###  5.2 Particle multiplicity

In this section, we examine the shower’s ability to reproduce the analytic
prediction of the average particle multiplicity given the shower’s transverse-
momentum cutoff. Although particle multiplicity is not an infrared-safe
quantity, the resummation structure of particle multiplicity, defined with an
infrared cutoff, can be related to that of subject multiplicity at NDL, which
is theoretically well-defined, and does not depend on the jet algorithm (as
long as it is infrared safe). The logarithmic accuracy of particle
multiplicity needs to be determined at the level of the distribution rather
than the logarithm of the distribution. For such non-exponentiating
observables, one writes

$\langle N_{\rm
NDL}(\alpha_{s},\xi)\rangle=h_{1}(\xi)+\sqrt{\alpha_{s}}h_{2}(\xi)\,,$ (54)

with $\xi=\alpha_{s}L^{2}$ and where $h_{1}$ collects double-logarithmic (DL)
$\mathcal{O}(\alpha_{s}^{n}L^{2n})$ terms, and $\sqrt{\alpha_{s}}h_{2}$ the
next-to-double logarithmic (NDL) $\mathcal{O}(\alpha_{s}^{n}L^{n})$ terms. At
the desired level of accuracy (NDL), the average multiplicity obtained by
counting the number of emissions generated by the shower with a transverse-
momentum cutoff $p_{t,\rm cut}=Qe^{L}$ (or equivalently, a shower that uses a
strong coupling equal to zero below a given value of $p_{t,\rm cut}$) is
equivalent to that obtained after using a well-defined jet algorithm. The
logarithm that is resummed then takes the form $L=\ln p_{t,\rm cut}/Q$.

For the process at hand, if we fix the Born flavour to be
${\tilde{\imath}}=q$, we have Catani:1991pm ; Catani:1993yx ; Medves:2022ccw

$\displaystyle h_{1}(\xi)$
$\displaystyle=\frac{2C_{F}}{C_{A}}\left[\cosh\nu-1\right]+N_{b}\,,$ (55a)
$\displaystyle h_{2}(\xi)$ $\displaystyle=\sqrt{\frac{1}{2\pi
C_{A}}}\sinh\nu\,\frac{Q^{2}{\partial}\ln
f_{q}(x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS},Q^{2})}{\partial
Q^{2}}+\frac{2C_{F}}{C_{A}}\Bigg{[}\frac{\beta_{0}}{2}\sqrt{\frac{\pi}{2C_{A}}}\left(\nu\cosh\nu+(\nu^{2}-1)\sinh\nu\right)$
$\displaystyle\hskip
28.45274pt+\sqrt{\frac{C_{A}}{2\pi}}B_{gq}\left(\frac{2C_{F}-C_{A}}{C_{A}}\nu\cosh\nu+\frac{5C_{A}-6C_{F}}{C_{A}}\sinh\nu+4\frac{C_{F}-C_{A}}{C_{A}}\nu\right)$
$\displaystyle\hskip
28.45274pt+\sqrt{\frac{C_{A}}{2\pi}}\left(2B_{qq}\sinh\nu+B_{gg}\left(\nu\cosh\nu-\sinh\nu\right)\right)\Bigg{]},$
(55b)

with $N_{b}=1$ the number of final-state partons at LO, and

$\displaystyle\nu=\sqrt{\frac{2C_{A}\xi}{\pi}}\,,\quad
B_{gg}=-\frac{11}{12}\,,\quad B_{qq}=-\frac{3}{4}\,,\quad
B_{gq}=\frac{2n_{F}}{3C_{A}}\,,\quad\beta_{0}=\frac{11C_{A}-2n_{F}}{12\pi}\,.$
(56)

As before, $f_{q}$ denotes the PDF of the incoming parton at LO, which carries
an energy fraction
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}$.

At DL accuracy, this test probes the soft-collinear nested structure of the
shower. To reproduce the average multiplicity at NDL, the parton shower must
also correctly incorporate hard-collinear corrections to the splitting
functions (corresponding to the second and third line of Eq. (55b)), the
running of the coupling constant (first line of Eq. (55b)), the DGLAP
evolution of the PDF (last line of Eq. (55b)), and the colour terms. The
treatment of these contributions is the same across all showers, hence we
expect to see agreement. The parton shower correctly reproduces the analytic
expectation at NDL if the ratio

$\frac{\langle N_{\rm shower}(\alpha_{s},\xi)\rangle-\langle N_{\rm
NDL}(\alpha_{s},\xi)\rangle}{\langle N_{\rm
NDL}(\alpha_{s},\xi)\rangle-\langle N_{\rm DL}(\alpha_{s},\xi)\rangle}\,,$
(57)

vanishes in the $\alpha_{s}\to 0$ limit, where $\xi$ is kept fixed and
$\langle N_{\rm DL}(\alpha_{s},\xi)\rangle$ corresponds to Eq. (54) with
$h_{2}$ set to $0$.

Figure 7: Extrapolation of $\frac{\langle N_{\rm
shower}(\alpha_{s},\xi)\rangle-\langle N_{\rm
NDL}(\alpha_{s},\xi)\rangle}{\langle N_{\rm
NDL}(\alpha_{s},\xi)\rangle-\langle N_{\rm DL}(\alpha_{s},\xi)\rangle}$ for
$\alpha_{s}\to 0$ at fixed value of $\xi=\alpha_{s}L^{2}=5$ for several
showers considering the DIS process $d\gamma^{*}\to d$ with $x_{\rm
DIS}=0.001$.

To extract the $\alpha_{s}\to 0$ limit, we run the PanScales showers for the
DIS process $q\gamma^{*}\to q$, where we set $q=d$,
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}=0.001$,
$\xi=5$ and sample $\alpha_{s}=\left\\{8\cdot
10^{-5},\,0.00128,0.003156\right\\}$, performing a quadratic polynomial
extrapolation to get the $\alpha_{s}\to 0$ result. Systematic uncertainties
are estimated performing an alternative extrapolation with
$\alpha_{s}=\,0.00512$ instead of $\alpha_{s}=0.003156$, and added in
quadrature to the statistical uncertainty. Like in the previous section, the
running of the coupling constant is performed at one loop, as the $2$-loop
running only enters at NNDL accuracy. Subleading colour corrections are
included using the NODS scheme. The obtained result is shown in Fig. 7. We
notice that all the showers are consistent with the NDL expectation, with an
uncertainty well below 0.1%.

###  5.3 Continuously-global observables

This section details tests on a range of global Banfi:2004yd  observables. The
cumulative distribution of such observables, i.e. the probability that the
observable $O$ takes a value smaller than $e^{L}$, can be written as a
function of $\lambda=\alpha_{s}L$, and takes the form

$\displaystyle\Sigma(O<e^{L})=H(\alpha_{s})\exp\left[-L\,g_{1}(\lambda)+g_{2}(\lambda)+\dots\right]+\dots\,,$
(58)

where $g_{1}$ contains the LL, and $g_{2}$ the NLL contribution. The function
$H(\alpha_{s})$ is the hard function, which can be set equal to $1$ at our
targeted NLL accuracy. To test the accuracy of the shower we examine

$\displaystyle\lim_{\alpha_{s}\to 0}\frac{\Sigma_{\rm PS}(\lambda)-\Sigma_{\rm
NLL}(\lambda)}{\Sigma_{\rm NLL}(\lambda)}\,,$ (59)

at a fixed value of $\lambda$, which should tend to $0$ if the shower is NLL
accurate. To extract the $\alpha_{s}\to 0$ limit, we need to run the shower at
very small $\alpha_{s}$ values. We are able to do so thanks to the numerical
techniques developed in Refs. Dasgupta:2020fwr ; Hamilton:2020rcu ;
vanBeekveld:2022ukn . These techniques assume that in the soft-collinear
limit, the observable scales as

$\displaystyle
O\sim\frac{p_{T}}{Q}e^{-\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}|\eta|}\,,$ (60)

where
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}\geq 0$ is a constant number. This implies that we can only consider
continuously-global observables Banfi:2004yd  in our tests.

Keeping this technical limitation in mind, to test our showers, we introduce
three sets of observables, parameterised by
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}\in[0,0.5,1]$:

$\displaystyle
S_{p,\,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$ $\displaystyle=\sum_{j\in{\rm
partons}}\frac{k_{t,j}e^{-\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}|y_{j}|}}{Q}\,,$ (61a) $\displaystyle
S_{j,\,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$ $\displaystyle=\sum_{j\in{\rm
jets}}\frac{k_{t,j}e^{-\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}|y_{j}|}}{Q}\,,$ (61b) $\displaystyle
M_{j,\,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$ $\displaystyle=\max_{j\in{\rm
jets}}\frac{k_{t,j}e^{-\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}|y_{j}|}}{Q}\,.$ (61c)

For the particle observable
$S_{p,\,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$, the sum over $j$ runs over all partons, and $k_{t,j}$, $y_{j}$
correspond to the transverse momentum and the rapidity of the parton in the
Breit frame. Jets are defined with the algorithm detailed in Appendix B. For
the two jet observables, the sum over $j$ runs over all jets found inside the
collection of beam jets and the final-state macro-jet. The $k_{t,j}$ and
$y_{j}$ are their primary Lund-plane coordinates. The analytic expectations
for the observables quoted in Eq. (61) are collected in Appendix C. Note that
$S_{p,\,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$ will only be computed for
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}>0$, i.e. when its NLL prediction corresponds to the one of a standard
additive observable and is equivalent to
$S_{j,\,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$. For
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}>0$ the contribution of the original hard final-state leg to
$S_{p,\,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$ is subleading and can be neglected, thus yielding the same result as
$S_{j,\,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$. This is no longer the case for
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=0$, as recoil effects contribute at NLL. We stress that these observables
are not directly measurable, because they are built from emissions in both the
remnant and the current hemisphere, but they have the property of having a
remarkably simple resummation structure at NLL. For this reason, they can be
easily used to test the logarithmic accuracy of showers, but they can also be
employed as resolution variables to build slicing methods Catani:2007vq ;
Stewart:2010tn  or NNLO+PS matching prescriptions Hamilton:2012rf ;
Alioli:2013hqa ; Monni:2019whf .

We further test our showers with observables that are phenomenologically
accessible. To remove contamination from the fragmenting beams in the remnant
hemisphere, allowing for a better experimental measurement, event shapes in
DIS are often defined in the current hemisphere $\mathcal{H}_{c}$. Examples of
such observables, which are continuously-global and satisfy the recursive
infrared-safety requirement of Ref. Banfi:2004yd , are

$\displaystyle B_{zE}=$
$\displaystyle\frac{\sum_{i\in\mathcal{H}_{c}}|\vec{p}_{\perp,i}|}{2\sum_{i\in\mathcal{H}_{c}}|\vec{p}_{i}|}\,,$
(62a) $\displaystyle B_{zQ}=$
$\displaystyle\frac{\sum_{i\in\mathcal{H}_{c}}|\vec{p}_{\perp,i}|}{Q}\,,$
(62b) $\displaystyle\tau_{zQ}=$ $\displaystyle
1-\frac{2\sum_{i\in\mathcal{H}_{c}}|p_{z,i}|}{Q},$ (62c)

where all the quantities are defined in the Breit frame, and the suffix $z$
implies they are measured with respect to the photon axis. The analytic
predictions for the two definitions of broadening, $B_{zE}$ and $B_{zQ}$, can
be found in Ref. Dasgupta:2001eq . They are examples of
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=0$ observables, and their NLL prediction is identical. For $\tau_{zQ}$,
which is a
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=1$ observable, these predictions were computed in Ref. Antonelli:1999kx .
These results are also summarised in Appendix C.

(a)

(b)

Figure 8: Cumulative distribution for the broadening normalised with respect
to the energy in the current hemisphere, $B_{zE}$, for the process
$q\gamma^{*}\to q$ with
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}=0.2$.
In the left panel we show
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$) results, obtained with finite values of $\alpha_{s}$, as well as the
$\alpha_{s}\to 0$ extraction. In the right panel, we illustrate the
$\alpha_{s}\to 0$ extraction for all the showers: the ratio with the analytic
prediction must be 1 if the shower is NLL. Note that the values for the
PanLocal shower have been shifted to the left to improve readability, i.e.
their probed $\lambda$ values coincide with those for
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$).

In Fig. 8(a) we show the ratio of the
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$) prediction to the NLL result for the cumulative distribution of the
broadening normalised with respect the energy in the current hemisphere,
$B_{zE}$, for increasingly smaller values of $\alpha_{s}$,

$\displaystyle\alpha_{s}=\\{0.0015625,0.003125,0.00625\\}$ (63)

and $\lambda\geq-0.5$. The coloured band represents only the statistical
uncertainty. In black we illustrate the result of the $\alpha_{s}\to 0$
extrapolation, which has been performed with a quadratic interpolation. In
Fig. 8(b) we summarise the $\alpha_{s}\to 0$ extractions for all the showers.
From this figure it is clear that all the new PanScales showers for DIS
processes reproduce the analytic expectation, while we observe deviations for
Dipole-$k_{t}$ reaching up to 15% for $\lambda=-0.5$. We obtain identical
results for $B_{zQ}$, not shown here.

Figure 9: The $\alpha_{s}\to 0$ extrapolation of the deviation of the shower
result from the NLL expectation for the process $q\gamma^{*}\to q$ with
$x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS}=0.2$
and $\lambda=-0.5$ for several continuously-global observables. The red colour
signals a deviation from the analytic expectation larger than two standard
deviations, while the green colour signals an agreement. When numerical
agreement between the shower and the analytic resummation is found, but fixed-
order issues are known, we use the amber colour.

The shower’s expectation for the other observables are summarised in Fig. 9,
where we show the ratio with the NLL result for the cumulative distribution
$\Sigma(O<e^{L})$ in the $\alpha_{s}\to 0$ limit with $\lambda=-0.5$. Here,
the central value is again obtained using the $\alpha_{s}$ values in eq. (63),
but the uncertainty is given by summing in quadrature the statistical
uncertainty and the difference between the central value and the one obtained
by performing the extrapolation on the set of $\alpha_{s}$ values where
$\alpha_{s}=0.003125$ is replaced with $\alpha_{s}=0.0125$.

The PanScales showers agree with the analytic expectations for all the
observables. Dipole-$k_{t}$ leads to manifestly wrong all-order results at NLL
for
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=0$ observables, with a 4% deviation for $S_{j,0}$ and a 12% deviation for
$M_{j,0}$. Like in Refs. Dasgupta:2020fwr ; vanBeekveld:2022ukn , we observe
that despite the fixed-order issue that we highlighted in Sec. 4.1,
Dipole-$k_{t}$ seems to reproduce the correct analytic expectations for
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}>0$. The NLL-violating terms manifest themselves as super-leading
logarithms that violate the exponentiation, as observed in Ref.
Dasgupta:2020fwr , which however resum to 0 in the all-orders limit for
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}>0$.

###  5.4 Non-global logarithms

Many observables are sensitive to radiation in a restricted portion of the
Lund plane. In the context of DIS, these observables are, for example, the
current jet mass, the $C$-parameter, the thrust with respect to the current-
hemisphere thrust-axis Dasgupta:2002dc  and $Q_{t}$, i.e. the transverse
momentum of the system comprising the partons in the current hemisphere
Dasgupta:2006ru . For what concerns Higgs production in VBF, isolation
criteria can be used to reduce this production mode from gluon fusion
Buckley:2021gfw . The resummation for these observables naturally involves
non-global logarithms (NGLs) Dasgupta:2001sh ; Dasgupta:2002bw , which can be
correctly reproduced at leading-colour single-logarithmic (SL) accuracy only
by dipole showers Banfi:2006gy .

To assess the ability of our showers to reproduce such NGLs, we consider the
scalar sum of the transverse momenta of the partons in a rapidity slice

$S_{\Delta}^{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}slice}=\frac{\sum_{i}k_{t,i}\Theta(|\eta_{i}|<\Delta)}{Q}\equiv\frac{k_{t,\text{slice}}}{Q}\,,$
(64)

where the transverse momentum $k_{t,i}$ and the rapidity $\eta_{i}$ of the
partons are defined in the Breit frame. NGLs for this observable are single-
logarithmic terms of the form $\lambda^{n}=\alpha_{s}^{n}L^{n}$, created by
soft large-angle emissions near the edge of the slice, i.e. near
$y=\pm\Delta$.

(a)

(b)

Figure 10:  (a) Cumulative distribution for the transverse momentum in a
rapidity slice of $|y|<1$ as a function of
$\lambda=\alpha_{s}\ln\frac{k_{t,\text{slice}}}{Q}$ for the PanScales showers
and Dipole-$k_{t}$. The top panel shows the expected (black dashed) and the
showers (solid) results, while the bottom panels show the ratio between the
shower and the analytic prediction for each of the showers. (b) Relative
difference between the shower $\Sigma_{\rm PS}$ and the expected single-
logarithmic result ($\Sigma_{\rm SL}$) for a fixed value of $\lambda=-0.5$ for
all the PanScales showers. Colour coding is like in Fig. 9.

Fig. 10(a) shows the comparison between the PanScales shower predictions with
the expected results for $-0.5\leq\lambda<0$ and $|y|<1$, while in Fig. 10(b)
we show the results for $\lambda=-0.5$. We generate our reference calculation
in the large-$N_{\text{{c}}}$ limit (with $C_{F}=C_{A}/2=3/2$) from the code
developed for Ref. Caletti:2021oor , which uses the strategy of Ref.
Dasgupta:2001sh . The shower predictions are obtained running with
$\alpha_{s}=10^{-9}$, so that NSLs are numerically negligible. To reduce the
parton multiplicity, without affecting the observable under consideration, we
impose a rapidity and a soft-emission cut, vetoing radiation with $|\eta|>13$
and $\ln k_{t}<\ln k_{t,\text{slice}}-18$. Furthermore, like in Secs. 5.1 and
5.2, the running of the coupling constant is performed at one loop, with
$K=0$, as these effects only enter at NSL/NNDL.

In all cases, we notice an excellent agreement between the PanScales showers
and the correct SL distribution. This is true also for Dipole-$k_{t}$, despite
fixed-order issues, exactly like in the case of continuously-global event
shapes with
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}>0$.

##  6 Phenomenological results Higgs production in VBF

We now move to the PanScales showers for VBF, presenting some exploratory
phenomenological results, and explaining more details on the implementation of
such process in our framework. This channel provides a clean experimental
signature and is therefore an ideal environment to study the Higgs
boson.999The VBF channel is also used to search for di-Higgs production at the
LHC ATLAS:2020jgy .  As already explained in Sec. 3.3, at NLO the VBF channel
can be seen as two independent DIS processes for each hadronic sector. The two
jets that are formed after the initial-state quark emits a vector boson are
typically produced with a large absolute rapidity. Colour coherence then
results in little jet activity in the the central rapidity region; radiation
will be primarily concentrated around the two hard jets at (opposite) large
rapidities.

For the phenomenological studies we produce the Higgs boson via the $ZZ$
channel in a VBF topology with a centre-of-mass energy of 13.6 TeV. We run all
showers using the NODS colour scheme, even though the CFFE scheme is the one
adopted by standard dipole showers, to more faithfully gauge the kinematic
differences between the LL and NLL showers. The hard process is obtained from
Pythia8.3 Bierlich:2022pfr  at LO accuracy, using the default values for the
electroweak parameters and the Higgs mass. We use the NNPDF 4.0 LO PDF set
with perturbative charm content NNPDF:2021njg  (LHAPDF label 332500
Buckley:2014ana ), corresponding to $\alpha_{s}(m_{Z})=0.118$. The default
factorisation and renormalisation scale used in Pythia8.3 to generate the hard
process is

$\mu_{F,h}=\mu_{R,h}\equiv\mu_{h}=\sqrt[1/3]{m_{T,H}m_{T,V_{1}}m_{T,V_{2}}},\quad
m_{T,i}=\sqrt{p_{\perp,i}^{2}+m_{i}^{2}}\,,$ (65)

where $V_{1,2}$ denote the vector-bosons exchanged in the $t$-channel
propagators for each of the two hadronic sectors, and $m_{H,V_{1},V_{2}}$ the
Higgs/vector-boson masses. However, as discussed in Sec. 3.3, another choice
would be to use the virtualities of the exchanged boson as two independent
scales for the two hadronic sectors. To this end, we apply the reweighting
factor

$\displaystyle
w=\frac{f_{i_{1}}\left(x_{1},\mu_{1}\right)\,f_{i_{2}}\left(x_{2},\mu_{2}\right)}{f_{i_{1}}\left(x_{1},\mu_{h}\right)\,f_{i_{2}}\left(x_{2},\mu_{h}\right)}\,,$
(66)

where $f_{i_{1,2}}$ is the PDF of the incoming quark $i_{1,2}$, which carries
an energy fraction $x_{1,2}$, and $\mu_{1,2}=\sqrt{Q_{1,2}^{2}}$. Secondly, we
start the parton shower at a distinct scale for each of the hadronic sectors,
i.e.

$v^{2}_{\max,i}=Q_{i}^{2}\frac{1-x_{i}}{x_{i}}\,.$ (67)

In practice, we use the maximum of the two values as common starting scale,
and we perform a veto to ensure that $v<v_{\max,i}$ in each hadronic section.

We estimate the uncertainty stemming from renormalisation scale variations
using a modified scheme for $\alpha_{s}$ Mrenna:2016sih ; vanBeekveld:2022ukn
, i.e.

$\displaystyle\alpha_{s}\left(\mu_{R}^{2}\right)\left(1+\frac{\alpha_{s}\left(\mu_{R}^{2}\right)}{2\pi}K+2\alpha_{s}\left(\mu_{R}^{2}\right)b_{0}(1-z)\ln
x_{R}\right)\,,\quad\mu_{R}=x_{R}\mu_{R,0}\,,$ (68)

with $\mu_{R,0}=\rho v{\rm
e}^{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}|\bar{\eta}_{Q_{i}}|}$ the central scale. The factor $z$ is the fraction of
the emitter-momentum carried away by the radiation. With this, $1-z$ ensures
that scale compensation at NLL is present for soft emissions, but not for hard
emissions. We omit this term for the LL-accurate shower Dipole-$k_{t}$, but
include the CMW factor $K$. Renormalisation-scale variations are probed taking
$x_{R}\in\\{0.5,1,2\\}$, and the infrared cutoff of the shower is implemented
such that $\alpha_{s}(\mu_{R}^{2})=0$ for $\mu_{R}<x_{R}\times 0.5$ GeV.
Factorisation-scale uncertainties are probed independently from the
renormalisation scale, and are assessed using

$\displaystyle\ln\mu_{F}\equiv\ln\left(x_{F}\mu_{F,0}\right)=\ln
Q_{i}+\frac{1}{1+\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}}\ln\frac{v_{i}}{Q_{i}}+\ln x_{F}\,,$ (69)

with $x_{F}\in\\{0.5,1,2\\}$. In our results we show the $7$-point scale
variation, obtained by taking
$(x_{R},x_{F})\in\\{(1,1),(0.5,1),(1,0.5),(0.5,0.5),(2,1),(1,2),(2,2)\\}$.

Figure 11: Distributions normalised to the total LO cross section $\sigma_{\rm
LO}$ before cuts for $m_{j_{1}j_{2}}$ (top) and $\Delta\eta_{j_{1}j_{2}}$
(bottom) after applying the minimal selection cuts (left) and the VBF cuts
(right). We show
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$), black dashed,
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$), blue dash-dotted,
PanLocal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$), green solid and Dipole-$k_{t}$, orange, dotted. The label $w$
indicates the total area under the curve. The vertical lines on the central
value of each bin indicate the statistical uncertainty, which is mostly
negligible and therefore not visible. The band indicates the
renormalisation/factorisation scale uncertainty, obtained as explained after
eq. (69). The bottom panels show the ratio with respect to the
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0$) result.

We use the anti-$k_{T}$ algorithm Cacciari:2008gp  with $R=0.4$, implemented
in FastJet Cacciari:2011ma  to cluster jets with the definition $p_{T,j}>25$
GeV, $|\eta_{j}|<4.5$, and consider two setups:

  * •

“minimal cuts”, where we require the presence of two resolved jets;

  * •

“VBF cuts”, where we require that the two leading jets (i.e. those two with
the largest transverse momenta) are separated by a rapidity
$\Delta\eta_{j_{1}j_{2}}>4.5$, have a dijet invariant mass of at least
$m_{j_{1}j_{2}}>600$ GeV and lie in opposite hemispheres
($\eta_{j_{1}}\eta_{j_{2}}<0$).

The results shown below do not include the simulation of beam remnants,
hadronisation or multi-parton interaction. We have the option to incorporate
them into our framework trivially using Pythia8.3, but we have made the
decision not to do so. This choice enables us to provide a clearer comparison
between the NLL and LL accurate showers.

Figure 12: As Fig. 11 but for $\eta_{j_{3}}$.

We first consider inclusive observables, i.e. those observables that are non-
vanishing at LO. The value of these observables is primarily set by the hard
scattering process, and the shower should impact these observables only
marginally. In Fig. 11 we show two such observables: the invariant mass of the
two leading jets $m_{j_{1}j_{2}}$ and the rapidity seperation between them
$\Delta\eta_{j_{1}j_{2}}$. As expected, one observes that uncertainties from
renormalisation/factorisation scale variations in the shower are small, as
these observables are largely insensitive to the additional radiation produced
by the shower. One exception is given by the regions where
$\Delta\eta_{j_{1}j_{2}}$ or $m_{j_{1}j_{2}}$ are small (after applying the
minimal selection cuts). These regions are contaminated by events where one of
the two tagged jets actually originates from the shower hardest emission.
Indeed, because of this effect, the scale uncertainties here are larger than
in the bulk of the distribution. With the exception of these regions, which
are excluded after applying the VBF cuts, it can be observed that differences
between the NLL showers are of the order of $\mathcal{O}(5-10\%)$, i.e.
commensurate with known NLO corrections, and the LL shower lies in between the
prediction of the NLL showers.

Figure 13: Normalised distributions for $p_{T,j_{3}}$ after applying the VBF
cuts. The right panel shows a magnification of the (relatively) small
transverse-momentum region. The lines, bands and ratio plots are the same as
in Fig. 11.

We now turn to observables that can only be defined in the presence of a third
jet, and therefore are a clean probe of the shower’s behaviour. In particular,
we consider the pseudo-rapidity of the third jet, $\eta_{j_{3}}$ (Fig. 12),
and the transverse momentum of the third jet $p_{T,j_{3}}$ (Fig. 13, after
applying the VBF cuts). Colour coherence predicts a suppression of radiation
in the central rapidity region. All the showers considered here indeed show
this behaviour, as can be observed in Fig. 12. It can be observed in the
high-$p_{T,j_{3}}$ region of Fig. 13 (left), and in Fig. 12 (right, after we
apply VBF cuts) that scale variation indeed cannot be used to accurately
capture the differences between different showers. We can see that
Dipole-$k_{t}$ lies outside the spread of NLL showers for these exclusive
observables for most of the shown phase space. Differences between
Dipole-$k_{t}$ and the NLL showers are most pronounced in the bulk of the
distributions: for $\eta_{j_{3}}$ the central values differ by
$\mathcal{O}(10-15\%)$ for $|\eta_{j_{3}}|$ around $3$, and by
$\mathcal{O}(10\%)$ for $p_{T,j_{3}}$ in the relatively small transverse-
momentum region. These normalisation differences are due to the fact that
Dipole-$k_{t}$ shows a higher rate of $3$-jet events than the NLL showers.
Conversely, for larger values of $p_{T,j_{3}}$, differences between the NLL
showers also start to grow. However, since we do not include matching, this
hard-emission region is not controlled in our showers, and we leave a further
detailed investigation of the differences between the LL and NLL showers for
future work.

##  7 Conclusions

In this work, we have introduced new NLL-accurate dipole showers for processes
involving the exchange of a colour-singlet in the $t$-channel, such as DIS,
VBF and VBS. The latter two processes are handled following a factorised
approach, i.e. neglecting non-factorisable corrections between the two
hadronic sectors. The main novelty of these showers, with respect to the
PanScales showers for hadron collisions introduced in Refs.
vanBeekveld:2022zhl ; vanBeekveld:2022ukn , is that the transverse-momentum
recoil due to initial-state radiation is smoothly redistributed primarily to
partons in the current hemisphere (i.e. anti-parallel to the direction of the
incoming proton in the Breit frame). This feature ensures that partons in the
remnant hemisphere remain mostly unaffected, which is required from colour
coherence. Furthermore, compared to standard showers for DIS and VBF/VBS, our
showers differ in the choice of the dipole partitioning: this is not done in
the dipole frame, but in the Breit frame. This, combined with a global recoil
scheme or with a careful choice of the ordering scale, prevents soft gluons
from taking unphysical recoil.

We have carried out a number of fixed-order tests, focusing on DIS, such as
analysing the phase-space contours for two emissions with commensurate
softness (Section 4.1), and colour/spin matrix-element comparisons (Section
4.2 and Appendix D), related to the PanScales conditions needed to achieve NLL
accuracy. All-order validations of our new showers for DIS have also been
carried out for a variety of observables. These include tests of the DGLAP
evolution (Section 5.1), jet-multiplicity (Section 5.2), DIS continuously-
global event-shapes (Section 5.3), and the scalar sum of transverse momenta in
a fixed rapidity slice (Section 5.4). In our comparisons, we introduced new
continously-global event shapes, which rely on the use of Lund Plane
coordinates and a jet algorithm specific for DIS. All these tests were carried
out including subleading colour corrections, except for those of the non-
global observable. The PanLocal shower with
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0.5$ and PanGlobal showers with
$\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
PS}=0,0.5$ sucessfully pass the fixed- and all-order NLL accuracy tests. We
have compared these showers to a ‘standard’ transverse-momentum ordered
shower, Dipole-$k_{t}$, which has fully local transverse-momentum recoil.
Differences between the LL-accurate Dipole-$k_{t}$ and the new NLL-accurate
showers can grow up to $15\%$ for phenomenologically relevant continuously-
global event-shape observables such as the current-hemisphere
broadening.101010The NLL test for the Dipole-$k_{t}$ shower has been performed
implementing subleading colour corrections with the NODS colour scheme,
despite the fact that the colour-factor-from-emitter scheme is the standard
choice for such showers.  The SL accuracy of our showers established for DIS
proceses, automatically propagates to VBF/VBS processes in the factorised
approach. This is the first time such accuracy is achieved for VBF/VBS global
and non-global observables. The neglected non-factorisable contributions are
SL subleading colour corrections, and are typically further surpressed after
applying VBF cuts.

In Section 6 we present an exploritory phenomenological application of our
newly developed showers for Higgs production in VBF at $\sqrt{s}=13.6$ TeV.
For each of the two hadronic sectors individually, we choose separate values
for the shower starting, renormalisation, and factorisation scales. While we
examined the impact of variations in the renormalization and factorization
scales, it is important to note that such variations typically do not capture
all sources of uncertainty. Indeed, the fact that we have developed not one,
but several showers, is important for a realistic estimate of shower
uncertainties. We find that for inclusive observables the differences among
the showers are of the order of 5-10%, and the LL shower is contained in the
spread of our NLL ones. This spread, which is not captured by scale-variation
uncertainties, demonstrates the importance of considering multiple showers to
obtain a realistic shower uncertainty. For exclusive observables, like the
rapidity or transverse momentum of the third jet, we find that the
Dipole-$k_{t}$ predictions typically overshoot the NLL showers by 5-15%.
Conversely, our NLL showers are in good agreement between themselves in the
bulk of the distributions. Next steps involve matching the NLL-accurate
showers to the NLO fixed-order results, where care needs to be taken not to
compromise the NLL accuracy Hamilton:2023dwb , as well as including heavy-
quark mass effects to handle $t$-channel single-top production.

## Acknowledgements

We are grateful to our PanScales collaborators (Mrinal Dasgupta, Frédéric
Dreyer, Basem El-Menoufi, Keith Hamilton, Jack Helliwell, Alexander Karlberg,
Rok Medves, Pier Monni, Gavin Salam, Ludovic Scyboz, Alba Soto-Ontoso, Gregory
Soyez, Rob Verheyen, and Scarlett Woolnough), for their work on the code,
comments on the manuscript, and the underlying philosophy of the approach. In
particular, we want to thank Pier Monni and Mrinal Dasgupta for having shared
with us their knowledge on resummation for DIS event shapes, Gavin Salam for
help with the implementation of jet algorithm for DIS, Gregory Soyez for
frequent discussions on the technical details of the NLL tests, and Alexander
Karlberg for pointing out relevant VBF literature. This work was supported by
a Royal Society Research Professorship (RP$\backslash$R1$\backslash$180112)
(MvB), by the European Research Council under the European Union’s Horizon
2020 research and innovation programme (grant agreement No. 788223, PanScales)
(MvB and SFR), and by the Science and Technology Facilities Council under
ST/T000864/1 (MvB).

##  Appendix A Description of the momentum-conservation restoring boost

This appendix details the derivation of the boost $\Lambda^{\mu\nu}$,
introduced in Sec. 3.2. For PanGlobal this boost acts on the collection of
final-state partons
$\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}$
to restore momentum conservation. For the PanLocal shower it in addition also
acts on the incoming parton $\bar{p}_{1}^{\mu}$ to restore momentum
conservation after the collection of partons are rotated so that
$\bar{p}_{1}^{\mu}$ is aligned with the direction of $n_{1}^{\mu}$. We demand
that the boost preserves the invariant mass of
$\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}$,
which implies

$\displaystyle\Lambda^{\mu\nu}\left(\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X,\nu}\right)\equiv
p_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{\mu}=n_{2}^{\mu}+\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}{Q^{2}}n_{2}^{\mu}\,.$
(70)

Using a Sudakov decomposition of
$\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}$
along the directions of $n_{1}^{\mu}$ and $n_{2}^{\mu}$, eq. (70) can be
written as

$\Lambda^{\mu\nu}\left(\beta
n_{2,\nu}+\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}+p_{t}^{2}}{\beta
Q^{2}}n_{1,\nu}+p_{\perp,\nu}\right)=n_{2}^{\mu}+\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}{Q^{2}}n_{1}^{\mu}\,,$
(71)

with $p_{\perp}^{2}=-p_{t}^{2}<0$, $n_{1}^{2}=n_{2}^{2}=0$, $2n_{1}\cdot
n_{2}=Q^{2}$ and $p_{\perp}\cdot n_{1,2}=0$. We aim to design a boost that
acts as a rescaling for momentum components along $n_{1}^{\mu}$, and assigns
all the transverse momentum recoil to the $n_{2}^{\mu}$ direction.111111The
construction of this boost shares similarities with the boost used in Deductor
for $pp\to Z/h$ collisions Nagy:2009vg . Inverting eq. (71), we have

$\left[\Lambda^{-1}\right]^{\mu\nu}\left(n_{2,\nu}+\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}{Q^{2}}n_{1,\nu}\right)=\beta
n_{2}^{\mu}+\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}+p_{t}^{2}}{\beta
Q^{2}}n_{1}^{\mu}+p_{\perp}^{\mu}\,.$ (72)

We now may examine the action of the inverse boost on the $n_{1}^{\mu}$ and
$n_{2}^{\mu}$ components individually. The direction $n_{2}^{\mu}$ absorbs all
the transverse-momentum recoil, but also must stay massless, meaning it also
must absorb a component in the $n_{1}^{\mu}$ direction. To satisfy these
constraints we infer that

$\left[\Lambda^{-1}\right]^{\mu\nu}n_{2,\nu}\equiv\beta
n_{2}^{\mu}+\frac{p_{t}^{2}}{\beta Q^{2}}n_{1}^{\mu}+p_{\perp}^{\mu}\,.$ (73)

To satisfy eq. (72) we then have to require

$\left[\Lambda^{-1}\right]^{\mu\nu}n_{1,\nu}\equiv\frac{Q^{2}}{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}}\left(\beta
n_{2}^{\mu}+\frac{\bar{p}_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}X}^{2}+p_{t}^{2}}{\beta
Q^{2}}n_{1}^{\mu}+p_{\perp}^{\mu}-\left[\Lambda^{-1}\right]^{\mu\nu}n_{2,\nu}\right)=\frac{1}{\beta}n_{1}^{\mu}\,.$
(74)

We now have established the action of the inverse boost on the $n_{1}^{\mu}$
and $n_{2}^{\mu}$ directions. What remains is to determine what happens when
the inverse boost acts on a generic transverse-momentum component
$q_{\perp}^{\mu}$ (with $q_{\perp}\cdot n_{1,2}=0$). We require that the
action of the inverse boost on such a perpendicular component is only allowed
to bring in a component in the direction of $n_{1}^{\mu}$. Furthermore, the
perpendicular component does not get rescaled, so that its norm is preserved:

$\left[\Lambda^{-1}\right]^{\mu\nu}q_{\perp,\nu}=q_{\perp}^{\mu}+An_{1}^{\mu}\,.$
(75)

To find the value of $A$, we require that the invariant mass of generic vector
$q^{\mu}=a_{q}n_{1}^{\mu}+b_{q}n_{2}^{\mu}+q_{\perp}^{\mu}$ is preserved by
the boost. This leads to

$A=-\frac{2q_{\perp}\cdot p_{\perp}}{\beta Q^{2}}\,.$ (76)

Note that we cannot introduce a component along $n_{2}^{\mu}$ in the action of
the inverse boost on a transverse component. To see this, consider the action
of the inverse boost on a generic transverse component $q_{\perp}^{\mu}$ where
we now assign a $n_{2}^{\mu}$ component,

$\displaystyle\left[\Lambda^{-1}\right]^{\mu\nu}q_{\perp,\nu}=Bn_{2}^{\mu}+q_{\perp}^{\mu}\,.$
(77)

With this definition, the action of the inverse boost on a generic four-
momentum $q^{\mu}$ would become

$\displaystyle\left[\Lambda^{-1}\right]^{\mu\nu}q_{\nu}=\left(\frac{a_{q}}{\beta}+\frac{b_{q}\,p_{t}^{2}}{\beta
Q^{2}}\right)n_{1}^{\mu}+\left(b_{q}\beta+B\right)n_{2}^{\mu}+b_{q}p_{\perp}^{\mu}+q_{\perp}^{\mu}\,.$
(78)

Requiring that the boost leaves the invariant mass of $q^{\mu}$ unchanged
gives us

$\displaystyle B=\frac{2p_{\perp}\cdot
q_{\perp}}{a_{q}Q^{2}+b_{q}p_{t}^{2}}b_{q}\beta\,.$ (79)

This is not linear in $q^{\mu}$, which means that if we would consider
$l^{\mu}=q^{\mu}+q^{\prime\mu}$, we would have
$\left[\Lambda^{-1}\right]^{\mu\nu}l_{\nu}\neq\left[\Lambda^{-1}\right]^{\mu\nu}q_{\nu}+\left[\Lambda^{-1}\right]^{\mu\nu}q^{\prime}_{\nu}$.
Therefore, $B$ needs to be set to zero.

The action of $\left[\Lambda^{-1}\right]^{\mu\nu}$ is now fully specified, and
we can invert it to find $\Lambda^{\mu\nu}$. We then obtain

$\displaystyle\Lambda^{\mu\nu}n_{1,\nu}$
$\displaystyle={\beta}\,n_{1}^{\mu}\,,$ (80a)
$\displaystyle\Lambda^{\mu\nu}q_{\perp,\nu}$
$\displaystyle=q_{\perp}^{\mu}+\frac{2q_{\perp}\cdot
p_{\perp}}{Q^{2}}n_{1}^{\mu}\,,$ (80b)
$\displaystyle\Lambda^{\mu\nu}n_{2,\nu}$
$\displaystyle=\frac{1}{\beta}\left(n_{2}^{\mu}+\frac{p_{t}^{2}}{Q^{2}}n_{1}^{\mu}-p_{\perp}^{\mu}\right),$
(80c)

so that the final form of the boost reads

$\displaystyle\Lambda^{\mu\nu}=$ $\displaystyle
g^{\mu\nu}+\frac{2n_{1}^{\mu}}{Q^{2}}\left[(\beta-1)n_{2}^{\nu}+\frac{p_{t}^{2}}{\beta
Q^{2}}n_{1}^{\nu}+p_{\perp}^{\nu}\right]+\frac{2n_{2}^{\mu}n_{1}^{\nu}}{Q^{2}}\frac{1-\beta}{\beta}-\frac{2p_{\perp}^{\mu}n_{1}^{\nu}}{\beta
Q^{2}}\,.$ (81)

##  Appendix B Cambridge/Aachen algorithm for DIS and Lund variables

In this section, we describe the exclusive Cambridge/Aachen algorithm for DIS.
The $e^{+}e^{-}$ variant was introduced in Ref. Dokshitzer:1997in , while a
DIS variant was proposed in Ref. Wobisch:1998wt . This algorithm is similar to
the $k_{\perp}$ algorithm for DIS introduced of Ref. Dokshitzer:1997in , but
using an angular distance to determine the cluster sequence.

Conversely to the original proposal, here we do not employ a resolution
variable (or equivalently, we set $y_{\rm cut}=0$). Instead, the clustering
procedure we use is

  * •

For every final-state parton $i$, define the angular distance with respect to
the beam $B$

$d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}iB}=1-\cos\theta_{i}\,,$
(82)

where $\theta_{i}$ is the angle between $i$ and $B$ in the Breit frame. In
addition, define the angular distance between each pair of partons $i$ and
$j$, again in the Breit frame,

$d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}ij}=1-\cos\theta_{ij}\,.$
(83)

  * •

Find the smallest distance among
$\left\\{d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}iB},d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}ij}\right\\}$.
If it is of the type
$d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}iB}$,
we remove $i$ from the list of partons, and make it a candidate jet. If
instead it is of the type
$d_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}ij}$,
we remove $i$ and $j$ from the list, and insert a new pseudo-parton with total
momentum $p_{i}^{\mu}+p_{j}^{\mu}$.

  * •

Repeat the procedure until only one pseudo-parton remains in the list: this
pseudo-parton forms its own candidate jet.

We have now a list of candidate jets, among which we need to find the final-
state macro jet, which contains all the radiation from the original final-
state quark that had a momentum $n_{2}^{\mu}$. Each jet momentum $p_{j}^{\mu}$
can be written as

$p_{j}^{\mu}=\alpha_{j}P^{\mu}+\beta_{j}n_{2}^{\mu}+p_{\perp,j}^{\mu}\,,$ (84)

where $P^{\mu}$ is the incoming proton momentum. We label as final-state macro
jet the one with the largest $\beta_{j}$ value, i.e. the one that retains the
largest fraction of the light-cone component of the original final-state quark
in the Breit frame. This jet can easily be found by searching for the jet
whose momentum $p_{j}^{\mu}$ yields the largest value of $p_{j}\cdot P$.
121212 In the runs of Sec. 5.3, we have tagged as final-state macro jet the
last pseudo-parton that becomes a candidate jet. We have verified that, at NLL
accuracy, this leaves the definition of the event shapes
$\mathcal{M}_{j,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$ and
$\mathcal{S}_{j,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$ (see Sec. 5.3 and C) unchanged, so it does not alter the results shown
in Fig. 9. We however believe that choosing this naive definition of the macro
jet, instead of looking for the jet with the largest $\beta_{j}$, can lead to
contributions which are not infrared safe beyond NLL, as the energy of the
final-state macro jet develops a soft tail, due to very energetic initial-
state radiation that moves the hard final-state quark to the remnant
hemisphere. It is easy to see that this soft tail is at the most an NNLL
correction for these event shapes.

After this clustering algorithm, the event contains a collection of initial-
state/beam jets and one final-state macro jet. To calculate our event shapes
we need to define the primary Lund-plane variables. For each initial-
state/beam jet, we define the primary Lund-plane coordinates as

$\eta_{i}=-\frac{1}{2}\ln\frac{1-\cos\theta_{i}}{1+\cos\theta_{i}}\,,\qquad
k_{t,i}=E_{i}\sin\theta_{i}\,.$ (85)

Note that these exactly correspond to the pseudorapidity and transverse-
momentum of $i$ in the Breit frame. The negative sign of the pseudorapidity is
due to the convention to orientate the incoming beam along the negative $z$
axis in the Breit frame.

For the final-state macro jet a different procedure is adopted. We consider
the clustering sequence starting from the last recombination, and consider the
energies of the two pseudo-jets that were combined, $E_{i}$ and $E_{j}$, where
$E_{i}>E_{j}$. The softer pseudojet, which we have labelled with $j$, is then
promoted to be a jet, with Lund coordinates

$\displaystyle\eta_{j}=\frac{1}{2}\ln\frac{1-\cos\theta_{ij}}{1+\cos\theta_{ij}}\,,\quad
k_{t,j}=E_{j}\sin\theta_{ij}\,.$ (86)

We iterate this procedure for the most energetic pseudojet.

The procedure is terminated when the most-energetic pseudojet has no further
children, i.e. it is a single parton (notice that this parton is not
associated with any jet). The secondary Lund plane is defined by the jet
(found inside the final-state macrojet, or among the collection of initial-
state jets) that has the largest value for the Lund variable $k_{t}$.

##  Appendix C Resummation formulae for continuously-global observables

We consider a continuously-global observable Banfi:2004yd  that in the soft-
collinear limit behaves as

$O\sim\frac{p_{T}}{Q}e^{-\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}|\eta|}.$ (87)

The cumulative cross section at NLL accuracy can be written as

$\Sigma_{\rm
NLL}(O<e^{L})=\exp\left[-Lg_{1}(\bar{\lambda})+g_{2}(\bar{\lambda})\right],$
(88)

with $\bar{\lambda}=-\alpha_{s}b_{0}L$ (and $L<0$) and

$\displaystyle b_{0}=\frac{11C_{A}-4n_{f}T_{R}}{12\pi}.$ (89)

The $g_{1}$ function contains the LL terms. For processes with two hard legs
it reads

$\displaystyle g^{\beta_{\text{obs}}=0}_{1}$ $\displaystyle=\frac{C_{i}}{\pi
b_{0}\bar{\lambda}}\left(2\bar{\lambda}+\ln(1-2\bar{\lambda})\right),$ (90a)
$\displaystyle
g^{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}\neq 0}_{1}$ $\displaystyle=\frac{C_{i}}{\pi
b_{0}\bar{\lambda}\beta_{\text{obs}}}\left((1+\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}-2\bar{\lambda})\ln\left(1-\frac{2\bar{\lambda}}{1+\beta_{\text{obs}}}\right)-(1-2\bar{\lambda})\ln(1-2\bar{\lambda})\right),$
(90b)

where $C_{i}$ is the Casimir factor ($C_{i}=C_{F}$ for quark radiators,
$C_{A}$ for gluon radiators). The NLL term $g_{2}$ can be written as

$g^{\beta_{\text{obs}}}_{2}=\ln\frac{f_{i}\left(x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS},Q^{2}e^{2L/(1+\beta_{\text{obs}})}\right)}{f_{i}\left(x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS},Q^{2}\right)}+\bar{g}^{\beta_{\text{obs}}}_{2}+\ln\mathcal{F}_{\text{obs}}\,,$
(91)

where $f_{i}$ is the PDF of the incoming parton. We use
$\bar{g}^{\beta_{\text{obs}}}_{2}$ to denote the universal
$\beta_{\text{obs}}$-dependent term originating from soft and hard-collinear
emissions, and $\mathcal{F}_{\text{obs}}$ is an observable-dependent
correction. We have

$\displaystyle\bar{g}^{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=0}_{2}=$ $\displaystyle\frac{C_{i}}{\pi
b_{0}^{2}}\Bigg{[}b_{0}B^{(1)}_{i}\ln(1-2\bar{\lambda})-\frac{K}{2\pi}\left(\frac{2\bar{\lambda}}{1-2\bar{\lambda}}+\ln(1-2\bar{\lambda})\right)$
(92a)
$\displaystyle\quad\quad\quad\quad\quad\quad+\frac{b_{1}}{b_{0}}\left(\frac{2\bar{\lambda}+\ln(1-2\bar{\lambda})}{1-2\bar{\lambda}}+\frac{1}{2}\ln^{2}(1-2\bar{\lambda})\right)\Bigg{]},$
$\displaystyle\bar{g}^{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}\neq 0}_{2}=$ $\displaystyle\frac{C_{i}}{\pi
b_{0}^{2}\beta_{\text{obs}}}\Bigg{[}\frac{K}{2\pi}\left(\ln(1-2\bar{\lambda})-(1+\beta_{\text{obs}})\ln\left(1-\frac{2\bar{\lambda}}{1+\beta_{\text{obs}}}\right)\right)$
(92b)
$\displaystyle\quad\quad\quad\quad+b_{0}\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}B^{(1)}_{i}\ln\left(1-\frac{2\bar{\lambda}}{1+\beta_{\text{obs}}}\right)-\frac{b_{1}}{b_{0}}\Bigg{(}\frac{1}{2}\ln^{2}(1-2\bar{\lambda})+\ln(1-2\bar{\lambda})$
$\displaystyle\quad\quad\quad\quad-\frac{1}{2}(1+\beta_{\text{obs}})\ln^{2}\left(1-\frac{2\bar{\lambda}}{1+\beta_{\text{obs}}}\right)-(1+\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs})\ln\left(1-\frac{2\bar{\lambda}}{1+\beta_{\text{obs}}}\right)\Bigg{)}\Bigg{]},$

with

$\displaystyle
b_{1}=\frac{17C_{A}^{2}-10C_{A}n_{f}T_{R}-6C_{F}n_{f}T_{R}}{24\pi^{2}}\,,$
$\displaystyle\qquad
K=\left(\frac{67}{18}-\frac{\pi^{2}}{6}\right)C_{A}-\frac{10}{9}n_{f}T_{R}\,,$
(93) $\displaystyle B^{(1)}_{q}=-\frac{3}{4}\,,$ $\displaystyle\quad
B^{(1)}_{g}=\frac{-11C_{A}+4n_{F}T_{R}}{12C_{A}}\,.$

Finally, we come to the observable-dependent correction. For the max-type
observable
$M_{j,\,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$ (defined in eq. (61)), we have $\ln\mathcal{F}_{\text{obs}}=0$, while
for
$S_{p/j,\,\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}}$, (eq. (61)), this correction reads

$\displaystyle\ln\mathcal{F}_{S_{j,\,\beta_{\text{obs}}}}=-\gamma_{E}R^{\prime}(\bar{\lambda})-\ln\Gamma\left(1+R^{\prime}(\bar{\lambda})\right),$
(94)

with $R^{\prime}(\bar{\lambda})$ defined as
$\partial_{L}\left(Lg_{1}(\bar{\lambda})\right)$,

$\displaystyle
R^{\prime}_{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=0}(\bar{\lambda})$ $\displaystyle=\frac{4C_{i}}{\pi
b_{0}}\frac{\bar{\lambda}}{1-2\bar{\lambda}}\,,$ (95a) $\displaystyle
R^{\prime}_{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}\neq 0}(\bar{\lambda})$ $\displaystyle=\frac{2C_{i}}{\pi
b_{0}\beta_{\text{obs}}}\left[\ln\left(1-\frac{2\bar{\lambda}}{1+\beta_{\text{obs}}}\right)-\ln(1-2\bar{\lambda})\right].$
(95b)

In Sec. 5.3 we also considered standard DIS event shapes, that have the
property of being continuously global, and can be defined considering only
partons in the current hemisphere: the thrust with respect to the photon axis
and normalised to $Q/2$, $\tau_{zQ}$ (eq. (62c)), and the broadening with
respect to the photon axis, normalised either to the energy in the current
hemisphere $B_{zE}$ (eq. (62a)) or to $Q$, $B_{zQ}$ (eq. (62b)). The NLL
analytic predictions for $\tau_{zQ}$ were computed for the first time in Ref.
Antonelli:1999kx , and they are identical to those for
$S_{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=1}$. The two definitions of broadening $B_{zE}$ and $B_{zQ}$ are
equivalent at NLL, and the analytic prediction can be read from Ref.
Dasgupta:2001eq  to read

$\displaystyle\ln\Sigma(B_{z}<e^{L})=-g_{1}^{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=0}(\bar{\lambda})\,L+\bar{g}^{\beta_{\text{obs}}=0}_{2}(\bar{\lambda})+\ln\mathcal{F}_{B}(\bar{\lambda})+\ln\frac{f_{i}\left(x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS},Q^{2}e^{2L}\right)}{f_{i}\left(x_{\rm\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}DIS},Q^{2}\right)}\,,$
(96)

where $B_{z}$ can be either $B_{zE}$ or $B_{zQ}$. The observable-dependent
factor $\ln\mathcal{F}_{B}$ had an analytic form that reads

$\displaystyle\ln\mathcal{F}_{B}=$ $\displaystyle-\frac{1}{2}(\ln
2+2\gamma_{E})R^{\prime}+\ln\Lambda(R^{\prime})-\ln\Gamma(1+R^{\prime})\,,$
(97a) $\displaystyle\Lambda(R^{\prime})=$
$\displaystyle\int_{0}^{\infty}\frac{ydy}{(1+y^{2})^{3/2}}\left(\frac{y(1+\sqrt{1+y^{2}})}{8}\right)^{-R^{\prime}/2}$
$\displaystyle=$ $\displaystyle
8^{R^{\prime}/2}\Gamma(1-R^{\prime}/4)\Gamma(1+R^{\prime})\frac{{}_{2}{F}_{1}(3R^{\prime}/4,1+R^{\prime},2+3R^{\prime}/4;-1)}{\Gamma(2+3R^{\prime}/4)},$
(97b)

where
$R^{\prime}=R^{\prime}_{\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=0}(\bar{\lambda})$ defined in eq. (95a) and ${}_{2}{F}_{1}(a,b,c;z)$ is
the hypergeometric function.

##  Appendix D Spin correlations

Spin correlations must be included to reproduce the correct azimuthal
structure of strongly angular-ordered collinear splittings. This can be
achieved by for example using the Collins-Knowles algorithm Collins:1987cp ;
Knowles:1987cu ; Knowles:1988vs ; Knowles:1988hu , which was also recently
applied to the PanScales final-state and initial-state showers
Karlberg:2021kwr ; vanBeekveld:2022zhl . In Ref. Hamilton:2021dyz  the
algorithm was further improved to include for the first time the treatment of
the dominant leading-colour soft azimuthal correlations. In this appendix, we
show that the algorithm introduced in Ref. vanBeekveld:2022zhl  can be applied
straightforwardly to DIS.

We consider the shower effective matrix element for the production of two
extra partons $i$ and $j$, and compare the result against the analytic
expectation as a function of the azimuthal angle $\Delta\psi_{ij}$ between the
planes spanned by two emissions $i$ and $j$. At $\mathcal{O}(\alpha_{s}^{2})$,
the differential cross section can be written as

$\displaystyle\frac{d\sigma}{d\Delta\psi_{ij}}\propto
a_{0}\left(1+\frac{a_{2}}{a_{0}}\cos(2\Delta\psi_{ij})\right)=a_{0}\left(1+A(z_{i})B(z_{j})\cos(2\Delta\psi_{ij})\right)\,,$
(98)

where $\psi_{ij}$ is the azimuthal difference between the plane defined by the
primary and secondary splittings with light-cone momentum fraction $z_{i}$ and
$z_{j}$, respectively. The values of $a_{0}$ and $a_{2}$ depend on the type of
branching, and are a function of $z_{i}$ and $z_{j}$. In the absence of spin
correlations, the ratio $a_{2}/a_{0}$ would be 0.

Figure 14: Size of the ratio between the two Fourier coefficients
$a_{2}/a_{0}$, defined in eq. (98) at $\mathcal{O}(\alpha_{s}^{2})$ for
collinear splittings, for the
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=0$) shower. The Feynman diagrams indicate the sequence of splittings
under consideration. The black lines indicate constant values for this ratio,
and are obtained using the analytic predictions. The (maximum) deviation of
the analytic prediction and the shower is given by ${\rm
D}\big{[}\frac{a_{2}}{a_{0}}\big{]}\leq{\rm
max}\big{|}\big{(}\frac{a_{2}}{a_{0}}\big{)}_{\rm
PS}-\big{(}\frac{a_{2}}{a_{0}}\big{)}_{\rm ME}\big{|}$, with
$\big{(}\frac{a_{2}}{a_{0}}\big{)}_{\rm PS}$ the shower prediction and
$\big{(}\frac{a_{2}}{a_{0}}\big{)}_{\rm ME}$ the analytic prediction for the
matrix element.

In Fig. 14 we illustrate the ratio $a_{2}/a_{0}$ across several values of
$z_{1}$ and $z_{2}$ for collinear-splittings involving an intermediate gluon.
In particular, in the left panel we show the case in which the first emission
is a final-state gluon, that further branches into a $q\bar{q}$ pair. In the
right panel we instead consider the backward evolution of the initial-state
quark into an initial-state gluon ($g_{I}$), and then an emission of a final-
state gluon very collinear to the incoming beam. We only show the predictions
obtained with the
PanGlobal($\beta_{\mathchoice{\displaystyle}{\scriptstyle}{\scriptscriptstyle}{\scriptscriptstyle}\rm
obs}=0$) shower, but identical results can be obtained by considering the
other PanScales showers. The differences between the shower predictions and
the analytic expectation were always smaller that $10^{-6}$, confirming our
implementation at fixed-order accuracy.

## References

  * (1) M. Bahr et al., _Herwig++ Physics and Manual_ , _Eur. Phys. J._ C58 (2008) 639–707, [0803.0883].
  * (2) J. Bellm et al., _Herwig 7.2 release note_ , _Eur. Phys. J. C_ 80 (2020) 452, [1912.06509].
  * (3) S. Gieseke, P. Stephens and B. Webber, _New formalism for QCD parton showers_ , _JHEP_ 12 (2003) 045, [hep-ph/0310083].
  * (4) G. Bewick, S. Ferrario Ravasio, P. Richardson and M. H. Seymour, _Logarithmic accuracy of angular-ordered parton showers_ , _JHEP_ 04 (2020) 019, [1904.11866].
  * (5) G. Bewick, S. Ferrario Ravasio, P. Richardson and M. H. Seymour, _Initial state radiation in the Herwig 7 angular-ordered parton shower_ , _JHEP_ 01 (2022) 026, [2107.04051].