Spin correlations in final-state parton showers and jet observables

As part of a programme to develop parton showers with controlled logarithmic accuracy, we consider the question of collinear spin correlations within the PanScales family of parton showers. We adapt the well-known Collins-Knowles spin-correlation algorithm to PanScales antenna and dipole showers, using an approach with similarities to that taken by Richardson and Webster. To study the impact of spin correlations, we develop Lund-declustering based observables that are sensitive to spin-correlation effects both within and between jets and extend the MicroJets collinear single-logarithmic resummation code to include spin correlations. Together with a 3-point energy correlation observable proposed recently by Chen, Moult and Zhu, this provides a powerful set of constraints for validating the logarithmic accuracy of our shower results. The new observables and their resummation further open the pathway to phenomenological studies of these important quantum mechanical effects.


Introduction
One of the most striking properties of Quantum Mechanics is that the spin angular momentum of two or more particles can be created in an entangled state [1]. As a consequence, when measuring the spin of the individual particles, or more generally the angular distributions of particle decays and branchings, long-distance correlations will be found depending on the degree of entanglement. At colliders, spin correlations are most widely studied in the context of heavy particle decays (see e.g. Refs. [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]). However they play a significant role also in the pattern of QCD branchings that occur in jet fragmentation, as studied for example at LEP [18,19]. The quantum mechanical nature of the problem is reflected in the need to sum coherently over the spin states of intermediate particles in the jet fragmentation, similarly to the need to sum coherently over the spins of an electron-positron pair in a e-mail: alexander@erikkarlberg.dk (corresponding author) an Einstein-Podolsky-Rosen (EPR) experiment [1]. It was recognised long ago [20][21][22] that the core tools for simulating jet fragmentation, i.e. parton showers, should incorporate such effects.
In this article we consider collinear spin correlations in the context of the PanScales programme [23][24][25] of QCD parton-shower development. The core aim of the programme is to develop parton showers with well-understood logarithmic accuracy. A first step on that path is to achieve socalled next-to-leading logarithmic (NLL) accuracy. There are many senses in which a shower can be NLL accurate and we choose two broad criteria [24]. The logarithmic phase-space for QCD branching involves two dimensions, corresponding to the logarithms of transverse momentum and of angle, which can conveniently be represented using Lund diagrams [26]. To claim NLL accuracy, we firstly require a shower to reproduce the correct matrix element for any configuration of emissions where all branchings are well separated from each other in the Lund diagram. Secondly, for all observables where suitable resummations exist (e.g. event shape distributions), the shower should reproduce the resummation results up to and including terms α n s L n . Here L is the logarithm of the value of the observable and terms α n s L n are NLL in a context where α n s L n+1 terms exponentiate [27]. Until now, the PanScales shower development has been based on unpolarised splitting functions, as is common in dipole and antenna showers. According to our accuracy criterion, however, spin correlations are a crucial part of NLL accuracy, because in configurations with successive branchings at disparate angles, they are required in order to reproduce the correct azimuthal dependence of the matrix elements. The core purpose of this article, therefore, is to start the implementation of spin correlations, specifically as concerns nested collinear emissions. The algorithm that we will adopt is based on the well-established proposal by Collins [21], used notably in the Herwig series of angular ordered showers [28][29][30][31]. Our adaptation for the PanScales antenna and dipole showers bears similarities with the implementation for the Herwig7 dipole shower by Richardson and Webster [32,33] (for work in other shower frameworks, see Refs. [34][35][36][37]). One class of configuration that is not addressed by the Collins algorithm is that of two or more commensurateangle energy-ordered soft emissions followed by a collinear splitting of one or more of them. Strictly these configurations should also be addressed for NLL accuracy, however we defer their study to future work.
Whilst the techniques to implement parton-shower spin correlations are relatively well established, our philosophy is that a parton shower implementation is only complete if one has a set of observables, resummations and associated techniques to validate the implementation. The PanScales shower development in Refs. [24,25] has been able to draw on a rich set of event shapes and associated resummed calculations. However for testing spin correlations there was no such set of observables. Recently, while our work was in progress, Chen et al. [38] introduced and resummed a 3-point energy-energy correlator that is spin-correlation sensitive. Here we introduce a new set of spin-correlation sensitive observables based on Lund declustering. We extend the MicroJets collinear resummation code [39,40] to allow for the treatment of azimuthal structure, so as to obtain α n s L n numerical predictions for all of these observables (in the process, confirming the analytic results of Ref. [38]). These new observables, and our study of their properties, are of potential interest also in their own right for practical measurements of spin correlations in jets.
The paper is structured as follows. In Sect. 2 we give details of the algorithm used to implement spin correlations in the PanScales showers. We introduce an azimuthal observable based on the Lund plane picture in Sect. 3, and discuss general features of spin correlations in the strictly collinear limit. 1 In Sect. 4 we validate our implementation and show that it achieves NLL accuracy. In Sect. 5 we conclude.

The Collins-Knowles algorithm and its adaptation to dipole showers
An efficient algorithm to include spin correlations in angularordered MC generators was proposed by Collins [21] for 1 Reference [24] also used azimuthal observables for testing showers, however for the case of emissions with commensurate k t that are widely separated in angle. That is a region dominated by independent soft emission and is free of spin correlations. The recoil effects in traditional dipole showers that cause that region to be incorrect at NLL accuracy can also introduce azimuthal correlations that obfuscate the genuine spin correlations that we discuss here, as is visible in Appendix D.1, though an understanding of the logarithmic structure associated with this effect would require further work. This issue was discussed also in Ref. [32].
final-state showers and subsequently extended by Knowles [22,41,42] to include initial-state radiation and backwards evolution. It is based on the factorisation of the tree-level matrix element in the collinear limit, and can be generalised to include spin correlations in the matching of the parton shower to a hard matrix element, between initial and final state radiation, as well as in decays [43]. The Collins-Knowles algorithm can be readily applied to showers that first generate a full branching tree with intermediate virtualities and momentum fractions at each stage of the splitting, and then only subsequently assign azimuthal angles and reconstruct the full event kinematics. However, in dipole showers the azimuthal angle needs to be chosen at each stage of the branching, because it affects the phase space for subsequent branchings. This requires a reordering of the steps in the Collins-Knowles algorithm. Furthermore, in an angular-ordered shower, it is straightforward to identify the mapping between the shower kinematics and the azimuthal angle as needed in the Collins-Knowles algorithm. In dipole (and antenna) showers the corresponding mapping can be less straightforward. In this section we therefore introduce a modified version of the Collins-Knowles algorithm that is applicable to any parton shower, including those of the dipole or antenna kind. We first develop a collinear branching formalism in terms of boost invariant spinor products and then show how to apply these in the context of dipole and antenna showers.
Our work here is not the first to adapt the Collins-Knowles algorithm to showers with dipole-like structures, see for example Refs. [32][33][34][35]. Our approach is inspired by and bears a number of similarities with that by Richardson and Webster [32,33], as used in the Herwig7 [31] framework. That algorithm continuously boosts between the lab frame and frames specific to each individual collinear splitting, where individual Collins-Knowles steps may be applied directly. In our implementation of the algorithm no such boost is necessary, as our expressions are formulated in terms of boost-invariant spinor products.

Collinear branching amplitudes
To determine the appropriate azimuthal distribution of shower branchings, the Collins-Knowles algorithm makes use of collinear branching amplitudes M λ a λ b λ c a→bc for a splitting a → bc with spin labels λ = ±1. As an example, consider the situation of Fig. 1, where an unpolarised parton 0 emits a gluon 2, which subsequently splits into a q q pair, 2 → 34. In the collinear limit, the azimuthal distributions are determined by the factorised matrix element where summation over repeated spin indices is implied. 2 Ultimately we will test our approach in the PanScales shower framework, which currently only supports massless particles in its kinematic maps. Accordingly, we work in the massless quark limit, though the extension to massive quarks is straightforward, as discussed in Ref. [32]. In the massless quark limit, the branching amplitudes may be written as where the functions F λ a λ b λ c a→bc (z), listed in Table 1, are the (colour-stripped) helicity-dependent Altarelli-Parisi splitting amplitudes that depend on the collinear momentum fraction z carried by parton b. They are related to the usual unregularised splitting functions through 3 The function S τ ( p b , p c ) is a spinor product, where the label τ = ±1 indicates the sign of the complex phase associated with the spinor product. It is given by The derivation of these branching amplitudes can be found in Appendix A, together with our convention for spinor products.
Inserting an explicit expression for the spinor product, Eq. (2) may alternatively be written as where φ is the azimuthal angle as defined in a reference frame where the parent is along a specific (e.g. z) direction. Eq. (5) is useful in the context of a parton shower only when this azimuthal angle is used to parameterise its phase space directly. This can be the case when the azimuthal variable φ for the shower's phase space generation is defined with respect to a fixed azimuthal reference direction. However, in dipole or antenna showers, φ usually represents the azimuthal direction of the transverse momentum with respect to the plane of the dipole or antenna parent partons. As the Pan-Scales showers are of the dipole/antenna type, we will use Eq.
(2) directly, and evaluate the spinor products numerically using Eq. (39). This choice guarantees that the implementation of the algorithm remains independent of the type of parton shower it is applied to. However one could also choose to track the relation between the shower azimuths and the collinear azimuths, as done in Refs. [32,33].

The algorithm
The purpose of the Collins-Knowles algorithm is to distribute the azimuthal degrees of freedom according to Eq. (1) for an arbitrary number of collinear branchings, while maintaining linear complexity in the number of particles. To that end, a binary tree is tracked, where the nodes correspond to collinear shower branchings. Figure 2 shows an example of a shower history and the corresponding binary tree. The algorithm is initialised from the hard scattering amplitude M λ 1 λ 2 hard , where the outgoing hard partons 1, 2 have spin indices λ 1 , λ 2 . The formalism is not limited to a particular number of final-state partons, but for readability we restrict ourselves to a two-body hard scattering. Each particle i (whether final or intermediate) is associated with a so-called decay matrix, D λ i λ i i with two spin indices. At the start of the shower the decay matrices are initialised to D The core of the algorithm consists of the rules for generating an azimuthal angle and adding nodes to and updating the binary tree, as well as a subset of the decay matrices, each time the shower adds an emission to the event. The shower first selects a branching dipole or antenna and generates an ordering scale and momentum-sharing variable 4 according to the regular spin-summed shower dynamics, in which the azimuthal dependence does not appear. 5 The branching parton will correspond to one of the terminal nodes in Fig. 2, as discussed in further detail below, and we refer to the node as a n , where the n index labels the depth of the node within the 4 The exact definitions of which depend on the shower implementation at hand. 5 This may no longer be the case when accounting for subleading colour effects, for example in the NODS method of Ref. [25]. That specific algorithm may reject an emission after its azimuth has been chosen. When combining it with our spin-correlation approach, we first choose the azimuth according to the spin-correlation algorithm, and then apply the colour rejection step. Table 1 The helicity-dependent Altarelli-Parisi splitting amplitudes F λa λ b λc tree. Then, to determine the azimuthal angle, the algorithm proceeds with the following rejection-sampling procedure: 1. Compute the spin-density matrix ρ λ an λ an a n as follows: (a) Starting from parton a n , trace the binary tree back to the root parton node a 0 ∈ [1,2]. This results in a sequence a 0 . . . a n of parton indices, along with a sequence of complementing parton indices b 0 . . . b n such that b 0 is the other root parton and a i and b i have common parent node a i−1 . (b) Compute the spin-density matrix for the root node, where the denominator is the trace of the numerator.
(c) Iteratively for i ∈ {1, . . . , n}, compute 2. Repeat the following until a value of φ is accepted: (a) Sample a value of φ uniformly and, using the other shower variables, construct the post-branching momenta p a n+1 and p b n+1 using the usual kinematic mapping of the parton shower.
(b) Compute the branching amplitude M λ an λ a n+1 λ b n+1 a n →a n+1 b n+1 using Eq. (2). (c) Compute the acceptance probability where N is a φ-independent normalisation factor to ensure p accept < 1. 6 (d) Accept φ with probability p accept .
3. Update the binary tree (a) Insert new nodes a n+1 and b n+1 with parent node a n and initialise their decay matrices to a Kronecker delta. (b) Store M λ an λ a n+1 λ b n+1 a n →a n+1 b n+1 in node a n . (c) Iteratively for i ∈ {n, . . . , 0}, recompute and store the updated decay matrices Note that the identification of the parton that branches is not without subtleties. We have assumed that it is always one of the terminal nodes in Fig. 2. To understand why this is a non-trivial choice, suppose that we have a dipole stretching between particles 3 and 7. When the left-hand end of the dipole emits a gluon (which we label 9), our spin-correlation algorithm always views this as a g → gg splitting of particle 3. However if the emitted gluon is at large angle relative to the 34 splitting, i.e. θ 93 θ 34 the gluon effectively sees the coherent charge of 3 and 4 and could more properly be viewed as being emitted from 2 (which unlike 3 is a quark). Because of the interplay between the shower ordering variable and 6 One can make use of the fact that the spin-density matrix is hermitian and has trace 1 to find, for instance emission kinematics, this occurs only for situations in which 9 is soft relative to particle 3, and also soft relative to any of the parents of 3. Inspecting Table 1, one sees that soft gluon emission (the z → 1 limit) leads to splitting amplitudes that are independent of the flavour of the parent, a, and that are non-zero only for λ a = λ b , i.e. they are diagonal in the spin space relating the parent and its harder offspring. This means that in the limit where emission 9 could conceivably have been emitted from 2, it is immaterial whether we actually view it as being emitted from 2 or instead organise the tree as if it had been emitted by 3. The latter is considerably simpler and so it is the solution that we adopt.

Collinear spin correlations: expectations and measurement strategy
In this section, we start (Sect. 3.1) by examining how the spin correlations translate into azimuthal correlations between the planes of separate collinear branchings, both within a single jet and across pairs of jets. We do so at fixed order, O α 2 s , where it is trivial to define the observables. We then propose (Sect. 3.2) a set of observables that are suitable for use at all orders. They exploit a Lund diagram [26] representation of individual jets [44]. Next (Sect. 3.3), we recall the definition of the EEEC spin-sensitive observable, which was proposed and resummed in Ref. [38]. Finally (Sect. 3.4), we use these observables to study the impact on the azimuthal correlations coming from the all-order resummation of collinear spincorrelation effects.

Azimuthal structure
Each collinear branching in an event can be associated with the plane that contains the momenta of the two offspring partons. The simplest observable one may think of to study spin correlations is the azimuthal difference, ψ, between the planes defined by two distinct branchings. Here we will consider two broad cases: intra-jet correlations, i.e. between the planes of two branchings within a single jet, for example between the plane of the 1 → 56 splitting and the plane of the 6 → 78 splitting in Fig. 2; and inter-jet correlation, i.e. between the planes of two splittings in separate jets, for example between the plane of the 1 → 56 splitting and the plane of the 2 → 34 splitting in Fig. 2. 7 We will refer to the two azimuthal differences as ψ 12 and ψ 11 where the 1 and 2 labels refer to the first and second splitting within a given jet and the 1 label refers to the first splitting in a distinct jet. The ψ 12 and ψ 11 observables are straightforward to define at O α 2 s relative to the hard scattering and it is this situation that we will concentrate on here. In the ψ 12 case, the splitting planes and the azimuthal angle between them are illustrated in Fig. 3. At second order in the coupling, and in the collinear limit, the cross sections differential in the intra-and inter-jet observables ψ 12 , respectively ψ 11 , take the simple form (see e.g. section 2.3 in Ref. [32], or the example calculated in Appendix B) where the coefficients a 0 and a 2 depend on the observable, the final state under consideration, and the momentum fractions associated with the first (z i ) and second splitting (z j ). If the splittings are restricted to have opening angles greater than e −|L| , with |L| 1, and ln z 1 and ln z 2 are both of order 1, the a 0 and a 2 coefficients are dominated by terms α 2 s L 2 , i.e. they belong to the singlelogarithmic set of terms that we aim to control for NLL accuracy. For large L, at this order, the ratio a 2 /a 0 is independent of α s L, and so Eq. (11) can always be written in terms of the functions A(z) and B(z) given in Table 2.
In Fig. 4 we show a contour plot of the ratio a 2 /a 0 = A(z 1 )B(z 2 ), for our intra-jet observable, as a function of z 1 (x-axis) and z 2 (y-axis) for a quark-initiated jet. In this case there are only two non-trivial final states to consider, namely q → qg(g → q q ) and q → qg(g → gg), while channels such as q → q(q → qg)g that do not involve an intermediate gluon have vanishing spin correlations. In the case of q → qg(g → q q ), we see that the coefficient a 2 is negative, corresponding to an enhancement when the q → qg and g → q q planes are perpendicular. The ratio a 2 /a 0 peaks at −100% around z 1 = 0 and z 2 = 0.5, i.e. when the gluon is soft and the quark-antiquark pair share the energy equally. 8 The spin correlations stay large even for moderate values of z 1 and z 2 although they vanish completely for z 1 → 1 (soft quark emerging from the q → qg splitting) and for z 2 → 0 or z 2 → 1 (either of the secondary q orq becomes soft). Similarly, for q → qg(g → gg), the ratio peaks at z 1 = 0 and z 2 = 0.5 but the correlation has the opposite sign to the g → qq case, which implies an enhancement when the q → qg and g → gg splittings are in the same plane. The magnitude of a 2 /a 0 is substantially smaller, with a peak at 1/9. The spin correlations fall off more sharply than in the g → q q case, in particular as a function of z 2 , and also vanish for z 1 → 1 and for z 2 → 0 or z 2 → 1.
The overall picture is very similar for a gluon-initiated jet, shown in Fig. 5, although the dependence on z 1 is moderately stronger for both g → q q and g → gg secondary splittings, which can be understood from Table 2.
Turning now to the correlations between azimuthal angles in two distinct jets, they are zero for a γ * → qq, but nonzero for H → g 1 g 1 , so we consider only the latter. In this case we examine the azimuthal difference ψ 11 between the splitting planes of the two primary emissions. There are three possible final states given by Fig. 6 we show the ratio a 2 /a 0 as a function of the energy fractions z 1 and z 1 . In all three cases the ratio is peaked when z 1 = z 1 = 0.5 and is largest in magnitude when both gluons split into quark-antiquark pairs. When only one gluon splits into a quark-antiquark pair the ratio becomes negative and −11.1% around the peak. When both gluons split into gluons the spin correlations almost vanish. In all three cases the ratio vanishes when either of the energy fractions approach 0 or 1.

Definition of spin-sensitive Lund observables
The observables discussed so far start at O(α 2 s ) relative to the hard scattering, but care has to be taken in order to define them in an infrared safe way beyond that order. To facilitate a definition which can also be directly applied in experimental analyses, we make use of the procedure from Ref. [44] to build Lund diagrams [26] from individual jets. 9 We remind the reader that the Lund jet plane is constructed through a declustering of a C/A-jet [45,46]. The primary Lund plane is constructed by following the harder branch at Table 2 The functions A(z) and B(z) entering Eq. (11) for the intraand inter-jet observables ψ 12 , ψ 11 , for all channels at second order. See e.g. Appendix B for an exemplified derivation. Here, for a 1 → 23 splitting, the variable z is the momentum fraction carried away by parton 2. Representative diagrams at O α 2 s are shown on the right. In those diagrams, the splitting in black is the one associated with the cor-responding function A(z) or B(z). The partons shown in grey serve as an example of what the rest of the branching history can look like, but they do not matter in choosing the function itself, since the contributions factorise in the final result, Eq. (11)

Δψ 12
Primary splitting Size of the spin correlations, at fixed order O(α 2 s ), for a quarkinitiated jet. We consider the azimuthal difference ψ 12 between the splitting planes of a primary (with gluon momentum fraction z 1 ) and secondary (z 2 ) branching, separated by channel. The colour scale indicates the relative size a 2 /a 0 of the correlations, where ψ 12 is dis-tributed proportionally to a 0 + a 2 · cos(2 ψ 12 ). Note the use of z 1 as the gluon momentum fraction in the q → qg splitting, while for amplitude expressions in the text, z often refers to the quark momentum fraction Δψ 12 H → g 1 g 1 , g 1 → gg (g → gg)  Δψ 11 H → g 1 g 1 , g 1 → gg, g 1 → gg largest for the quark-only final-state (top left) which peak at +100%, mid-range for the mixed quark-gluon final state (top right) peaking at −11.1%, and much smaller for the all-gluon final state (bottom), with a peak at +1.23% each step of the declustering. If at some point one instead follows a softer branch and then all its subsequent harder branches, those subsequent branches populate a secondary Lund plane. The procedure that we describe here can be applied at hadron colliders and e + e − colliders (or also in the decay of heavy particles at hadron colliders, for example in their rest frame). The first branching within a jet is identified by taking all declusterings on the primary Lund plane that satisfy z larger than some z cut , and selecting the one with largest relative k t . 10 The second branching within a jet is identified by taking all declusterings on the secondary Lund plane associated with the first branching, and then again identifying those that satisfy z larger than some z cut and selecting the one with largest relative k t . One may choose to use different z cut values for the first and second branchings within a single jet. 11 Given a choice of z cut , one may then study the results differentially in the two z values, so as to obtain plots similar to Figs. 4, 5 and 6. Here, however, we will consider just results integrated over z. Note that in Sect. 3.1, we investigated the structure for 0 < z < 1, while in the Lund construction we always have z ≤ 1/2, because z is the momentum carried by the softer offspring (and the secondary Lund plane is the one associated with further emission from the softer offspring). So, for example, in Sect. 3.1 we could meaningfully discuss a q → qg splitting where the quark was soft (z 1 close to 1), followed by a splitting of the gluon. In contrast, in the Lund declustering picture, the quark in such a case will be assigned to the secondary plane, and the second splitting that we study will necessarily be a further q → qg splitting. The reason for adopting this procedure in Lund declustering is that experimentally one cannot observe the flavour of the underlying parton. 12 10 In the pp case, z, k t and the azimuth φ for a branching were defined in Ref. [44]; in our e + e − studies here, for an i → jk branching, we use z = min(| p j |, | p k |)/(| p j | + | p k |), k t = min(| p j |, | p k |) sin θ jk and we discuss the definition of the azimuths below. 11 Instead of selecting the highest-k t emission that satisfies z > z cut , one could instead take the first that appears in the declustering. That would have made the procedure more similar to the modified mass-drop tagger [47] or SoftDrop with β = 0 [48]. At the logarithmic accuracy that we discuss here, both procedures should yield equivalent results. Selecting the highest-k t emission without a z cut would be similar to Dynamical Grooming [49] and would involve a double rather than single logarithmic resummation structure [50]. The use of two declusterings is also part of a range of top taggers, e.g. Refs. [51][52][53][54][55], theoretical aspects of which are further discussed in Ref. [56]. 12 We will still show results below classified according to the flavour channel. This is useful in terms of diagnostics, and it is possible because the parton shower does contain information about flavours. Strictly speaking, that flavour information for massless quarks is infrared unsafe [57] within the Cambridge/Aachen algorithm. The infrared-safe flavour algorithm of Ref. [57] cannot be applied in conjunction with Lund declustering, because it adapts the k t jet algorithm rather than the Cam-To track the azimuthal angles in the e + e − case within the Lund declustering, we adopt a procedure that differs from that introduced in Ref. [24] (supplemental material), because we have found that the latter results in differences in azimuthal angles that are not invariant under rotations of the event. We start with an initial azimuthal angle ψ 0 = 0 (the particular choice has no impact on the differences of azimuthal angles that we eventually study). Then, working through the declusterings in the Cambridge/Aachen sequence, for each declustering i we: 1. Compute the normalised cross product, and negative otherwise, andn i−1 is the normalised cross product obtained for the splitting that produced the parent parton.
The variable ψ 12 is now defined as the difference between the ψ obtained for the primary splitting and secondary splittings as selected above (which may not follow in immediate sequence in the C/A declustering). For the case of two successive splittings, this is equivalent to the ψ 12 illustrated in Fig. 3. Likewise we define ψ 11 as the ψ difference between the highest-k t primary passing the z > z cut requirement inside each of the two different jets. The Lund diagrams are illustrated at second order in Figs. 7 and 8.
In practice, in our e + e − implementation of the ψ 12 observable, we cluster the event back to two jets and analyse each jet independently. Our final distributions will be normalised to the number of events, rather than the number of jets.
For the intra-jet definition of the Lund observable ψ 12 , we consider the primary q → qg 1 and secondary splittings g 1 → g 1 g 2 (left) and g 1 → q q (right). We apply a cut (dashed lines), z cut , on the momentum fraction used to identify the primary and secondary split-tings. The red (blue) dots indicate that the azimuthal angles in those splittings are correlated, with the preferred value of ψ 12 being inplane (out-of-plane) where σ tot is the total cross section, Q is the event centre-ofmass energy, θ mn is the opening angle between two emissions m and n and φ (i j)k is the angle between the plane that contains the p i + p j and p k directions and the plane that contains the p i and p j directions. The average is carried out across events, and in each event the sum over each of the i, j and k runs over all particles in the event, 1 . . . N . One may also apply the definition to particles in a single jet of energy E jet , in which case one would replace 8/Q 3 with 1/E 3 jet . Refs. [38,58] provided techniques for resumming such observables, and below we will compare our resummation results for the EEEC to theirs.

MicroJets resummation of spin correlations and comparison to fixed order
Although the above fixed-order analysis gives a sense of the overall structure of spin correlations in quark-and gluoninitiated jets, there are important effects which can only be captured by all-order resummation.
To our knowledge, no analytical result exists for the logarithmic structure of observables sensitive to spin correlations, except for the recently computed all-order result [38] for the 3-point energy correlator, reproduced in Sect. 3.3. In order to enable comparisons to other observables we have implemented a numerical resummation, henceforth called the toy shower, based on the MicroJets code [39,40,59].
Fig. 8 For the inter-jet definition of the Lund observable ψ 11 , to be used in H → gg events, we consider two primary splittings in different hemispheres The toy shower is ordered in an angular-type evolution variable t, where E is the energy of the hard parton initiating the shower, and θ is an angular scale. For a 1-loop running of the strong coupling α s ( p t θ), the scale t is related to the opening angle of the splitting θ by where β 0 = 1 6 (11C A − 4T R n f ), and where we have introduced λ = α s ln θ . Single-logarithmic terms of the form (α s ln(1/θ )) n can then be resummed by solving correspond-ing DGLAP-style equations. Unlike in a real parton shower, there is no kinematic map associated with the emissions, and crucially it thus does not generate any spurious higher-order terms. Since the toy shower is angular ordered, the Collins-Knowles algorithm can readily be applied to it, and we can make predictions for angular observables correct at NLL in the strongly ordered limit. The toy shower can also provide fixed order predictions at O(α 2 s ) in the collinear limit.

Results for Lund declustering observables
We start by evaluating the impact of the single-logarithmic resummation on the size of the spin correlations, a 2 /a 0 , for the observables ψ 12 and ψ 11 introduced above. In order to do so, we consider the distributions of ψ 12 and ψ 11 generated by the toy shower for a value of t max = 0.5 π , which corresponds to λ ≈ −0.3743 for a 1-loop running of the strong coupling, see Eq. (14), which is compatible with the range of α s and L accessible at the LHC (cf. table 1 of Ref. [25]). We then apply the Lund declustering procedure: we identify a splitting as primary or secondary only if it passes the cut z > z cut . We then determine the coefficients a 0 and a 2 from the distribution of the observable. Figure 9 depicts the ratio a 2 /a 0 at fixed second order (FO) and all orders (AO) for the case of ψ 12 in γ * → qq events (quark-initiated intra-jet azimuthal angle), for five values of z cut ∈ {0.05, 0.1, 0.2, 0.3, 0.4} (the fixed order extends to z cut = 0, however the resummation cannot be extended to that region without also addressing soft logarithms). The size of the spin correlations is shown separately for the two different branching histories q → qg(g → gg) in Fig. 9a, and q → qg(g → q q ) in Fig. 9b. The rest channel, characterised by the absence of an intermediate gluon (e.g. q → gq(q → gq) at second order, and other possible histories at all orders), is not shown as it does not produce correlations. 13 The ratio a 2 /a 0 is also given for the sum of all channels in Fig. 9c.
Turning on the IR cut, z cut > 0, first leads to an increase in the absolute size of |a 2 /a 0 |, with the large-z 2 contributions in Fig. 4 driving the spin correlations. As z cut continues to increase, the intermediate gluon becomes harder and spin correlations decrease again up to z cut = 0.5. 14 We observe that the size of the spin correlations agrees between the fixed order and the resummed, all-order case, for small values of z cut ∼ 0, in the separate channels. Once we turn z cut on, the resummation starts diluting the spin correlations, as intermediate, soft gluons can be emitted with values of z < z cut , transporting away some of the spin information before it is propagated to the identified secondary splitting. Thus, the size of the correlations in the resummed prediction decreases with respect to the fixed-order calculation, with the resummed |a 2 /a 0 | being about 88% of the fixed-order value at z cut = 0.4.
If one considers the sum of all channels, given in Fig. 9c, i.e. without the inclusion of any flavour-tagging, the relative size of the resummed correlations starts at about 95% of their fixed-order counterpart at small values of z cut ∼ 0. This decrease finds its source in the higher relative contribution from possible branching histories without correlations -what we call the rest channel -in all-order events: indeed, there is only one such possible history at second order (q → q(q → qg)g), where the emission of the first gluon is hard and the declustering follows the soft quark branch to the secondary splitting). This rest channel contributes to a 0 but not to a 2 , thus the additional damping of |a 2 /a 0 | at small values of z cut in a 2 a 0 all = a 2,q q + a 2,gg a 0,q q + a 0,gg + a 0,rest .
Similar results are shown in Fig. 10 for the gluon-initiated (H → gg) intra-jet observable ψ 12 . While the spin correlations are smaller when compared to the quark-initiated case presented above, we find that they are also less sensitive to resummation: the all-order correlations |a 2 /a 0 | are about O(95%) of their fixed-order counterpart at z cut = 0.4 for the separate channels, and O(92.5%) for the sum of all channels.
Finally, fixed-and all-order predictions are shown in Fig. 11 for the inter-jet observable ψ 11 in H → g 1 g 1 events. The correlations |a 2 /a 0 | are very small in each separate channel, with the exception of the all-quark final state g 1 → qq, g 1 → q q . Because the largest cross section comes from the all-gluon final state, where the correlations |a 2 /a 0 | 1% are extremely small, and because of partial cancellations between different channels, the correlations in the sum of all channels are almost vanishing. Since the effect is so small, there is a large statistical uncertainty on the value of the coefficient a 2 fitted from the Monte-Carlo runs of the toy shower. It still gives results consistent with the observations made above. Additionally, we note that the inter-jet spin correlations are also somewhat more sensitive to resummation than for the intra-jet correlations in ψ 12 studied above.

Validation and results for energy correlators
As part of the validation of the MicroJets code, we also compare it to the analytic resummation of the EEEC presented in Ref. [38]. While in the latter reference, the spin correlations in the EEEC were shown differentially as a function of the opening angle θ S of the secondary (small) splitting and fixed opening angle θ L of the primary (large) splitting, here we show the results integrated over angles to enhance the statistics. Inspired by Ref. [38] we choose the following integration bounds on the opening angles of the primary and secondary branchings 15 : and take α s = 0.0868 corresponding to a hard scale of roughly 1 TeV. The toy shower and analytic resummation results [38] are shown in Fig. 12, summed over all final branching flavour channels, demonstrating good agreement. The figure also shows the second-order expansion, normalised so that its mean value coincides with the resummed result, illustrating a modest reduction in the degree of spin correlations from the resummation, which is similar to our findings above with the Lund declustering observables. 15 Cf. figure 3 in Ref. [38].

Numerical validation of spin correlations within PanScales showers
In this section, we validate the PanScales showers against various numerical predictions. In particular we want to demonstrate that the algorithm described above reproduces fixedorder matrix elements in the strongly ordered limit and that it produces the correct NLL distributions at all orders. Here, we first provide a brief summary of the PanScales showers. A comprehensive description can be found in Ref. [24]. The PanScales showers fit into two categories, namely PanLocal which features local recoil, and PanGlobal which features global recoil. The PanLocal mapping for the emission of p k from a dipole {p i ,p j } is given by where k ⊥ = k ⊥,1 cos(ϕ) + k ⊥,2 sin(ϕ) and −k 2 ⊥ = k 2 t . The coefficients of the map are functions of the ordering variable,  All-order EEEC, λ = −0.4 Fig. 12 All-order comparison of the toy shower and the analytic resummation performed in Ref. [38], for a quark-initiated (left) and a gluon-initiated jet (right). The resummation is performed using α s = 0.0868 and restricting the opening angles as in Eq. (16). The result from a fixed-order expansion, normalised so that its mean coincides with the mean value of the analytic curve, is also shown for comparison v, and an auxiliary variableη, where sĩj = 2 p i · p j , sĩ = 2 p i · Q, and Q is the total event momentum. The quantity β parameterises the choice of ordering variable. The PanGlobal map drops the k ⊥contributions to p i and p j , as well as the terms b i and a j , and instead boosts and rescales the full event. The PanLocal shower comes in a dipole variant, where f = 1 and every soft eikonal is reproduced by the sum of two opposite branching kernels, and in an antenna variant, where and every eikonal is instead contained in a single contribution. We will show results for the PanGlobal shower with β = 0 16 and the PanLocal dipole shower with β = 1/2 and in some places also the PanLocal antenna shower with β = 1/2.

Validation at second order
The first essential test to carry out is at fixed O(α 2 s ) since this is the first order at which spin correlations are non-vanishing. We compare the collinear branching amplitudes presented in Sect. 2.1 against the exact O(α 2 s ) tree-level matrix element, taking the limit of strongly ordered angles for fixed z 1 and z 2 . We then proceed to validate our adaptation of the Collins-Knowles algorithm by applying it at fixed order in the PanScales showers and comparing against the second order expansion of the toy shower and the exact O(α 2 s ) cross section.

Matrix element validation of collinear branching amplitudes
Matrix elements constructed using collinear branching amplitudes are expected to reproduce the exact amplitudes only in the strongly ordered limit. In order to demonstrate that this is the case, we consider the angular correlations ψ 12 and ψ 11 as a function of 12 = max(θ 1 , θ 2 /θ 1 ), 16 We have performed validations also with the PanGlobal β = 1/2 shower but do not show them here to avoid cluttering already busy plots.
where θ i is the opening angle of the ith branching. The strongly ordered limit is then approached as → 0. In practice we achieve this by fixing the energy fractions z 1 and z 2 . We then generate angles θ 1 , θ 2 , φ 1 , and φ 2 from which we can deduce the full kinematics of our event using the PanLocal dipole β = 0.5 map. 17 The full kinematics are then passed to the exact matrix element, obtained using amplitudes from Ref. [60], to be compared directly with the Collins-Knowles weights, Eqs. (1) and (2). In Fig. 13 we show this comparison as a function of ψ 12 and 12 for the processes γ * → qqq q and γ * → qqgg for z 1 = 0.01 and z 2 = 0.3. On the left we show the comparison between the exact matrix element and the shower matrix element without spin correlations, while they are included on the right. We show the absolute difference between the exact matrix element and the shower matrix element normalised to the exact matrix element integrated in a slice of ψ 12 . When spin correlations are not included in the branching amplitudes the exact matrix element and the shower matrix element disagree for all values of 12 and ψ 12 in both processes. The band structure that shows up for ψ 12 = { −3π 4 , −π 4 , π 4 , 3π 4 } corresponds to the points in ψ 12 where the azimuthal modulation of the full matrix element intersects the prediction by the shower, which is flat in the collinear limit. When spin correlations are included we see that the discrepancy between the exact matrix element and the shower matrix element only persists away from the collinear limit and we find that the residual differences vanish as 12 becomes small. The picture is very similar if we instead focus on H → gg processes, as can be seen in Fig. 14. In this case we show the three subprocesses H → gggg, H → ggqq, and H → qqq q as a function of ψ 11 and 11 , using z 1 = z 2 = 0.3. When spin correlations are not included in the shower branching amplitudes, the shower matrix element and the exact matrix element again disagree for all values of 11 and ψ 11 . The same band structure can be seen as in Fig. 13 which is again due to the intersection of the cosine and an almost flat curve. When spin correlations are included in the shower we again find excellent agreement as we approach the strongly ordered limit. Some band structure remains, associated with the existence of azimuths for which the exact and parton shower results happen to coincide exactly.

Tests with full shower phase space and observables
We now consider integrated results over a portion of the Lund plane, i.e. using the analysis approach discussed in −π −π/2 0 π/2 π Δψ 12  Sect. 3.2. This serves to test the full combination of shower phase space generation, branching probabilities and observable implementation (which comes in independent variants for the toy shower and the full shower). We still work at fixed (second) order, and compare the resulting distribution of ψ 12 for several of the PanScales showers, the toy shower, and the exact tree-level matrix element. We produce results for a fixed strong coupling at Q = 100 TeV, setting the following cuts on the primary and secondary splittings: Figure 15 depicts the second-order (normalised) differential cross-section for different values of 12 , where X may be qq or gg. As in the previous section, the latter variable serves as a measure of ordered collinearity for the two successive splittings. We show the distributions obtained by the toy shower, by the PanLocal shower (run with a value of β = 0.5, see Eq. (18)) in its dipole formulation, and by the exact tree-level matrix element. In the bulk of the phase space, where opening angles are not strongly restricted, the distribution of ψ 12 is not expected to agree across all considered setups: the parton shower will feature nonnegligible recoil effects, and the toy shower is applying the strictly-collinear limit of the Collins algorithm irrespective of the opening angles. However, we do expect agreement between those setups, and the exact matrix element, in the limit where the angles are small and strongly ordered.
To quantify the agreement between the different setups in the strongly-ordered collinear limit, lim 12 →0  tions that fulfil the strong angular ordering requirement for a given value of 12 . Specifically, for a distribution with ψ 12 bins running from i = 0, . . . , n − 1 we define the kth DCT coefficient as follows The A 0 and A 2 coefficients can be related to the Fourier coefficients a 0 and a 2 used above, This definition exclusively picks out the even (cosine-like) Fourier components, whereas the odd (sine-like) modes are zero by definition. The results are summarised in Fig. 16, for γ * → qqq q (left column) and γ * → qqgg (right column). We show the first four DCT coefficients, normalised by the toy shower A ts 0 ( 12 ), for two PanScales showers, the exact matrix element, and the toy shower. For the latter, since the toy shower is free of recoil effects, the correlations introduced by the spin algorithm take the form given by Eq. (11) irrespective of the value of 12 . Thus, for the toy shower, A 0 and A 2 are the only non-zero coefficients, and 2 A 2 /A 0 = a 2 /a 0 gives the relative size of the integrated spin correlations, which are of the order O(−77%) for the g → q q channel, and O(+7%) for the g → gg channel in most of the phase space. 18 The parton showers and the matrix element partly generate non-zero coefficients at large values of 12 , all of which become consistent with zero in the stronglyordered limit at ln 12 ∼ −4, with the exception of the coefficient A 2 encoding the spin correlations. All setups produce compatible values of A 2 at small 12 . We also observe that the PanScales showers display distortions that are of same sign, and similar in size, as the exact matrix element at larger values of 12 . The difference between the PanLocal dipole and the PanGlobal showers is a consequence of the different kinematic map rather than the difference in β values that we have used (0.5 and 0.0 respectively).

Validation at all orders in α s
We now turn to the validation of our spin-correlation algorithm in the PanScales parton showers. To test specifically the single logarithmic (NLL), α n s L n , terms generated by the shower, we run the PanGlobal (β = 0) shower and the dipole and antenna versions of the PanLocal (β = 0.5) showers with asymptotically small α s → 0, keeping λ = α s L fixed. For the Lund declustering observable (cf. Sect. 3.2), e −|L| is the smallest value that we allow for k t,2 /Q, or equivalently at α n s L n accuracy, the smallest value for θ 2 , while θ 1 and k t,1 are allowed to take on any value (typically, the Lund declustering procedure ensures θ 1 > θ 2 and k t,1 > k t,2 ). For the EEEC variable (cf. Sect. 3.3), e −|L| is the smallest value that we allow for θ S , while θ L is allowed to take on any value larger than θ S . For the purpose of the following tests we choose a value of λ = −0.5, which corresponds roughly to the range of α s and L accessible at the LHC.
For results generated by the toy shower, this translates to setting the toy shower's cutoff scale to t max = t (λ) as given in Eq. (14). We use a 1-loop running of α s for all results shown below (in our α s → 0 but fixed λ = α s L limit, higher-loop running leaves the results unchanged).
In the following, we consider both γ * → qq and H → gg initial setups. A variety of issues arise from the use of small values of α s and large values of L. The approaches that we take to dealing with them include truncating the shower to avoid large numbers of soft gluons, tracking differences in angles between particles as well as their 4-momenta, the use of a numerical type that extends the available exponent beyond that normally accessible in double precision, and a stand-in analysis using the shower tree structure rather than the full Lund declustering analysis. Some of the techniques were introduced in Refs. [24,25] and Appendix D gives further details about the techniques that are new here, including a number of tests to verify that our conclusions remain robust even with these techniques. We apply a cut z cut = 0.1 in the identification of the primary and secondary (respectively, both primary) splittings used in reconstructing the Lund observable ψ 12 (respectively, ψ 11 ). We also use a stand-in analysis for the EEEC observable, presented in Appendix D.2, which again is valid in the extreme collinear limit. Figures 17 and 18 depict the distribution of our new Lund observables, ψ 12 and ψ 11 , as well as the EEEC, for γ * → qq, respectively for H → gg initial events. The three showers PanGlobal (β = 0), PanLocal (β = 0.5, dipole and antenna versions) are compared to the numerically resummed result obtained from the toy shower. In all cases, we show the contributions stemming from the different channels to the full observable. The relative deviation between the PanScales showers and the toy shower is shown on the right, separately for each channel, and is compatible with zero with statistical uncertainties below the 5 permille level.

Phenomenological remarks
We comment on three aspects here that are potentially relevant for phenomenological applications.
Our first comment concerns the relative size of spin correlations in the EEEC and the ψ 12 Lund declustering observable. The EEEC has the advantage of not requiring a z cut , reducing the number of parameters that need to be chosen for the observable. However its weighting with the energies in Eq. (12) tends to favour configurations where a q → qg(g → x y) splitting shares energy equally between the three final particles. In the notation of Figs. 4 and 5, this corresponds to z 1 2/3 and z 2 1/2. While z 2 1/2 acts to enhance the spin correlations, z 1 2/3 tends to reduce them. In contrast, with the Lund declustering ψ 12 one can adjust the cuts on the z 1 and z 2 values so as to max- Exact ME γ * → qq → qqg 1 g 2 (b) qqg 1 g 2 Fig. 15 The second-order distribution of ψ 12 from the toy shower, the dipole PanLocal shower, and the tree-level matrix element, shown for increasingly collinear configurations, ln 12 ∈ {−1.3, −1.7, −2.1, −2.5, −2.9, −3.3}. Results are given for a quark-only and b two quarks-two gluons final states imise the azimuthal modulations. 19 Table 3 summarises the degree of azimuthal modulation for different observables in γ * → qq events. With our default (non-optimised) cuts of z 1 and z 2 > 0.1, we see substantially larger azimuthal modulations than in the EEEC variables, both in individual flavour channels and in their sum. The potential for further enhancement of the modulations is made evident by the results obtained with the z 2 > 0.3 requirement.
Our second comment concerns the sum over all flavour channels. The results shown here have been obtained with n f = 5 light flavours. The final magnitude of the spin correlations after the sum over flavour channels is quite sensitive to the cancellation between g → qq and g → gg splittings and the degree of cancellation is strongly influenced by the value of n f . At the scales where one might aim to probe spin correlations, the c-and especially b-quark masses are not entirely negligible. A full phenomenological study of the flavoursummed structure of azimuthal correlations might, therefore, needs to take into account finite quark-mass effects. Note that effects related to k t values in the neighbourhood of a heavy-quark threshold are formally suppressed by a logarithm. For a complete understanding of phenomenological expectations one would also want to examine the impact of other subleading logarithmic effects, as well as contributions suppressed by powers of k t /Q, and possibly also non-perturbative corrections. It would clearly also be of interest to find ways of carrying out measurements with flavour tagging, given the strong effects to be seen with g → qq splittings. While b and c flavour tagging are the most obviously robust starting points in this respect, one may also wish to consider s tagging [61] and generic light-quark versus gluon discrimination observables.
Our final comment concerns the inter-jet ψ 11 observable. If Higgs boson decay to two gluons is eventually observed and if one can carry out flavour tagging so as to have visible net modulation, this observable could provide an interesting example of an EPR type measurement at colliders (a constraint on a gluon's splitting angle would then translate to a constraint on the distance travelled by the gluon before it splits). and their associated coefficients suitably converge to identical values in the strongly angular-ordered limit 12 → 0

Conclusions
The developments shown in this paper are among the last remaining aspects needed in order to claim complete NLL accuracy in the sense of Ref. [24] for the PanScales family of parton showers. The one component that now remains (at leading colour) is a treatment of the spin correlations for soft emissions, an aspect that we leave to future work. The approach we have taken is largely based on that proposed by Collins and Knowles, and in part also on the extension to dipole showers proposed by Richardson and Webster. The main difference relative to the Richardson and Webster work is our use of spinor products, avoiding the need for correlating branching variables across different steps and associated boosts. The spin correlation treatment is then expected to be correct at single-logarithmic accuracy for showers that adhere to the general requirements set out in Ref. [24].
Relative to earlier work on spin correlations in parton showers, one of the main novelties of this paper is the framework for validating the implementation of spin correlations. As part of our testing framework we extended the MicroJets code to enable single-logarithmic resummation for a range of spin-correlation observables. Figures 13, 14, 15 and 16 show Table 3 The relative magnitude of the azimuthal modulation, a 2 /a 0 (cf. Eq. (11)), for the EEEC and Lund intra-jet ψ 12 observables, the latter for two sets of cuts on z 1 and z 2 . The results are shown for γ * → qq events for n f = 5, separately for two specific flavour channels, as well as the sum over all flavour channels (including the channel without spin correlations, q → qg). As in Fig. 17, the results are obtained in the limit α s → 0 for fixed λ = α s L = −0.5 and for the Lund declustering ψ 12 we consider events with k t,2 /Q > e −|L| , while for the EEEC ψ we consider events with θ S > e −|L| Flavour channel for 2nd splitting g → qq g→ gg All EEEC −0.36 0.026 −0.008  12)). The results are obtained in the limit α s → 0 for fixed λ = α s L = −0.5. For the Lund declustering ψ 12 we consider events with k t,2 /Q > e −|L| and for the EEEC ψ we consider events with θ S > e −|L| is more differential in momentum fractions of the different branchings. This makes it possible to enhance the relative magnitude of the spin-correlation signal, cf. Table 3. Spincorrelation effects can be large for g → qq splittings, while they are smaller, with an opposite sign for g → gg splittings (cf. Figs. 4, 5 and 6). All-order resummation has a modest effect on them (cf. Figs. 9, 10 and 11). Phenomenologically, the opposite signs for g → gg and g → qq lead to a par-tial cancellation of the spin correlation effects in flavourinsensitive observables, although our studies indicate that some effect remains visible. Further study of the potential for measuring these effects, possibly in combination with flavour tagging, would clearly be of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

A Deriving the branching amplitudes
In this appendix we collect some details regarding the computation of branching amplitudes in terms of spinor products and their numerical evaluation. The calculations here are done following the conventions of Ref. [62]. For two light-like momenta p and q, we define the spinor product where λ = ±1 is the Dirac spinor helicity. Spinor products have the properties Following the HELAS conventions [63], the polarisation vector of a gluon with momentum p is defined as * μ where r is a light-like vector specifying the gauge. For calculational purposes, the Chisholm identity [62] is often useful. We define the momenta of a collinear splitting as p a → p b + p c , such that p b = zp a and p c = (1 − z) p a . In this limit, the gauge vector r cancels from the collinear branching amplitudes, and they may be written in terms of a single spinor product S λ ( p b , p c ) using the identities Next, we compute the branching amplitudes for all colourstripped collinear QCD branchings, dropping any overall phase factors for convenience.

q → q g
The branching amplitude is The non-zero helicity configurations are The branching amplitude is The non-zero helicity configurations are For this case, we first derive some general properties. If we consider the collinear limit and pick all gluons to have the same gauge vector r , we find * where p i , p j ∈ {p a , p b , p c }. The branching amplitude is

B An example of calculating the functions A(z), B(z)
We illustrate the computation of the spin-correlated matrix element squared, taking the example of the observable ψ 12 at second order for a quark-initiated jet, q → qg(g → q q ), of Fig. 1. We then cast the result into the form of Eq. (11). The matrix element squared, summing over all helicities, is given by where the last line includes the flipping of all helicities, and we have already excluded forbidden helicities. Using Eq. (5) and inserting the functions F λ a λ b λ c a→bc (z) given in Table 1, the above reduces to where each of the four terms in the parenthesis of Eq. (41a) comes from one of the four lines in Eq. (40). Finally, Eq. (41b) makes the form of A(z 1 ) and B(z 2 ) explicit.

C Rotational invariance in the collinear limit
As part of the numerical evaluation of spinor products in Eq. (39), a set of reference vectors must be selected. Any dependence on these reference vectors vanishes in the calculation of a full squared matrix element, but the Collins-Knowles algorithm only accounts for the contributions that dominate the matrix element in the collinear limit. Consequently, any dependence on the reference vectors should also vanish from the Collins-Knowles result as long as the shower emissions are strongly ordered and sufficiently collinear.
In this appendix, we validate this property in the Pan-Scales implementation of the Collins-Knowles algorithm. To that end, like in Sect. 4.1.1, the Collins-Knowles weight as a function of an angular resolution is considered. Rather than comparing to a matrix element directly, the weight is computed for an event, and for the same event that is rotated randomly. As the reference vectors stay fixed, rotating the event is equivalent to rotating the reference vectors in the original event. Figure 19 shows the relative difference of these weights for two shower emissions off γ * → qq and H → gg hard scatterings. The relative difference in weights decreases as a power of the angular resolution, and is thus a power correction that vanishes in the collinear limit. Numerically, these effects are also small compared to the size of the spin correlations themselves.

D Technical details of the all-order comparisons
To facilitate the isolation of the NLL structure of the parton shower, the all-order comparisons of Sect. 4.2 are performed at extremely small values of α s and extremely large values of the logarithm. Several techniques need to be employed to maintain numerically feasible analyses in these extreme regimes. In this appendix we detail these techniques.

D.1 Removal of soft radiation
The particle multiplicity generated by the parton shower scales like α s L 2 . This means that, because the product α s L is kept constant as α s decreases and L increases, the multiplicity also increases. At the logarithmic values used in Sect. 4.2, the multiplicity has increased to levels that cause the event generation to become numerically unfeasible. However, the all-order observables considered here are insensitive to soft gluon radiation, and the spin correlations incorporated by the Collins-Knowles algorithm are also unaffected. As a result, as long as any radiation that is removed is sufficiently soft that its recoil has no noticeable impact on the momenta of the partons that dominate the Lund declustering and energy-correlator observables, the all-order tests should remain unaffected by the removal of soft gluon radiation. Such a cut is implemented in the PanScales shower by limiting the sampling range of the auxiliary parton shower variable that controls the collinear momentum fraction. This strategy is more efficient than the alternative of vetoing soft branchings, as the shower would spend much of its time generating and vetoing soft emissions at low scales. An illustration of this procedure is shown in Fig. 20.
To validate the legitimacy of the application of this cut on soft radiation, we compare the predictions of the PanScales showers and our implementation of Dire v1 for ψ 12 with z cut = 0.1 at α s = 0.01 and L = −27.5. Figure 21 shows the difference between the distributions with and without the application of the collinear cut ln z > ln z PS cut = −10. This value of z PS cut is significantly below the z values that dominate in our observables. For the PanScales showers, we see that the cut has no statistically significant effect on the azimuthal distribution and thus we conclude that we can safely use a z PS cut in our logarithmic accuracy tests. For other showers, the cut can have an impact on the results, as illustrated with the Dire v1 dipole shower, where the cut induces a change in a 2 that is a significant fraction of the actual a 2 . We attribute this to transverse recoil effects between gluons of commen-surate k t values, of the kind discussed in Refs. [23,24]. For such showers, a test of the logarithmic structure of the spin correlations would need to be carried out without a z PS cut . We suspect that there are observables for which such tests would reveal problems in the α n s L n logarithmic structure, associated with the following type of configuration: consider a hardcollinear g → qq splitting with a transverse momentum k t,1 , followed by a soft-collinear emission with transverse momentum k t,2 k t,1 from the dipole that contains the q. Within standard dipole-shower recoil schemes, the softcollinear emission can take its recoil from the quark, altering its 2-dimensional vector transverse momentum, effectively smearing out the azimuthal angle of the qq pair. 20 In certain circumstances (for example if one considers events where the g → qq splitting is the highest-k t splitting in the event), the product of squared-matrix element and phase space where the recoil can occur will lead to a factor α s L for this smearing to occur, thus spoiling single-logarithmic accuracy. Note, however, that in the observables that we actually consider, we do not require either of the collinear splittings to be the hardest in the event, and we expect that feature to bring further complications, of the kind discussed in section 3 of the supplemental material of Ref. [24], whereby additional double-logarithmic effects would be expected to (at least) partially mask these issues in the all-order tests, potentially associated with superleading logarithms, α n s L m with m > n, for observables that should only have m ≤ n. Based on these arguments, we consider the presence of any effect from emissions below z cut , as seen in Fig. 21, to be a sign of potential danger.

D.2 Stand-in observables
While performing tests at small α s and large logarithms, it quickly becomes challenging to accurately evaluate the tiny angles between collinear momenta. To avoid large numerical cancellations in the evaluation of inner products, the PanScales showers have the option to explicitly track the difference in direction between dipole components. This technique, as well as the use of a new floating-point type (double_exp) that increases the maximum size of the exponent of a regular double-precision type, was already used in [25], and was equally crucial for the all-order tests performed here. While it allows the shower to be run at asymptotically small values of the coupling, the evaluation of the observables is not as straightforward.
ln k t η Fig. 20 An example of the shower Lund plane after the first branching, with a cut on soft emissions. The red shaded area indicates the part of phase space removed by the soft-emission cut, which is applied relative to the total event energy. The black dashed line represents a possible observable cut on z, such as that applied by the Lund declustering observables. It is applied relative to the parent momentum. To avoid removing parts of the phase space above the dashed line due to recoil effects, the soft-emission cut is applied well below the observable cut The key also shows the result for a 2 /a 2 where a 2 is the difference between the a 2 Fourier coefficient as obtained with and without the cut. The result labelled "Dipole" corresponds to our implementation of the Dire v1 [64] shower algorithm supplemented with our spin-correlation algorithm and we expect it to be representative of a variety of standard dipole-type showers. It indicates that the use of a ln z cut would not be safe for testing spin correlations in such standard dipole-type showers In the case of Lund-declustering observables, a regular analysis is technically limited by the absence of directional difference information for parton pairs that are not in the same dipoles. In events with at most a Born qq pair, we have implemented functionality whereby we use the full dipole chain and in-dipole direction differences to construct a lookup table for the direction differences between any pair of partons. For N particles this has an O N 2 time and memory cost, while normal e + e − clustering in FastJet also carries an O N 2 time cost, but an O (N ) memory cost. Used together with a winner-takes-all recombination scheme [65][66][67] in a version of the FastJet [68] code adapted to work with the double_exp type, we can carry out the Lund declustering analysis at asymptotic values of the logarithm. However the current technical limitation of having at most the Born qq pair in the event means that we need an alternative strategy for complete events.
Accordingly, we have developed a stand-in analysis that uses the binary tree generated by the Collins-Knowles algorithm, even though it is normally unobservable. During the showering, the Lund structure is determined using explicit shower information, where at every collinear branching that appears in the binary tree, the appropriate azimuthal information is stored if the collinear momentum fraction cuts are passed. This procedure is equivalent to the evaluation of the exact observable if the jet clustering produces the same binary tree (for the hard emissions that we ultimately consider in the observable) as in our implementation of the Collins-Knowles algorithm and if recoil effects that can affect the momentum fraction cuts are absent. At asymptotically small α s , for the subset of PanScales showers that passed the NLL tests in Ref. [24], and with the condition that one works with a hardness cut z > z cut such that α s ln z cut 1, these are both valid assumptions. We have verified this by considering events without g → qq splittings and establishing that the full Lund declustering and the stand-in analysis produced event-by-event identical results. Note that the stand-in analysis also allows easy access to the flavour separation of the spin-correlation effects, without the need to consider flavoured jet clustering.
In the case of the EEEC, a major issue it its O(N 3 ) time complexity, which quickly becomes prohibitive at the multiplicities under consideration. The following procedure, which again makes direct use of the binary tree, is again equivalent to the full observable in the asymptotic limit. For all nodes a n that are not terminal, and where n is the depth of that node, do as follows: 1. Retrieve the opening angle θ S and the normal vector n S of the branching of node a n . 2. Retrieve the global momentum fractions z 1 = 2E 1 /Q, z 2 = 2E 2 /Q of the outgoing momenta 1 and 2 of node a n . 3. Iteratively for i ∈ {1, . . . , n − 1}: • Retrieve the opening angle θ L and the normal vector n L of the branching of node a i . • Retrieve the global momentum fraction z 3 = E 3 /2E CM , where 3 is the outgoing momentum of a i that is not a i+1 . • Compute the signed angle ψ between n S and n L .
• Add to the ψ-θ S -θ L -bin with weight 2z 1 z 2 z 3 . This procedure has time-complexity O(N log N ) and is again equivalent to the exact observable if θ S θ L and in the absence of recoil effects, which are both the case for asymptotically small α s for the showers we study. The above algorithm was validated by verifying that it yields identical results to the exact observable evaluation for a representative number of events.