What is a parton shower?

TL;DR

This paper explores the theoretical structure of idealized parton shower generators in QCD, proposing a formulation that relates shower operators to infrared singularities and extends to any perturbative order.

Contribution

It introduces a formulation for parton showers at any order in _s, linking shower operators to QCD infrared singularities, and clarifies the theoretical foundation of parton showers.

Findings

Shows how shower operators relate to QCD infrared singularities

Proposes a formalism for order _s^k showers

Highlights the gap between theoretical formulation and practical algorithms

Abstract

We consider idealized parton shower event generators that treat parton spin and color exactly, leaving aside the choice of practical approximations for spin and color. We investigate how the structure of such a parton shower generator is related to the structure of QCD. We argue that a parton shower with splitting functions proportional to $\alpha_s$ can be viewed not just as a model, but as the lowest order approximation to a shower that is defined at any perturbative order. To support this argument, we present a formulation for a parton shower at order $\alpha_s ^k$ for any $k$. Since some of the input functions needed are specified by their properties but not calculated, this formulation does not provide a useful recipe for an order $\alpha_s ^k$ parton shower algorithm. However, in this formulation we see how the operators that generate the shower are related to operators that specify the infrared singularities of QCD.

Tables

Table 1. Table 1: New and old notation

New	Old
${𝒮}_{\mathrm{pert}}^{(1,0)}$	${\mathscr{H}}_I^{\mathrm{pert}}$
${𝒮}_{\mathrm{pert}}^{(0,1)}$	$-{𝒮}^{\mathrm{pert}}$
$\left[ℱ\circ {\overline{𝒮}}^{(1,0)}\right]{ℱ}^{-1}$	$𝒱$
${𝒮}^{\left(1\right)}-{𝒮}_{i\pi }^{(0,1)}$	${\mathscr{H}}_I-𝒱$
${𝒮}_{𝒱}^{\left(1\right)}+{𝒮}_{i\pi }^{(0,1)}$	$𝒱-𝒮$

Equations

$$\begin{split}\sigma[J]={}&\big{(}{1}\big{|}\mathcal{O}_{J}\,\mathcal{U}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})\,\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})\,\mathcal{F}(\mu_{\textsc{h}}^{2})\big{|}{\rho_{\textsc{h}}}\big{)}+\mathcal{O}(\alpha_{\mathrm{s}}^{k+B+1})+\mathcal{O}(\mu_{\mathrm{f}}^{2}/Q[J]^{2})\;.\end{split}$$

$${\mathcal{U}}_{\mathcal{V}}({\mu }_{\mathrm{f}}^2,{\mu }_{\text{textsc}h}^2)=\mathbb{T}\exp \left({\int }_{{\mu }_{\mathrm{f}}^2}^{{\mu }_{\text{textsc}h}^2}\frac{d{\mu }^2}{{\mu }^2}{\mathcal{S}}_{\mathcal{V}}({\mu }^2)\right),$$

$$\mathcal{U}({\mu }_{\mathrm{f}}^2,{\mu }_{\text{textsc}h}^2)=\mathbb{T}\exp \left({\int }_{{\mu }_{\mathrm{f}}^2}^{{\mu }_{\text{textsc}h}^2}\frac{d{\mu }^2}{{\mu }^2}\mathcal{S}({\mu }^2)\right),$$

$$\begin{array}{rl}{P}_{\mathrm{A}}=& \left(P_{\mathrm{A}}^{+},0,0\right),\\ P_{\mathrm{B}}=& \left(0,P_{\mathrm{B}}^{-},0\right).\end{array}$$

$$p_{\text{textsc}h}=(e^{y_{\text{textsc}h}}\sqrt{(m_{\text{textsc}h}^2+{\mathbf{p}}_{\text{textsc}h,\perp }^2)/2},e^{-y_{\text{textsc}h}}\sqrt{(m_{\text{textsc}h}^2+{\mathbf{p}}_{\text{textsc}h,\perp }^2)/2},{\mathbf{p}}_{\text{textsc}h,\perp }).$$

$$p_i=(e^{y_i}\sqrt{{\mathbf{p}}_{i,\perp }^2/2},e^{-y_i}\sqrt{{\mathbf{p}}_{i,\perp }^2/2},{\mathbf{p}}_{i,\perp }).$$

$$\sigma [J]\approx \sum _{a,b}\int d{\eta }_{\mathrm{a}}\int d{\eta }_{\mathrm{b}}{f}_{a/A}({\eta }_{\mathrm{a}},{\mu }^2)f_{b/B}({\eta }_{\mathrm{b}},{\mu }^2)\hat{\sigma }[J].$$

$$\begin{array}{rl}\sigma [J]=& \int d{y}_{\text{textsc}h}\frac{d{\sigma }_0}{d{y}_{\text{textsc}h}}{J}_0(p_{\text{textsc}h})+\int d{y}_{\text{textsc}h}d{y}_{1}d{\mathbf{p}}_{1,\perp }\frac{d{\sigma }_1}{d{y}_{\text{textsc}h}d{y}_{1}d{\mathbf{p}}_{1,\perp }}{J}_1(p_{\text{textsc}h},p_1)\\ & +\frac{1}{2!}\int d{y}_{\text{textsc}h}d{y}_{1}d{y}_{2}d{\mathbf{p}}_{1,\perp }d{\mathbf{p}}_{2,\perp }\frac{d{\sigma }_2}{d{y}_{\text{textsc}h}d{y}_{1}d{y}_{2}d{\mathbf{p}}_{1,\perp }d{\mathbf{p}}_{2,\perp }}{J}_2(p_{\text{textsc}h},p_1,p_2)\\ & +\cdots .\end{array}$$

$$\{p{\}}_m=\{p_{\text{textsc}h},p_1,\dots ,p_m\}.$$

$$\begin{array}{rl}& p_{m+1}\to z{\tilde{p}}_m,\\ & p_m\to (1-z){\tilde{p}}_m.\end{array}$$

$$J_{m+1}(\{p_{\text{textsc}h},p_1,\dots ,p_{m-1},p_m,p_{m+1}\})\to J_m(\{p_{\text{textsc}h},p_1,\dots ,p_{m-1},{\tilde{p}}_m\})$$

$$p_{m+1}\to \xi P_{\mathrm{A}}$$

$$p_{m+1}\to \xi P_{\mathrm{B}}.$$

$$J_{m+1}(\{p_{\text{textsc}h},p_1,\dots ,p_m,p_{m+1}\})\to J_m(\{p_{\text{textsc}h},p_1,\dots ,p_m\})$$

$${\mu }_{\mathrm{split}}^2=\frac{\sqrt{Q_{\text{textsc}H}^2}}{E_m+E_{m+1}}(p_m+p_{m+1}{)}^2.$$

$${\mu }_{\mathrm{split}}^2=-\frac{\sqrt{Q_{\text{textsc}H}^2}}{E_{\mathrm{a}}-E_{m+1}}(p_{\mathrm{a}}-p_{m+1}{)}^2.$$

$$⟨\frac{\delta \sigma [J]}{{\mu }_{\mathrm{split}}^2}⟩=\frac{\sigma [J]}{Q^2[J]}.$$

$$Q^2[J]\approx \frac{\sqrt{Q_{\text{textsc}H}^2}}{P_{\mathrm{J}}}\frac{(R{P}_{\mathrm{J}}{)}^2}{N}.$$

$$N=\left|\frac{P_{\mathrm{J}}}{\sigma (P_{\mathrm{J}})}\frac{d\sigma (P_{\mathrm{J}})}{d{P}_{\mathrm{J}}}\right|.$$

$$\delta \sigma [J]\sim \frac{{\mu }_{\mathrm{split}}^2}{Q^2[J]}\sigma [J].$$

$$\delta \sigma [J]\sim \frac{(1\mathrm{GeV}{)}^2}{Q^2[J]}\sigma [J].$$

$$\sigma [J]=\sum _{a,b}\int d{\eta }_{\mathrm{a}}\int d{\eta }_{\mathrm{b}}{f}_{a/A}({\eta }_{\mathrm{a}},{\mu }^2)f_{b/B}({\eta }_{\mathrm{b}},{\mu }^2){\hat{\sigma }}_{a,b}[{\eta }_{\mathrm{a}},{\eta }_{\mathrm{b}},{\mu }^2;J]+\mathcal{O}\left({\left[m/Q\right]}^n\right).$$

$$\begin{array}{rl}{\hat{\sigma }}_{a,b}[{\eta }_{\mathrm{a}},{\eta }_{\mathrm{b}},{\mu }^2;J]=& {\hat{\sigma }}_{a,b}^{(0)}[{\eta }_{\mathrm{a}},{\eta }_{\mathrm{b}},{\mu }^2;J]+\left[\frac{{\alpha }_{\mathrm{s}}({\mu }^2)}{2\pi }\right]{\hat{\sigma }}_{a,b}^{(1)}[{\eta }_{\mathrm{a}},{\eta }_{\mathrm{b}},{\mu }^2;J]\\ & +{\left[\frac{{\alpha }_{\mathrm{s}}({\mu }^2)}{2\pi }\right]}^2{\hat{\sigma }}_{a,b}^{(2)}[{\eta }_{\mathrm{a}},{\eta }_{\mathrm{b}},{\mu }^2;J]+\cdots .\end{array}$$

$$\{p,f{\}}_m=\{{\eta }_{\mathrm{a}},a,{\eta }_{\mathrm{b}},b,p_{\text{textsc}h},p_1,f_1,\cdots ,p_m,f_m\}.$$

$$\rho(\{p,f\}_{m})=\sum_{\{s,s^{\prime},c,c^{\prime}\}_{m}}\rho(\{p,f,s,s^{\prime},c,c^{\prime}\}_{m})\big{|}{\{s,c\}_{m}}\big{\rangle}\big{\langle}{\{s^{\prime},c^{\prime}\}_{m}}\big{|}\;.$$

$$dP\times\langle O\rangle=[d\{p\}_{m}]\sum_{\{s,s^{\prime},c,c^{\prime}\}_{m}}\rho(\{p,f,s,s^{\prime},c,c^{\prime}\}_{m})\,\big{\langle}{\{s^{\prime},c^{\prime}\}_{m}}\big{|}O\big{|}{\{s,c\}_{m}}\big{\rangle}\;,$$

$$\begin{split}\big{[}d\{p\}_{m}\big{]}\equiv{}&\frac{d^{d}p_{\textsc{h}}}{(2\pi)^{d}}\,2\pi\delta_{+}(p_{\textsc{h}}^{2}-m^{2}_{\textsc{h}})\prod_{i=1}^{m}\left\{\frac{d^{d}p_{i}}{(2\pi)^{d}}\,2\pi\delta_{+}(p_{i}^{2})\right\}d\eta_{\mathrm{a}}\,d\eta_{\mathrm{b}}\\ &\times(2\pi)^{d}\delta\bigg{(}p_{\mathrm{a}}+p_{\mathrm{b}}-p_{\textsc{h}}-\sum_{i=1}^{m}p_{i}\bigg{)}\ \;.\end{split}$$

$$\big{(}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{|}{\rho}\big{)}=\rho(\{p,f,s,s^{\prime},c,c^{\prime}\}_{m})\;.$$

$$1=\sum_{m}\frac{1}{m!}\int\![d\{p\}_{m}]\sum_{\{f\}_{m}}\sum_{\{s,s^{\prime},c,c^{\prime}\}_{m}}\big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)}\big{(}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{|}\;.$$

$$\big{(}{1}\big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)}=\big{\langle}{\{s^{\prime},c^{\prime}\}_{m}}\big{|}{\{s,c\}_{m}}\big{\rangle}\;.$$

Full text

perturbative QCD, parton shower

††preprint: DESY 17-063

I Introduction

Parton shower event generators for hadron collisions, such as Herwig Herwig , Pythia Pythia , and Sherpa Sherpa , perform calculations of cross sections according to an approximation to the standard model or its possible extensions. They are essential for the analysis of experiments at the Large Hadron Collider. The main ideas behind these generators were developed in the 1980s EarlyPythia ; Gottschalk ; EarlyHerwig . There has been extensive development of the algorithms since then dipolesG ; lambdaR ; dipolesGP ; Herwig1992 ; Pythia1994 ; SherpaAlpha ; SjostrandSkands ; NSI ; NSRingberg ; Geneva . The successor programs Herwig ; Pythia ; Sherpa ; DinsdaleTernickWeinzierl ; SchumannKrauss ; PlatzerGiesekeI ; PlatzerGiesekeII ; Deductor ; HochePrestel ; Vincia2016 , are quite sophisticated. Useful reviews of the field can be found in review2011 ; SjostrandReview . One of the available successor programs is our own, Deductor NSI ; NSII ; NSspin ; NScolor ; Deductor ; ShowerTime ; PartonDistFctns ; ColorEffects ; NSThreshold ; NSThresholdII . This paper concerns the perturbative part of these parton shower generators, leaving aside models for the underlying event and hadronization. Furthermore, we consider an idealized version of a parton shower generator in which one accounts exactly for spin and color. Approximations for spin and color are a separate issue, which we do not discuss here.

Our aim in this paper is to investigate how the structure of a shower that treats spin and color exactly is related to the structure of QCD. In particular, we ask whether a shower with splitting functions proportional to $\alpha_{\mathrm{s}}$ can be the leading order approximation to something that is defined at any order in $\alpha_{\mathrm{s}}$. We find an affirmative answer to this question. Specifically, we find that there is a construction for defining a parton shower that generalizes current showers at any order of perturbation theory. We find also that the problem of relating the structure of a parton shower to the structure of QCD is not as straightforward as one might have naively guessed. First, the construction makes use of functions analogous to the Catani-Seymour dipole splitting functions CataniSeymour that specify the infrared behavior of QCD, but beyond leading order the shower splitting functions are not related to the functions that specify the infrared behavior of QCD by anything so simple as just changing their sign. Second, the formulas for the shower automatically includes factors that sum threshold logarithms Sterman1987 ; AppellStermanMackenzie ; Catani32 ; CataniWebberMarchesini ; Magnea1990 ; Korchemsky1993 ; CataniManganoNason32 ; SudakovFactorization ; KidonakisSterman ; KidonakisOderdaSterman ; KidonakisOderdaStermanColor ; LaenenOderdaSterman ; CataniManganoNason ; LiJointR ; KidonakisOwensPhotons ; StermanVogelsang2001 ; KidonakisOwensJets ; LSVJointR ; KSVJointR ; Kidonakis2005 ; MukherjeeVogelsang ; Bolzoni2006 ; deFlorianVogelsang ; RavindranSmithvanNeerven ; RavindranSmith ; BonviniForteRidolfi ; BFRSaddlePt ; deFlorianJets ; Catani:2014uta ; Ahmed:2014uya ; BonviniMarzani ; Bonvini:2015ira ; ManoharSCET ; IdilbiJiSCET ; BecherNeubert ; BecherNeubertPecjak ; BecherNeubertXu ; BecherSchwartz ; Stewart:2009yx ; Beneke:2010da ; Bauer:2011uc ; Wang:2014mqt ; Lustermans:2016nvk ; Bonvini:2012az ; StermanZeng ; Bonvinietalcompare . These factors are not included in current parton shower generators at even leading order, except for Deductor NSThreshold ; NSThresholdII . Third, the formulas automatically include matching of the parton shower to a perturbative calculation beyond leading order of the hard scattering that starts the shower. This is fairly straightforward MCatNLO ; PowHeg1 ; PowHeg2 ; MENLOOPS ; PowHegBox ; FxFx ; MINLO ; MINLO2 ; MEPSatNLO ; UNLOPS ; PlatzerMatching ; NNLOPS ; MG5aMCNLO ; GENEVANNLO ; Czakon ; KrkNLO2015 ; KrkNLO2016 ; MINLO3 ; HocheMatching1 ; HocheMatching2 ; CKKW ; CKKWL ; MLM for a leading order shower, but not for a shower beyond leading order.

We believe that is important to understand that a lowest order parton shower generator can represent the lowest order in a systematically improvable approximation. However, the construction in this paper is not a useful recipe for actually creating a parton shower algorithm beyond the leading order: some of the components of the recipe are specified by their properties but not explicitly constructed. A complete construction will require specifying such choices as a shower ordering variable and momentum mappings. Such a specification will require considerable effort, which lies beyond the scope of this paper.

The papers JadachJHEP ; JadachPolonica ; SkandsNLO ; HocheNLO1 ; HocheNLO2 present treatments of a parton shower at order $\alpha_{\mathrm{s}}^{2}$. Ref. SkandsNLO attempts to extend the dipole splittings often used in a leading order shower to a higher order analogue for the case of $e^{+}e^{-}$ annihilation. This approach is similar in spirit to what we do in this paper. An alternative approach JadachJHEP ; JadachPolonica ; HocheNLO1 ; HocheNLO2 , concentrates on the next-to-leading order (NLO) Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) kernel for the evolution of parton distribution functions. Although the parton evolution kernels play a role in our formalism, it is not a central role.

II Overview

The construction of a parton shower at any perturbative order, presented in Sec. VI, is rather abstract. In this section, we attempt to provide an overview of what the later mathematics is intending to do.

In our view, it is most useful to think of a parton shower algorithm as beginning with the theorem factorization that allows us to write a cross section for an infrared safe observable as a convolution of a hard scattering factor with parton distribution functions. Then the parton shower fills in more detail by using the renormalization group. The parton shower develops with decreasing values of a parameter that is a measure of the hardness of interactions.111Herwig then rearranges the ordering of splittings in its shower so that larger angle splittings come first. The essential insight is that the scattering process appears differently depending on the hardness scale at which one examines it. At the hardest scale, the scale of the hard interaction, there are just a few partons (typically quarks and gluons). Then, as the hardness scale at which we examine the process decreases, these partons split, making more partons in a parton shower.222Thus, with respect to initial state partons, the shower evolution starts from the hard interaction and moves backward in time to softer initial state interactions. At any stage, a certain amount of structure has emerged, while softer structure remains unresolved.

In this paper, we start with the principle that a parton shower should fully reflect the infrared singularity structure of Feynman graphs for QCD and also the role of parton distribution functions in absorbing initial state singularities. Thus, we start with the infrared sensitive operator associated with the parton distributions and with a perturbative operator $\mathcal{D}(\mu^{2})$ that represents the infrared singularities of QCD Feynman diagrams. We connect $\mathcal{D}(\mu^{2})$ to both the shower splitting kernels and to the subtractions Daleo:2009yj ; Boughezal:2010mc ; Gehrmann:2011wi ; GehrmannDeRidder:2012ja ; Currie:2013vh ; Somogyi:2006da ; Somogyi:2006db ; Somogyi:2013yk needed to calculate a perturbative cross section beyond the leading order (LO). We work at arbitrary perturbative order. That is, we consider a hard scattering cross section calculated, with subtractions, at $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$ and a parton shower with $\mathrm{N}^{k-1}\mathrm{L}\mathrm{O}$ splitting functions. (This counts a LO shower as having splitting functions proportional to $\alpha_{\mathrm{s}}$.)

The construction that we present is based on the operator $\mathcal{D}(\mu^{2})$. This operator is to contain the infrared singularity structure of Feynman graphs for QCD. There is no unique recipe for constructing the $\alpha_{\mathrm{s}}^{n}$ contribution, $\mathcal{D}^{(n)}(\mu^{2})$, to $\mathcal{D}(\mu^{2})$. As described in Sec. V, one needs to specify a definition of hardness associated with the integrations in graphs, one needs a momentum mapping, and one needs to specify the form of the functions used as one moves away from the strict soft and collinear limits. At first order, we have made these choices, so that $d\mathcal{D}^{(1)}(\mu^{2})/d\log(\mu^{2})$ is part of Deductor. At higher orders, we do not attempt to construct the $\mathcal{D}^{(n)}(\mu^{2})$. Rather, we provide formulas for what to do once one has $\mathcal{D}^{(n)}(\mu^{2})$ for $n\leq k$.

The formalism uses another operator $\mathcal{V}(\mu^{2})$. This operator is obtained from $\mathcal{D}(\mu^{2})$ and factors associated with the parton distributions but it is obtained by integrating over all of the parton splitting variables, so that it is infrared finite. In a standard first order parton shower, the Sudakov exponent is quite directly related to the part of $\mathcal{V}(\mu^{2})$ that comes from one real parton emission. There is some freedom in setting the color and spin structure of $\mathcal{V}(\mu^{2})$. Thus we leave $\mathcal{V}(\mu^{2})$ partly unspecified.

The parton shower defined here needs to respect the structure of quantum field theory. Thus it includes quantum interference and maintains an exact accounting for the quantum spins and colors of the partons in the shower. The formulation is based on what we call the statistical space, introduce in Ref. NSI . It consists of states that describe the momenta, flavors, colors, and spins of any number of partons. The colors and spins are treated as fully quantum mechanical. This means that the statistical states are density matrices in the quantum color and spin space. We review the statistical space in Sec. IV.

Using the statistical space, we maintain an exact accounting for the quantum spins and colors. It is not known how to make this practical in computer code, particularly for color. Thus one needs separate approximations, such as the LC+ approximation for color that is used in Deductor. We view the issue of approximations for spin and color as separate from the construction of the shower with exact spin and color. We do not discuss spin and color approximations in this paper.

The final result of the construction, presented in Eq. (134), is

[TABLE]

Here \big{|}{\rho_{\textsc{h}}}\big{)} is a statistical state representing the hard scattering, calculated at order $\alpha_{\mathrm{s}}^{k+B}$, where $\alpha_{\mathrm{s}}^{B}$ is the order of the Born hard scattering process. The hard scattering statistical state \big{|}{\rho_{\textsc{h}}}\big{)} includes subtractions, as in a normal $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$ perturbative calculation. Then $\mathcal{F}(\mu_{\textsc{h}}^{2})$ is an operator on the statistical space that multiplies by appropriate parton distribution functions and a parton luminosity factor.

The next operator in Eq. (1),

[TABLE]

is a process independent operator on the statistical space that leaves the number of partons unchanged and provides perturbative corrections needed to keep the measured cross section correct to order $\alpha_{\mathrm{s}}^{k+B}$. This factor also sums threshold logarithms associated with the hard scattering statistical state. The threshold logarithms are an essential part of the construction. As we will discuss, they are included in Deductor, but they are not part of other leading order shower generators. We presented an earlier formulation of threshold summation in a leading order shower in Ref. NSThreshold . This formulation turned out to have certain flaws. In a companion paper NSThresholdII , we exhibit the practical effects of the threshold summation according to Eq. (1) (but with \big{|}{\rho_{\textsc{h}}}\big{)} evaluated at leading order only).

The next operator in Eq. (1),

[TABLE]

is a process independent operator that creates more partons in a parton shower. In the first order case, $k=1$, this is a rather standard parton shower if we average over spins and take the leading color approximation. In general, the splitting generator $\mathcal{S}(\mu^{2})$ consists of terms that are of order $\alpha_{\mathrm{s}}^{n}$ with $1\leq n\leq k$. The shower starts at a hardness scale $\mu_{\textsc{h}}^{2}$ appropriate for the hard scattering and ends at smaller hardness scale $\mu_{\mathrm{f}}^{2}$ that should be large enough so that perturbation theory at this scale can still be trusted.

The final operator in Eq. (1), $\mathcal{O}_{J}$, specifies the infrared safe measurement that one wants to make on the parton state after the shower. The hardness scale associated with this measurement is $Q[J]^{2}$. Finally, \big{(}{1}\big{|} is an instruction to integrate over all of the parton variables.

The cross section $\sigma[J]$ is then correct to order $\alpha_{\mathrm{s}}^{k+B}$ and includes a version of the cross section beyond this order, within the approximations of a parton shower. Notice that the property that the cross section including showering, $\sigma[J]$, is correct to order $\alpha_{\mathrm{s}}^{k+B}$ means that the shower is matched to an order $\alpha_{\mathrm{s}}^{k+B}$ perturbative calculation of $\sigma[J]$. This matching is an intrinsic part of the shower formulation.

We discuss a very general formulation of parton showers. However, we want to keep the notation simple, so, without loss of generality, we use Higgs boson production as an example. We use five flavors of quarks. In practical applications, one uses a variable flavor number scheme in which non-zero values of $m_{\mathrm{b}}$ and $m_{\mathrm{c}}$ appear. However, this creates complications, especially if we want to work at an arbitrary order of perturbation theory. Thus in this paper we set $m_{\mathrm{b}}=m_{\mathrm{c}}=0$.

We have found it useful to change some of the notation that we used in our previous papers in order to address a much more general problem. We hope that this does not cause confusion. We provide a translation in table 1.

Following this brief overview, we include a brief Sec. III on factorization, which plays an important role in the conceptual development. Then we devote Sec. IV to partons and the spin and color density operator, which we use to define the statistical space. This leads us describe the perturbative cross section and the infrared sensitive operator $\mathcal{D}(\mu^{2})$ in Sec. V. Then in Sec. VI we manipulate the perturbative cross section to define the parton shower. Sec. VII presents a summary and outlook.

There are four appendices. Appendix A presents a toy model for the operators used in the construction. We hope that this concrete model will prove instructive. Appendix B discusses the definition of parton distribution functions needed for a shower. Appendix C discusses how certain scale parameters can be chosen dynamically instead of statically, as in the main text. Appendix D discusses the relation of $\overline{\mathrm{MS}}$ renormalization to the definition of the parton shower.

III Factorization

We consider a hard scattering process in the collisions of two high energy hadrons, A and B. The hadrons carry momenta $P_{\mathrm{A}}$ and $P_{\mathrm{B}}$. The hadron energies are high enough that we can simplify the equations describing the collision kinematics by treating the colliding hadrons as being massless. Then with a suitable choice of reference frame, the hadron momenta are

[TABLE]

Here we use momentum components $(p^{+},p^{-},\bm{p}_{\perp})$ with $p^{\pm}=(p^{0}\pm p^{3})/\sqrt{2}$. We then imagine a parton level process in which a parton from hadron A, with flavor $a$ and momentum $p_{\mathrm{a}}=\eta_{\mathrm{a}}P_{\mathrm{A}}$ collides with a parton from hadron B, with flavor $b$ and momentum $p_{\mathrm{b}}=\eta_{\mathrm{b}}P_{\mathrm{B}}$.

We will be interested in an inclusive cross section to create some hard state, for instance, a Z boson plus possibly jets, or just jets. We will use the production of a Higgs boson, $A+B\to\textsc{H}+{\it QCD\ partons}$ as our principle example. At the Born level, it is produced via the partonic process $\mathrm{g}+\mathrm{g}\to\textsc{H}$. We treat the Higgs boson as being stable and on shell. We denote the momentum of the Higgs boson by

[TABLE]

The collision also produces QCD partons with flavors $f_{i}$ and momenta $p_{i}$, with $i=1,\dots,m$. In this paper, we consider the QCD partons to be massless. Each final state parton has rapidity $y_{i}$ and transverse momentum $\bm{p}_{i,\perp}$, so that the components of its momentum are

[TABLE]

It is up to us to decide what we want to measure about the final state of our process. We can consider many cases at once by simply saying that we are interested in a cross section $\sigma[J]$ to measure an observable quantity $J$, leaving the definition of $J$ unspecified. We will see in the following subsection how $\sigma[J]$ can be specified for a general observable $J$. Then parton distribution functions relate $\sigma[J]$ to an analogous cross section $\hat{\sigma}[J]$ for the collision of two partons. In its briefest form, the relation is

[TABLE]

III.1 Infrared safety

We demand that the observable $J$ be infrared safe. To specify what that means, we write Eq. (7) in more detail:

[TABLE]

Here we start with the cross section to produce the Higgs boson plus $m$ partons with momenta

[TABLE]

We multiply the cross section by a function $J_{m}(\{p\}_{m})$ that specifies the measurement that we want to make on the final state partons. These functions are taken to be symmetric under interchange of the QCD momentum arguments $\{p_{1},\dots,p_{m}\}$. Accordingly, we divide by the number $m!$ of permutations of the QCD parton labels. We integrate over the momenta of the final state partons. The transverse momentum of the Higgs boson and the needed momentum fractions for the incoming partons are determined by momentum conservation. Finally, we sum over the number $m$ of final state QCD partons.

Infrared safety is a property of the functions $J_{m}$ that relates each function $J_{m+1}(\!\{p\}_{m+1}\!)$ to the function $J_{m}(\{p\}_{m})$ with one fewer parton. There are two requirements needed for $J$ to be infrared safe.

First, consider the limit in which partons $m+1$ and $m$ become collinear:

[TABLE]

Here $\tilde{p}_{m}$ is a lightlike momentum and $0<z<1$. We can concentrate on just partons with labels $m+1$ and $m$ because the functions $J$ are assumed to be symmetric under interchange of the parton labels. In order for $J$ to be infrared safe, we demand that

[TABLE]

in the collinear limit (LABEL:eq:collinearlimit).

Second, consider also the limit in which parton $m+1$ becomes collinear to one of the beams,

[TABLE]

or

[TABLE]

Here $0\leq\xi$. When $\xi=0$, parton $m+1$ is simply becoming infinitely soft. In order for $J$ to be infrared safe, we demand that

[TABLE]

in either limit (12) or (13).

Briefly, then, infrared safety means that the result of the measurement is not sensitive to whether or not one parton splits into two almost collinear partons and it is not sensitive to any partons that have very small momenta transverse to the beam directions.

III.2 A more quantitative view of infrared safety

We now discuss infrared safety a little more quantitatively. Consider, as above, two final state partons that are nearly collinear. This is modeled in a parton shower algorithm as a splitting of a parton with momentum $\tilde{p}_{m}$ into two partons with momenta $p_{m}$ and $p_{m+1}$. We can measure how close we are to the collinear limit by calculating333One could choose other hardness measures $\mu^{2}$. This one is based on the parton shower formulation in Deductor ; ShowerTime .

[TABLE]

Here $Q_{\textsc{H}}$ is a momentum vector that describes the hard scattering, for instance the momentum $p_{\textsc{H}}$ of the produced Higgs boson in our example. We define the parton energies in the rest frame of $Q_{\textsc{H}}$. The limit expressed in Eq. (LABEL:eq:collinearlimit) is $\mu_{\rm split}^{2}\to 0$.

Alternatively, we can consider a splitting of an initial state parton as modeled in a parton shower algorithm. It suffices to consider the splitting of an initial state parton in hadron A. The initial state parton with momentum $\tilde{p}_{\mathrm{a}}$ becomes a new initial state parton with momentum $p_{\mathrm{a}}$ and a new final state parton with momentum $p_{m+1}$.444This is in the “backwards evolution” picture. Going forward in time, the parton with momentum $p_{\mathrm{a}}$ splits. We can measure how close we are to the collinear limit by calculating

[TABLE]

Again, the limit expressed in Eq. (12) is $\mu_{\rm split}^{2}\to 0$.

We now suppose that we measure the cross section $\sigma[J]$ corresponding to an infrared safe measurement function $J$, using Eq. (8). When $\mu_{\rm split}^{2}=0$ in either example, application of Eq. (8) gives a cross section that we can call $\sigma_{0}[J]$. In applying Eq. (8), we can use the term with $m+1$ final state partons, with parton $m+1$ exactly collinear with either parton $m$ or parton “a.” Equivalently, we can use the $m$ parton state before the splitting. Because of Eq. (11) or Eq. (14), the result is exactly the same. Now, when $\mu_{\rm split}^{2}$ is small but not zero, we get a slightly different result, $\sigma[J]$. Let $\delta\sigma[J]=\sigma_{0}[J]-\sigma[J]$. We assume that $J_{m+1}(\{p\}_{m+1})$ is a smooth function of the parton momenta, at least near the soft or collinear limits. Then we will have $\delta\sigma[J]\to 0$ as $\mu_{\rm split}^{2}\to 0$. A typical case is $\delta\sigma[J]\propto\mu_{\rm split}^{2}$ as $\mu_{\rm split}^{2}\to 0$. Then we can define a scale $Q^{2}[J]$ that is characteristic of the observable by

[TABLE]

The ratio in Eq. (17) can be sensitive to the parton configuration, so we average over all configurations with the same $\mu_{\rm split}^{2}$. The scale $Q^{2}[J]$ measures how sensitive the cross section is to parton splittings.

An example may be helpful. Let $\sigma[J]$ measure the one-jet-inclusive cross section for jet transverse momentum $P_{\mathrm{J}}$ at small rapidity, using the anti-$k_{\mathrm{T}}$ jet algorithm with radius parameter $R$. Consider a very narrow jet in which one parton splits into two, one of which is soft. Then a simple estimate using the definitions above is

[TABLE]

Here $N$ measures how fast the jet cross section falls with increasing $P_{\mathrm{J}}$:

[TABLE]

This is a fairly large number, so that the effective scale $Q^{2}[J]$ for the jet measurement is smaller than $P_{\mathrm{J}}^{2}$. Additionally, $R$ is typically chosen to be less than one. The smaller $R$ is, the smaller $Q^{2}[J]$ is.

We see that we can understand infrared safety in terms of measurements that we might make in a parton shower simulation of a high energy scattering event. As the shower progresses, there are splittings. The corresponding values of $\mu_{\rm split}^{2}$ get smaller and smaller. If we measure $\sigma[J]$ at each stage of the shower, we will see that $\sigma[J]$ changes as the shower develops, but the changes get smaller and smaller according to

[TABLE]

This leads us to a limit on the possible accuracy of a perturbative approach to calculating $\sigma[J]$. When the parton shower has evolved to a 1 GeV scale, then we have reached the limits of perturbation theory. The fractional uncertainty associated with one more splitting is then

[TABLE]

The accuracy of the perturbative calculation can never be better than this.

III.3 Factorization

We can now state how one calculates the cross section for whatever observable $J$ we want – as long as $J$ is infrared safe. The formula we use was stated in Eq. (7) and we restate it here in a slightly more detailed form factorization :

[TABLE]

The intuitive basis for this is very simple. The factor $f_{a/A}(\eta_{\mathrm{a}},\mu^{2})\,d\eta_{\mathrm{a}}$ represents the probability to find a parton of flavor $a$ in a hadron of flavor $A$. For the other hadron, the corresponding probability is $f_{b/B}(\eta_{\mathrm{b}},\mu^{2})\,d\eta_{\mathrm{b}}$. Then $\hat{\sigma}[J]$ is the cross section to obtain the observable $J$ from the scattering of these partons, as given in Eq. (8). Naturally, this parton level cross section depends on the parton variables $a,b,\eta_{\mathrm{a}},\eta_{\mathrm{b}}$. Here the differential cross sections to produce $m$ final state partons contain delta functions that relate the momentum fractions $\eta_{\mathrm{a}}$ and $\eta_{\mathrm{b}}$ to the final state parton momenta. The parton distributions depend on a scale $\mu$. This is often called the factorization scale $\mu_{\mathrm{F}}$ and distinguished from the argument of $\alpha_{\mathrm{s}}$ and other running couplings, which is called the renormalization scale $\mu_{\mathrm{R}}$. In order to keep our notation simple, we set $\mu_{\mathrm{F}}=\mu_{\mathrm{R}}=\mu$.

The cross section $\hat{\sigma}[J]$ has a perturbative expansion in powers of $\alpha_{\mathrm{s}}(\mu^{2})$. That is

[TABLE]

Here we do not display the factors of $\alpha_{\mathrm{s}}$ or $\alpha_{\rm ew}$ that appear in the Born level cross section $\hat{\sigma}^{(0)}$. Perturbative calculations can be at lowest order (LO), corresponding to one term in the expansion, next-to-lowest order (NLO) with two terms, sometimes NNLO, and, in general, $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$.

One useful property is that the dependence of the calculated cross section on $\mu^{2}$ diminishes as we go to higher orders. Indeed, the cross section in nature, $\sigma[J]$, does not depend on $\mu^{2}$. Thus if we calculate to order $\alpha_{\mathrm{s}}^{k}$, the derivative of the calculated cross section with respect to $\mu^{2}$ will be of order $\alpha_{\mathrm{s}}^{k+1}$.

There is an error term in Eq. (22). No matter how many terms are included in $\hat{\sigma}$, there are contributions that are left out. These terms are suppressed by a power of $m\sim 1\ \mathrm{GeV}$ divided by a large scale parameter $Q$ that characterizes the hard scattering process to be measured. These contributions arise from the approximations needed to derive Eq. (22). For instance if a loop momentum $l$ flows through the wave function of quarks in a proton, we have to neglect $l$ compared to the hard momenta, say the transverse momentum of an observed jet. Not much is known about the general form of the power corrections for hadron-hadron collisions. In the rest of this paper, we will assume that the power $n$ in Eq. (22) is $n=2$ and we use $\sqrt{Q^{2}[J]}$ from Eq. (20) for $Q$. However, even if we lack a good estimate of the power corrections, it is important that they are there. If $Q$ is of order 100 GeV, then the power corrections are completely negligible. However, if $Q$ is of order 5 GeV, then we ought not to claim 1% accuracy in the calculation of $\sigma[J]$, no matter how many orders of perturbation theory we use.

IV Partons and the density operator

We will describe perturbative calculations of cross sections, how these are connected to the parton shower description of these same cross sections, and how this is connected to factorization. We begin in this section with definitions that we need to describe the evolution of a parton shower. We follow the framework of Ref. NSI .

IV.1 Amplitudes and the density operator in spin and color

In a perturbative calculation of a cross section, one constructs an amplitude \big{|}{M(\{p,f\}_{m})}\big{\rangle}. This amplitude depends on the momenta and flavors of two initial state partons, whatever outgoing electroweak partons there are, and $m$ outgoing QCD partons. For our example of Higgs boson production, the momentum and flavor observables are

[TABLE]

The partons carry spin and color, so the amplitude is a vector in the partonic spin and color space for $m$ final state QCD partons plus the two incoming partons, as indicted by the representation of the amplitude as a ket vector \big{|}{M(\{p,f\}_{m})}\big{\rangle}. We use basis vectors \big{|}{\{s\}_{m}}\big{\rangle} for the partonic spin space and basis vectors \big{|}{\{c\}_{m}}\big{\rangle} for the partonic color space.555The spin basis vectors can be chosen in a very simple way, but it is not so trivial to choose useful basis vectors for the color space. The choice that we make for Deductor is specified in Ref. NSI . With this choice, the color basis vectors are not exactly normalized and the basis vectors for different colors are not exactly orthogonal. To describe the evolution of a parton shower, it is useful to use quantum statistical mechanics, keeping the full quantum nature of the colors and spins. Thus we use the density operator in the color and spin space. The density operator is a linear combination of basis operators \big{|}{\{s,c\}_{m}}\big{\rangle}\big{\langle}{\{s^{\prime},c^{\prime}\}_{m}}\big{|}. Thus to describe the system at a certain stage of evolution, we use a function $\rho$ of the number $m$ of final state QCD partons and of the momenta and flavors $\{p,f\}_{m}$, where the value of $\rho$ is a color-spin density operator. That is

[TABLE]

The interpretation of this is that the differential probability $dP$ for there to be $m$ final state QCD partons with momenta and flavors $\{p,f\}_{m}$, times the expectation value of an operator $O$ on the color and spin space, is

[TABLE]

where

[TABLE]

Here we use dimensional regularization with $d=4-2\epsilon$.

The set of all such functions $\rho$ constitutes a vector space, which we call the statistical space. We represent the vector $\rho$ as a ket vector, \big{|}{\rho}\big{)}. The rounded end of the ket is meant to distinguish a vector in the statistical space from a vector in the quantum spin$\otimes$color space.

Notice that we use the symbol $\rho$ for four different but related concepts. First, for each choice of $m$ and $\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}$, $\rho(\{p,f,s,s^{\prime},c,c^{\prime}\}_{m})$ is a complex number. Second, for each choice of $m$ and $\{p,f\}_{m}$, $\rho(\{p,f\}_{m})$ is a linear operator on the quantum spin$\otimes$color space for the $m+2$ partons. Third, $\rho$ is a linear map from $m$ and $\{p,f\}_{m}$ to the space of operators on the quantum spin$\otimes$color space. Fourth, \big{|}{\rho}\big{)} is this linear map considered as an element of a vector space, the statistical space. This may seem complicated, but in the end we use almost entirely the statistical vectors \big{|}{\rho}\big{)}. This then gives us what we think is a compact and powerful notation.

We can define basis vectors \big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} in the statistical space in such a way that

[TABLE]

The completeness relation for the basis vectors is

[TABLE]

IV.2 Making an inclusive measurement

There is a special vector \big{(}{1}\big{|} defined by

[TABLE]

With this definition,

[TABLE]

is the total probability associated with the statistical state \big{|}{\rho}\big{)} as defined in Eq. (26) with $O=1$.

With this notation, we begin with perturbatively calculated amplitudes \big{|}{M(\{p,f\}_{m})}\big{\rangle} for just a few partons. Thus we begin with a perturbatively calculated vector \big{|}{\rho}\big{)} in the statistical space. Then we use perturbative operations that are represented as linear operators on the statistical space. Similarly, the measurement $J$ in Sec. III.1 is represented as a linear operator on the statistical space. Finally, multiplication by \big{(}{1}\big{|} allows us to obtain the expectation value of the measurement operator.

IV.3 Scales

The formalism also uses a reference vector $Q_{\textsc{H}}$ and several scales: a renormalization scale $\mu_{\mathrm{R}}^{2}$, a factorization scale $\mu_{\mathrm{F}}^{2}$, and an ultraviolet cutoff scale $\mu_{\mathrm{s}}^{2}$. For simplicity, all of these scales are set to a single scale $\mu^{2}$. The vector $Q_{\textsc{H}}$ is used to set the value $\mu_{\textsc{h}}^{2}$ of the common scale associated with the hard state \big{|}{\rho_{\textsc{h}}}\big{)}, defined later in Eq. (58): $\mu_{\textsc{h}}^{2}=Q_{\textsc{H}}^{2}$. A parton shower needs a measure of hardness of parton splittings. We also use $Q_{\textsc{H}}$ as a vector to help define one possible measure of hardness, $\Lambda^{2}$ defined in Eq. (55). For this purpose, $Q_{\textsc{H}}$ should be roughly in the direction of $p_{\mathrm{a}}+p_{\mathrm{b}}$ in an imagined initial hard scattering. The simplest way to define $Q_{\textsc{H}}$ is to use the hard core part of the intended measurement. For instance, if we are looking for the cross section to produce a Higgs boson with rapidity near zero, we can take $Q_{\textsc{H}}^{2}=m_{\textsc{h}}^{2}$ with zero rapidity and zero transverse part for $Q_{\textsc{H}}$. We will mention another, more dynamic, way to define $Q_{\textsc{H}}$ later, in Appendix C.

IV.4 Multiplying by parton distribution functions

In order to turn matrix elements into cross sections, we divide by a parton luminosity factor666We have noticed that this factor was too small by a factor 2 in Eq. (3.15) of NSI . $n_{\mathrm{c}}(a)n_{\mathrm{s}}(a)\,n_{\mathrm{c}}(b)n_{\mathrm{s}}(b)\ 4p_{\mathrm{a}}\cdot p_{\mathrm{b}}$. Here the color counting factor is $n_{\mathrm{c}}(f)=3$ for a quark flavor $f$ and $n_{\mathrm{c}}(f)=8$ for $f=\mathrm{g}$. The spin counting factor is $n_{\mathrm{s}}(f)=2$ for a quark flavor $f$ and $n_{\mathrm{s}}(f)=2(1-\epsilon)$ for $f=\mathrm{g}$ when we work in $4-2\epsilon$ dimensions in order to regularize infrared divergences. Then we need to multiply by a parton distribution factor

[TABLE]

The parton distribution functions here could be the five flavor $\overline{\mathrm{MS}}$ parton distribution functions, or they could have a different definition. We combine the two parton distribution functions as one operator $\mathcal{F}(\mu^{2})$ that acts on the statistical space:

[TABLE]

We sometimes need a more general parton factor in which the parton distributions are convolved with a function that is also a matrix in the parton flavors. For instance, the evolution equation for the product of parton distribution functions is

[TABLE]

where the prime denotes differentiation and $F^{\prime}$ is given by the convolution product of $F$ with an evolution kernel $P$,

[TABLE]

The precise definition is777Based on the order of the flavor indices, it would be more conventional to write this as $P\circ F$, but we believe that the notation in Eq. (35) better expresses the physics in the context of this paper.

[TABLE]

The evolution kernel for the product of parton distributions is the sum of parton evolution kernels for each of the two parton distribution functions:

[TABLE]

This gives us an operator $\mathcal{F}^{\prime}(\mu^{2})$ defined by

[TABLE]

In place of $\mathcal{F}^{\prime}(\mu^{2})$, we use a notation that directly displays that $F^{\prime}_{a,b}(\eta_{\mathrm{a}},\eta_{\mathrm{b}},\mu^{2})$ is constructed according to Eq. (36),

[TABLE]

The circle indicates the convolution and the square brackets $[\cdots]$ indicate what is included in the convolution.

V The perturbative cross section

We now consider the cross section for some hard process. We can use any hard process that can lead to an infrared safe cross section, but the details of the notation depend on what hard process we consider. In order to keep the notation simple, we consider a specific process, the cross section make a Higgs boson plus QCD partons, calculated at $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$. We let $J$ represent an infrared safe measurement of interest as in Sec. III.1. The hardest scale associated with this measurement is $\mu_{\textsc{h}}^{2}$. The scale $Q^{2}[J]$ introduced in Eq. (17) could be much smaller, as long as it is large compared to $1\ \mathrm{GeV}^{2}$. In that case, the sort of fixed order calculation discussed in this section is not so useful because there are large logarithms, $\log(\mu_{\textsc{h}}^{2}/Q^{2}[J])$. Thus in this section, it may be helpful to imagine that $Q^{2}[J]$ is not very different from $\mu_{\textsc{h}}^{2}$.

V.1 The Born cross section

Let us begin at the Born level. We call the statistical state corresponding to the Born level matrix element \big{|}{\rho^{(0,0)}(\mu^{2})}\big{)}. The measurement $J$ can be represented as an operator $\mathcal{O}_{J}$. To make a cross section, we need a luminosity factor and parton distribution functions. We define an operator \mathcal{F}_{\!\raisebox{-2.0pt}{{\overline{\phantom{\raisebox{-0.6pt}{$\scriptstyle{\textsc{ms}}$}}}\hskip-9.5pt\scriptstyle{\textsc{ms}}}}}(\mu^{2}) as in Eq. (33) by

[TABLE]

We use $\overline{\mathrm{MS}}$ parton distributions for five flavors. Then the cross section for measurement operator $\mathcal{O}_{J}$ is

[TABLE]

V.2 The higher order cross section

Now we consider higher orders of perturbation theory. We define $\alpha_{\mathrm{s}}$ and the parton distribution functions that we start with using $\overline{\mathrm{MS}}$ renormalization, as reviewed in Appendix D. Including terms up to order $\alpha_{\mathrm{s}}^{k}$, the $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$ cross section can be written as

[TABLE]

Here we convolve \mathcal{F}_{\!\raisebox{-2.0pt}{{\overline{\phantom{\raisebox{-0.6pt}{$\scriptstyle{\textsc{ms}}$}}}\hskip-9.5pt\scriptstyle{\textsc{ms}}}}}(\mu^{2}) with the renormalization factor $\mathcal{Z}_{F}(\mu^{2})$ for five flavor $\overline{\mathrm{MS}}$ parton distribution functions.888This is for massless partons. Conceptually, $\mathcal{Z}_{F}(\mu^{2})$ should be understood as the inverse, in the sense of convolutions, of the product of two parton-in-a-parton distribution functions with on-shell massless incoming parton states. Then this factor removes infrared poles from the cross section. At the level of bare operators, $f^{\rm bare}_{a/b}(\xi)=\delta_{ab}\delta(1-\xi)$. Convolving the renormalized parton-in-a-parton distribution functions with $\mathcal{Z}_{F}$ gives the bare distribution functions, leading to Eq. (42). See Appendix D. This factor is just 1 at order $\alpha_{\mathrm{s}}^{0}$. Its higher order contributions are defined by working in $4-2\epsilon$ dimensions and consist of the $1/\epsilon^{n}$ pole terms needed to remove the ultraviolet poles from the renormalized operator that defines parton distributions. Eq. (42) is a perturbative formula. We are to expand all of the factors up to the desired order $\alpha_{\mathrm{s}}^{k}$ and neglect the remainder, indicated by the error estimate $\mathcal{O}(\alpha_{\mathrm{s}}^{k+1})$.

The statistical state \big{|}{\rho(\mu^{2})}\big{)} has a perturbative expansion

[TABLE]

The order $n$ contribution to the statistical state is a sum of terms,999It may be helpful to note that we define \big{|}{\rho^{(n)}(\mu^{2})}\big{)} using on-shell matrix elements \big{|}{M(\{p,f\}_{m})}\big{\rangle} and their complex conjugates, including the factors needed to make \big{|}{M(\{p,f\}_{m})}\big{\rangle} into an S-matrix element. Then the right hand side of Eq. (43) is invariant under changes of the renormalization scale $\mu^{2}$ up to order $\alpha_{\mathrm{s}}^{k+1}$.

[TABLE]

In \big{|}{\rho^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})}\big{)}, there are $n_{\textsc{r}}$ final state partons and $n_{\textsc{v}}$ virtual loops. This is for Higgs boson production as the Born process. If we had chosen two jet production as the Born level hard process, then there would be $2+n_{\textsc{r}}$ partons in the final state.

The $n_{\textsc{v}}$ virtual loops in \big{|}{\rho^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})}\big{)} can each produce $1/\epsilon$ and $1/\epsilon^{2}$ poles. The $n_{\textsc{r}}$ partons in the final state can give soft and collinear singularities. The statistical state vector is singular when any of these partons become soft or collinear with the beam directions or collinear with each other. In the case of a Born process that has final state QCD partons, there are also singularities when any of the $n_{\textsc{r}}$ additional partons becomes collinear with starting final state partons. In these singular regions, the infrared safe measurement operator $\mathcal{O}_{J}$ sees the partons that are collinear to a given direction as equivalent to a single parton and it does not see partons that are soft or collinear to the beam directions at all. Thus, it is as if we had a completely inclusive measurement as defined by left multiplying by \big{(}{1}\big{|}. Then we again get $1/\epsilon$ and $1/\epsilon^{2}$ poles. Many of the poles cancel between real and virtual graphs. There are, however, some poles associated with momenta that are collinear with the initial state parton momenta. These cancel the poles in $\mathcal{Z}_{F}(\mu^{2})$. We are left with a finite result.

V.3 Introduction of the infrared sensitive operator $\mathcal{D}(\mu^{2})$

The formula (42) is not completely practical because each term \big{|}{\rho^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})}\big{)} generates infrared singularities or poles. Only at the end of a calculation, which includes some complicated integrations that include the measurement function, do the poles produced by the infrared singularities cancel. To make this more practical, we define a certain operator $\mathcal{D}(\mu^{2})$ and insert a factor $\mathcal{D}(\mu^{2})\mathcal{D}^{-1}(\mu^{2})$ into Eq. (42), giving

[TABLE]

The operator $\mathcal{D}(\mu^{2})$ depends on the dimensional regularization parameter $\epsilon$, but we do not display this dependence explicitly. It depends on a second scale $\mu_{\mathrm{s}}^{2}$ along with $\mu^{2}$, but we set $\mu^{2}=\mu_{\mathrm{s}}^{2}$.

The idea behind $\mathcal{D}(\mu^{2})$ is that a contribution to \big{|}{\rho(\mu^{2})}\big{)} has poles from virtual loops and has singularities when some of its external lines become collinear or soft. This is illustrated in Fig. 1. It is simplest to think of the graphs depicted as being in a physical gauge. There are two hard subgraphs, represented as yellow blobs, one for the amplitude and one for the conjugate amplitude. The subgraphs can be tree graphs or can contain virtual loops. We suppose that everything inside the hard subgraphs is harder than a reference scale $\mu_{\rm hard}^{2}$. That is, all of the internal propagators are far off shell. Two initial state lines and $m$ final state parton lines emerge from the hard subgraph. Here $m=3$. Then there are additional interactions. Some number $n_{\textsc{r}}$ of additional partons are emitted and $n_{\textsc{v}}$ parton lines are exchanged. Here $n_{\textsc{r}}=1$ and $n_{\textsc{v}}=1$. The external parton momenta are labeled $\{\hat{p}\}_{m+n_{\textsc{r}}}$. There is an infrared divergence when the the virtual gluon becomes soft and there are singularities when the real gluon momentum becomes soft or collinear to the antiquark line while the virtual gluon momentum is becoming soft.

We want to capture the structure of these infrared singularities, as illustrated in Fig. 2. We note that near the singularities, the partons emerging from the hard subgraphs are nearly on shell and their momenta lie in the directions of the external parton momenta. The momenta carried by lines internal to the hard subgraphs are almost unchanged. Therefore, we can approximate the graph by letting the momenta $\{p\}_{m}$ of the partons emerging from the hard subgraphs be exactly on shell. Their momenta are given as functions of the momenta $\{\hat{p}\}_{m+n_{\textsc{r}}}$ of the external partons. Here, we need to define a momentum mapping $\{\hat{p}\}_{m+n_{\textsc{r}}}\to\{p\}_{m}$. Then we can approximate the original graph by a hard part \big{(}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{|}{\rho_{\rm hard}(\mu^{2})}\big{)} and a singular factor. We call the singular factor \big{(}{\{\hat{p},\hat{f},\hat{s},\hat{s}^{\prime},\hat{c},\hat{c}^{\prime}\}_{m+n_{\textsc{r}}}}\big{|}\mathcal{D}(\mu^{2})\big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)}. The singular factor is derived just from the singular part of the graph. It is independent of what is in the hard part.

We thus assert that any such set of poles and singularities can be organized into a hard subgraph, \big{|}{\rho_{\rm hard}(\mu^{2})}\big{)}, convolved with a singular factor:

[TABLE]

The division between singular and hard factors depends on the singularity to be examined. In the hard factor, the external momenta $\{p\}_{m}$ and any internal loop momenta are to be hard at some scale that we can call $\mu^{2}_{\rm hard}$. This means that they are not closely collinear to each other or soft at scales softer than $\mu^{2}_{\rm hard}$.

The singular factor is typically represented as separate factors labeled ${\it soft}$ and ${\it jet}_{i}$, where ${\it soft}$ can include Glauber exchanges and two of the jets are in the beam directions. However, we do not need to separate the singular factor into separate subfactors. An early and instructive analysis of the singularities of QCD was given by Libby and Sterman LibbySterman . An extensive modern analysis can be found in Collins JCCbook . For one real gluon emission, an example of $\mathcal{D}(\mu^{2})$ can be defined from the Catani-Seymour dipole splitting functions CataniSeymour . The operators $\bm{Sp}$ of Catani, de Florian, and Rodrigo CataniFactViolated are also closely related to $\mathcal{D}(\mu^{2})$ for certain cases. At one loop, our version of $d\mathcal{D}(\mu^{2})/d\log\mu^{2}$ is implemented in Deductor.

The operator $\mathcal{D}(\mu^{2})$ has a perturbative expansion

[TABLE]

The order $n$ contribution, $\mathcal{D}^{(n)}(\mu^{2})$, is a sum of infrared sensitive operators,

[TABLE]

Acting on a state \big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} with $m$ final state QCD partons, $\mathcal{D}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})$ produces a state with $m+n_{\textsc{r}}$ final state QCD partons with momenta $\{\hat{p}\}_{m+n_{\textsc{r}}}$. There are integrations over the loop momenta $\{\ell_{1},\dots,\ell_{n_{\textsc{v}}}\}$ of $n_{\textsc{v}}$ virtual loops.

The operator $\mathcal{D}(\mu^{2})$ depends on a hardness scale $\mu^{2}_{\mathrm{s}}$ that defines an infrared sensitive region $R(\mu^{2}_{\mathrm{s}})$ in the space of the momenta $\{\hat{p}\}_{m+n_{\textsc{r}}}$ and $\{\ell_{1},\dots,\ell_{n_{\textsc{v}}}\}$. We can think of $\mu^{2}_{\mathrm{s}}$ as being comparable to the scale $\mu^{2}_{\rm hard}$ of \big{|}{\rho_{\rm hard}(\mu^{2})}\big{)}. The infrared sensitive region $R(\mu^{2}_{\mathrm{s}})$ surrounds the leading singularity, at which each of the momenta $\{\hat{p}\}_{m+n_{\textsc{r}}}$ and $\{\ell_{1},\dots,\ell_{n_{\textsc{v}}}\}$ is soft or collinear to one of the input momenta $\{p\}_{m}$. In the case of the output momenta $\{\hat{p}\}_{m+n_{\textsc{r}}}$, this means that, at the leading singularity, these momenta form $m$ infinitely narrow jets with momenta $\{p\}_{m}$. Of course, if we look just a little bit away from the limit of infinitely narrow jets, we see that the jets can have subjets. The singularity structure of the subjets, including both real and virtual momenta, is included in $\mathcal{D}(\mu^{2})$.

We need the infrared sensitive region $R(\mu^{2}_{\mathrm{s}})$ because, when we form $\mathcal{D}(\mu^{2})$ by making approximations that apply near the leading singularity, we necessarily simplify the behavior away from this singularity. We introduce cuts such that $\mathcal{D}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})$ gets contributions only from inside $R(\mu^{2}_{\mathrm{s}})$. The leading singularity is inside the region for any $\mu^{2}_{\mathrm{s}}$, but for larger $\mu^{2}_{\mathrm{s}}$, the region is larger. Of course, there is more than one way to introduce $\mu_{\mathrm{s}}^{2}$.

For simplicity, we set the renormalization and factorization scale $\mu^{2}$ equal to $\mu_{\mathrm{s}}^{2}$.

We will say more about the $\mathcal{D}^{(n)}(\mu^{2})$ later, although we do not construct them. For now we simply assume that they are available and investigate how they can be used to construct subtractions for a fixed order calculation and splitting functions for a parton shower. See Appendix A for an example of the operators $\mathcal{D}^{(n)}(\mu^{2})$ in a toy model.

V.4 Subtractions for the perturbative cross section

Given $\mathcal{D}(\mu^{2})$, we can construct

[TABLE]

The perturbative coefficients $\widetilde{\mathcal{D}}^{(n)}(\mu^{2})$ are defined by

[TABLE]

This gives, for instance,

[TABLE]

or for the higher orders

[TABLE]

The order $n$ contribution, $\widetilde{\mathcal{D}}^{(n)}(\mu^{2})$, is a sum of operators,

[TABLE]

Acting on a state \big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} with $m$ partons, $\widetilde{\mathcal{D}}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})$ produces a state with $m+n_{\textsc{r}}$ partons while integrating over $n_{\textsc{v}}$ virtual loops. One constructs $\widetilde{\mathcal{D}}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})$ using Eq. (50).

The operator $\mathcal{D}^{-1}(\mu^{2})$ is very useful. In Eq. (45), we have the factor \mathcal{D}^{-1}(\mu^{2})\mathcal{O}_{J}\big{|}{\rho(\mu^{2})}\big{)}. We are to expand this product in powers of $\alpha_{\mathrm{s}}$, keeping the terms up to order $\alpha_{\mathrm{s}}^{k}$. The factor \mathcal{O}_{J}\big{|}{\rho(\mu^{2})}\big{)} has infrared singularities and poles, but the operator $\mathcal{D}^{-1}(\mu^{2})$ removes them. The simple argument is that \mathcal{D}(\mu^{2})\big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} contains the infrared singularities and poles produced by QCD from a hard parton state \big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} but, according to Eq. (50), when we apply $\mathcal{D}^{-1}(\mu^{2})$ to this state, the poles and singular terms are cancelled. That is, $\mathcal{D}^{-1}(\mu^{2})$ provides the subtraction terms that we need to remove the singularities and poles from an $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$ perturbative calculation.

To understand what the operator $\mathcal{D}^{-1}(\mu^{2})$ does, it is helpful to examine the familiar case of an NLO calculation. At this order, Eq. (45) becomes

[TABLE]

We have defined $\mathcal{D}^{(0,1)}(\mu^{2})$ so that it leaves the number of partons unchanged and so that it has infrared poles. Furthermore, the poles in \mathcal{D}^{(0,1)}(\mu^{2})\big{|}{\rho^{(0,0)}}\big{)} should directly cancel those of \big{|}{\rho^{(0,1)}(\mu^{2})}\big{)}. We have defined $\mathcal{D}^{(1,0)}(\mu^{2})$ so that, acting on the state \big{|}{\rho^{(0,0)}(\mu^{2})}\big{)}, it adds one parton and so that when this parton is soft or nearly collinear with one of the existing partons (in our example, the initial state partons) \mathcal{D}^{(1,0)}(\mu^{2})\mathcal{O}_{J}\big{|}{\rho^{(0,0)}(\mu^{2})}\big{)} approaches \mathcal{O}_{J}\big{|}{\rho^{(1,0)}(\mu^{2})}\big{)}. In a standard application, one performs the integrations over the momentum of the emitted parton numerically. The integrand in the subtraction cancels the integrand in \big{|}{\rho^{(1,0)}(\mu^{2})}\big{)} in the infrared region, so that one obtains a convergent integration. Having subtracted the operators $\mathcal{D}^{(1,0)}$ and $\mathcal{D}^{(0,1)}$, we add them back as part of $\mathcal{D}$. Now, in a standard application, all of the integrations corresponding to the first line of Eq. (54) are performed analytically. All of the $1/\epsilon^{2}$ and $1/\epsilon$ poles cancel and we are left with a completely finite order $\alpha_{\mathrm{s}}$ contribution to the cross section. Note that the contribution from the first line beyond just the parton distribution functions is infrared finite, but it is not zero. It forms a significant part of the NLO calculation. In Eq. (54), we have a product of one infrared finite object times another, each expanded to order $\alpha_{\mathrm{s}}$. In the customary NLO calculation, one expands the product to order $\alpha_{\mathrm{s}}$ and drops the $\alpha_{\mathrm{s}}^{2}$ term, but one could keep the $\alpha_{\mathrm{s}}^{2}$ term if desired.

We note that if we had wanted to use $\mathcal{D}(\mu^{2})$ for the single purpose of defining subtractions for the hard scattering, \big{|}{\rho(\mu^{2})}\big{)}, we could have used a fixed scale. We need an adjustable scale $\mu_{\mathrm{s}}^{2}$ to use $\mathcal{D}(\mu^{2})$ in a shower algorithm because the hardness scale of the shower changes as the shower progresses.

V.5 Properties of the infrared sensitive operator

We can now say a little more about the infrared sensitive operator $\mathcal{D}(\mu^{2})$, without giving a detailed specification. This operator is decomposed into operators $\mathcal{D}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})$ according to Eqs. (47) and (48). Acting on a state \big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} with $m$ partons, $\mathcal{D}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})$ produces a state with $m+n_{\textsc{r}}$ partons while adding $n_{\textsc{v}}$ virtual loops. The resulting state can be expanded in basis states for $m+n_{\textsc{r}}$ partons, \big{|}{\{\hat{p},\hat{f},\hat{s},\hat{s}^{\prime},\hat{c},\hat{c}^{\prime}\}_{m+n_{\textsc{r}}}}\big{)}. There is then an invertible mapping between the new momenta $\{\hat{p}\}_{m+n_{\textsc{r}}}$ and the starting momenta $\{p\}_{m}$ together with a set of splitting variables $\zeta_{p}({n_{\textsc{r}}})$. One can choose what this mapping is.

To make this a little more concrete, consider $\mathcal{D}^{(1,0)}(\mu^{2})$ with one parton emitted. There is one term in the emission probability for each final state parton $l\in\{1,\dots,m\}$ and one for each initial state parton $l\in\{\mathrm{a},\mathrm{b}\}$. We think of $l$ as the emitting parton and let the emission probability be singular when $\hat{p}_{m+1}$ becomes collinear with $\hat{p}_{l}$. The emission probabilities are also singular in the limit in which $\hat{p}_{m+1}$ becomes soft, $\hat{p}_{m+1}\to 0$. Deductor uses something similar to the Catani-Seymour CataniSeymour dipole splitting functions to construct $\mathcal{D}^{(1,0)}(\mu^{2})$. The splitting variables $\zeta_{p}$ are taken to be an azimuthal angle $\phi$, a momentum fraction $z$, and a measure of the hardness of the splitting. In Deductor, the hardness variable is the virtuality of the splitting divided by the energy of the mother parton,

[TABLE]

where the vector $Q_{\textsc{H}}$ is defined globally as described in Sec. IV.3. Alternatively, it can be defined dynamically as described in Appendix C.

There is freedom to choose the functional form of $\mathcal{D}^{(1,0)}(\mu^{2})$ away from the limits of soft and collinear emissions.

There is also freedom to choose the momentum mapping. The simplest case is a splitting of a parton $l$ into two partons $l$ and $m+1$. Then we cannot have $p_{l}$ be the same as $\hat{p}_{l}+\hat{p}_{m+1}$ with $p_{l}^{2}=\hat{p}_{l}^{2}=\hat{p}_{m+1}^{2}=0$, so the momentum mapping has to take a some momentum from the other partons and supply it to $\hat{p}_{l}+\hat{p}_{m+1}$. In Deductor, we use a global mapping, taking a small amount of momentum from each of the other partons. A second possibility beyond a simple splitting is interference between emission of a gluon $m+1$ from parton $l_{\mathrm{L}}$ in the ket state and emission from another parton $l_{\mathrm{R}}$ in the bra state. In this case, $\mathcal{D}^{(1,0)}(\mu^{2})$ in Deductor is a linear combination of contributions that use the momentum mappings for a splitting of parton $l_{\mathrm{L}}$ and for a splitting of parton $l_{\mathrm{R}}$. The coefficients in the linear combination are a “dipole partitioning” function $A^{\prime}$ that is specified in Deductor.

The operator $\mathcal{D}^{(1,0)}(\mu^{2})$, acting on a state \big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)}, produces a state with one more parton, parton $m+1$. It is crucial that there be an ultraviolet cutoff for $\hat{p}_{m+1}$. The cutoff is specified by a parameter that we call $\mu_{\mathrm{s}}^{2}$. In Deductor, we use $\Lambda^{2}$ as given in Eq. (55) to define the cutoff. In $\mathcal{D}^{(1,0)}(\mu^{2})$, we require

[TABLE]

A similar cutoff applies inside the integration for a virtual loop in $\mathcal{D}^{(0,1)}(\mu^{2})$. Defining this cutoff is more involved than we can review here. The calculations are described in Refs. NSThreshold ; NSThresholdII .

VI From the perturbative cross section to a parton shower

In this section, we begin with Eq. (45) for the perturbative cross section. We set the scale to $\mu_{\textsc{h}}^{2}$, which we take to be equal to $Q_{\textsc{H}}^{2}$. Here $Q_{\textsc{H}}^{2}$ is defined in Sec. IV.3 to be a fixed vector, although we can use a dynamical definition as described in Appendix C. We now seek a more powerful formulation that will enable us to use more general measurement operators $\mathcal{O}_{J}$ for which a perturbative expansion of the cross section might contain large logarithms of the generic form $\alpha_{\mathrm{s}}^{n}\log^{2n}(k^{2}/Q^{2})$. Often, a parton shower can approximately sum such logarithms.

As just stated, we set the scale $\mu^{2}$ in Eq. (45) to $\mu_{\textsc{h}}^{2}$. In our example in which the Born process is Higgs boson production, we might choose $\mu_{\textsc{h}}^{2}=m_{\textsc{h}}^{2}$. This affects the scale at which $\alpha_{\mathrm{s}}$ and the parton distribution functions are evaluated. It also affects the upper cutoff on the scale of emissions in the subtraction terms in $\mathcal{D}^{-1}(\mu_{\textsc{h}}^{2})$ and in the parton shower that we will discuss below. However, in the multi-parton matrix elements in \big{|}{\rho(\mu_{\textsc{h}}^{2})}\big{)}, the partons can have any momenta. Parton emissions with small transverse momenta need subtractions, but parton emissions with very large transverse momenta do not need subtractions. If we work at order $\alpha_{\mathrm{s}}^{k}$, then we can have up to $k$ high transverse momentum jets in addition to the Higgs boson in \big{|}{\rho(\mu_{\textsc{h}}^{2})}\big{)}.101010If we are limited to a $k=1$ shower, then Higgs plus one jet is LO in \big{|}{\rho(\mu_{\textsc{h}}^{2})}\big{)}. Then one might want define a lower cutoff on the $P_{\mathrm{T}}$ of the jet and use a calculation with $p+p\to H+J$ as the Born process, calculated at NLO. In that case, one has two calculations and one may want to define a procedure to merge them. If a $k=2$ shower is available, then \big{|}{\rho(\mu_{\textsc{h}}^{2})}\big{)} can include inclusive Higgs production at NNLO and Higgs plus one jet at NLO. Then merging different calculations is less needed.

VI.1 Moving the measurement operator

The first step towards a more general formulation is to interchange the order of the measurement operator $\mathcal{O}_{J}$ and the operators $\mathcal{D}^{-1}(\mu_{\textsc{h}}^{2})$ and $\mathcal{D}(\mu_{\textsc{h}}^{2})$. This doesn’t change the result, since $\mathcal{D}\,\mathcal{D}^{-1}=1$ and $\mathcal{O}_{J}$ commutes with 1. With the $\mathcal{O}_{J}$ moved, we have

[TABLE]

Here we have denoted

[TABLE]

where the product is expanded to $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$ using Eqs. (43), (44), (49), and (53):

[TABLE]

where

[TABLE]

These quantities \big{|}{\rho_{\textsc{h}}^{(n_{\textsc{r}},n_{\textsc{v}})}}\big{)} are finite without dimensional regularization.

VI.2 Introducing shower oriented parton distribution functions

We can do a little more by introducing an operator $\mathcal{F}(\mu_{\textsc{h}}^{2})$ that multiplies by parton distribution functions $f_{a/A}(\eta_{\mathrm{a}},\mu^{2})\,f_{b/B}(\eta_{\mathrm{b}},\mu^{2})$ and a luminosity factor. However, these are not the five-flavor $\overline{\mathrm{MS}}$ parton distribution functions that we used in \mathcal{F}_{\!\raisebox{-2.0pt}{{\overline{\phantom{\raisebox{-0.6pt}{$\scriptstyle{\textsc{ms}}$}}}\hskip-9.5pt\scriptstyle{\textsc{ms}}}}}(\mu^{2}). Rather, they are adapted to the choice of the definition for $\mathcal{D}(\mu^{2})$ that we use.

The shower oriented parton operator $\mathcal{F}(\mu^{2})$ is related to \mathcal{F}_{\!\raisebox{-2.0pt}{{\overline{\phantom{\raisebox{-0.6pt}{$\scriptstyle{\textsc{ms}}$}}}\hskip-9.5pt\scriptstyle{\textsc{ms}}}}}(\mu^{2}) by factor $\mathcal{K}(\mu^{2})$,

[TABLE]

This is a rather compact notation, so it is worthwhile to write it in more detail. The left hand side is defined by

[TABLE]

The right hand side is

[TABLE]

Here we take the kernel $K$ to be a product, so that each of the two parton distributions is transformed separately. We allow each kernel to depend on the momentum and flavor variables of the parton state to which $[\mathcal{F}(\mu^{2})\circ\mathcal{K}(\mu^{2})]$ is applied. The kernels each have a perturbative expansion beginning with

[TABLE]

The choice of $\mathcal{K}(\mu^{2})$ defines the shower-oriented parton distribution functions. The evolution of these parton distribution functions needs to be matched to the parton splitting functions introduced in the following sections. In particular, the choice of $\mathcal{K}(\mu^{2})$ is largely determined by the definition of the cutoff $\mu_{\mathrm{s}}^{2}$ that we use for the shower. We provide an example in Appendix B.

The parton operators \mathcal{F}_{\!\raisebox{-2.0pt}{{\overline{\phantom{\raisebox{-0.6pt}{$\scriptstyle{\textsc{ms}}$}}}\hskip-9.5pt\scriptstyle{\textsc{ms}}}}}(\mu^{2}) and $\mathcal{F}(\mu^{2})$ obey evolution equations

[TABLE]

Using Eq. (61), we see that the evolution kernels are related by

[TABLE]

With the transformation from \mathcal{F}_{\!\raisebox{-2.0pt}{{\overline{\phantom{\raisebox{-0.6pt}{$\scriptstyle{\textsc{ms}}$}}}\hskip-9.5pt\scriptstyle{\textsc{ms}}}}}(\mu^{2}) to $\mathcal{F}(\mu^{2})$, we write

[TABLE]

Thus our cross section is

[TABLE]

We will introduce $\mathcal{F}(\mu^{2})$ into another place in the formalism shortly.

VI.3 Changing the scale of the subtraction operators

Next, we would like to change the scale of the operators in Eq. (68) from a large scale $\mu_{\textsc{h}}^{2}$ to something smaller. We let $\mu_{\textsc{i}}^{2}$ be an “intermediate” size scale that is much smaller than the scale $Q^{2}[J]$ associated with the operator $\mathcal{O}_{J}$ but is nevertheless large compared to $1\ \mathrm{GeV}^{2}$ and is certainly large enough to allow the use of perturbation theory in $\alpha_{\mathrm{s}}(\mu_{\textsc{i}}^{2})$. We can change the parton factor to be evaluated at scale $\mu_{\textsc{i}}^{2}$ because this factor is a renormalization group invariant:

[TABLE]

The operator $\mathcal{D}(\mu^{2})$ is not invariant under changes of scale. However, we can write

[TABLE]

where

[TABLE]

Here we note that $\mathcal{D}(\mu_{1}^{2})$ generates $1/\epsilon$ poles and infrared singularities, but $\mathcal{D}^{-1}(\mu_{2}^{2})$ provides the proper subtractions to remove the poles and infrared singularities when we expand the product of operators to a fixed order of perturbation theory. Thus we can evaluate $\mathcal{U}_{\rm pert}(\mu_{2}^{2},\mu_{1}^{2})$ in four dimensions instead of $4-2\epsilon$ dimensions. The perturbative evolution operator $\mathcal{U}_{\rm pert}(\mu^{2},\mu^{\prime\,2})$ obeys the differential equation

[TABLE]

where

[TABLE]

Since $\mathcal{U}_{\rm pert}(\mu_{2}^{2},\mu_{1}^{2})$ is infrared finite, so is $S_{\rm pert}(\mu^{2})$. We can write the solution of Eq. (72) as

[TABLE]

where $\mathbb{T}$ indicates $\mu^{2}$ ordering of the exponential with smaller $\mu^{2}$ to the left. Working to order $\alpha_{\mathrm{s}}$ with use of Eqs. (47) and (48), we have

[TABLE]

where $\alpha_{\mathrm{s}}(\mu^{2})$ is the running coupling in the four dimensional theory and

[TABLE]

Recall that $\mu^{2}=\mu_{\mathrm{s}}^{2}$. Thus $\mathcal{S}^{(1,0)}_{\rm pert}(\mu^{2})$ is the derivative of (approximated) real emission graphs with respect to the ultraviolet cutoff that we impose. Similarly, $\mathcal{S}^{(0,1)}_{\rm pert}(\mu^{2})$ is the derivative of approximated one loop virtual graphs with respect to the ultraviolet cutoff. The subscript “pert” emphasizes that only perturbative Feynman diagrams are used to obtain $\mathcal{S}^{(1,0)}_{\rm pert}(\mu^{2})$ and $\mathcal{S}^{(0,1)}_{\rm pert}(\mu^{2})$.

With these changes, we have

[TABLE]

Now we note that very soft or collinear splittings at scales much smaller than $Q[J]^{2}$ are not resolved by the measurement operator $\mathcal{O}_{J}$. The operator $\mathcal{D}(\mu_{\textsc{i}}^{2})$ generates splittings at scales $\mu_{\textsc{i}}^{2}$ and smaller. Since we have chosen $\mu_{\textsc{i}}^{2}\ll Q[J]^{2}$, the operator $\mathcal{O}_{J}$ commutes with $\mathcal{D}(\mu_{\textsc{i}}^{2})$ to a good approximation, with an error of order $\mu_{\textsc{i}}^{2}/Q[J]^{2}$. Thus Eq. (68) can be written as

[TABLE]

VI.4 The inclusive infrared finite operator $\mathcal{V}(\mu^{2})$

We now introduce an operator $\mathcal{X}(\mu^{2})$ defined by

[TABLE]

Using $\mathcal{X}(\mu^{2})$, Eq. (78) is more compact:

[TABLE]

The operator $\mathcal{X}(\mu^{2})$ involves parton distribution functions and purely perturbative operators. If we evaluate the perturbative operators at order zero, we get simply $\mathcal{F}(\mu^{2})\mathcal{F}^{-1}(\mu^{2})$. Thus

[TABLE]

The operator $\mathcal{X}(\mu^{2})$, when expanded in powers of $\alpha_{\mathrm{s}}$, creates partons, up to $k$ partons at order $\alpha_{\mathrm{s}}^{k}$.

The operator $\mathcal{X}(\mu^{2})$ is infrared sensitive. When we apply $\mathcal{X}(\mu^{2})$ to a state \big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)}, we get a state \mathcal{X}(\mu^{2})\big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} containing poles $1/\epsilon$ and singularities when the partons that $\mathcal{X}(\mu^{2})$ creates become soft or collinear with other partons or with each other. In that sense, $\mathcal{X}(\mu^{2})$ is like $\mathcal{D}(\mu^{2})$. However, $\mathcal{X}(\mu^{2})$ contains the parton factor $\left[\mathcal{F}(\mu^{2})\circ\mathcal{K}(\mu^{2})\circ\mathcal{Z}_{F}(\mu^{2})\right]$. This factor gives $\mathcal{X}(\mu^{2})$ a property not shared by $\mathcal{D}(\mu^{2})$. If we integrate over the momenta of the partons created by $\mathcal{X}(\mu^{2})$ and sum over their colors and flavors by forming the inclusive sum \big{(}{1}\big{|}\mathcal{X}(\mu^{2})\big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)}, then the singularities cancel and we obtain a finite result.

In fact, we need to ensure that \big{(}{1}\big{|}\mathcal{X}(\mu^{2})\big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} is not only finite after dimensional regularization is removed but that it vanishes in the limit $\mu^{2}\to 0$. For example, if this quantity arises from an integration

[TABLE]

then, with subtractions included, $G(k^{2})$ needs to be a smooth function that is well enough behaved for $k^{2}\to 0$ that the integral is convergent. This property is needed later in Eq. (133).

Suppose for a moment that we worked in a modified theory, denoted by subscripts M, in which partons carried only momenta and flavors, but not color and spin. Then from the inclusive sum \big{(}{1}\big{|}\mathcal{X}_{\mathrm{M}}(\mu^{2})\big{|}{\{p,f\}_{m}}\big{)} we could define another operator $\mathcal{V}_{\mathrm{M}}(\mu^{2})$ that leaves the number of partons, their momenta and flavors unchanged:

[TABLE]

Then \lambda(\{p,f\}_{m},\mu^{2})=\big{(}{1}\big{|}\mathcal{V}_{\mathrm{M}}(\mu^{2})\big{|}{\{p,f\}_{m}}\big{)}. We define $\mathcal{V}_{\mathrm{M}}(\mu^{2})$ by Eq. (83) and

[TABLE]

Now return to QCD. With spin and color, we can define an operator $\mathcal{V}(\mu^{2})$ that satisfies \big{(}{1}\big{|}\mathcal{V}(\mu^{2})=\big{(}{1}\big{|}\mathcal{X}(\mu^{2}). However, its structure is more complex. The operator $\mathcal{X}(\mu^{2})$ can be expanded in powers of $\alpha_{\mathrm{s}}$:

[TABLE]

The order $n$ contribution, $\mathcal{X}^{(n)}(\mu^{2})$, is a sum of infrared sensitive operators,

[TABLE]

Acting on a state \big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} with $m$ final state partons, $\mathcal{X}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})$ produces a state with $m+n_{\textsc{r}}$ final state partons with momenta and flavors $\{\hat{p},\hat{f}\}_{m+n_{\textsc{r}}}$. There are integrations over the loop momenta $\{\ell_{1},\dots,\ell_{n_{\textsc{v}}}\}$ of $n_{\textsc{v}}$ virtual loops.

We need to understand the color and spin structure of $\mathcal{X}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})$. Suppose that we have constructed a basis of operators that act on the quantum spin$\otimes$color space and create a quantum spin$\otimes$color state for $n_{\textsc{r}}$ more partons. We label the basis operators by an index $i$. A convenient choice would be

[TABLE]

Then we could let

[TABLE]

Using these basis operators, we can expand $\mathcal{X}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu^{2})$ as

[TABLE]

Here $\mathcal{X}^{(n_{\textsc{r}},n_{\textsc{v}})}_{i,j}(\mu^{2})$ is still an operator on the momentum and flavor part of the statistical space, which has basis vectors \big{|}{\{p,f\}_{m}}\big{)}. In the case $n_{\textsc{r}}=0$, this operator adds no partons and leaves the parton momenta and flavors $\{p,f\}_{m}$ unchanged.

Now we wish to define another operator $\mathcal{V}(\mu^{2})$ with an expansion

[TABLE]

The order $n$ contribution is to add no partons and leave the parton momenta and flavors $\{p,f\}_{m}$ unchanged, but it still can be a non-trivial operator on the spin$\otimes$color space

[TABLE]

Here $\mathcal{V}^{(n)}_{i,j}(\mu^{2})$ is still an operator on the momentum and flavor part of the statistical space. We want $\mathcal{V}(\mu^{2})$ to be related to $\mathcal{X}(\mu^{2})$ by

[TABLE]

Thus we want

[TABLE]

This is the same as

[TABLE]

This needs to work for any choice of $m$-parton spin$\otimes$color states \big{\langle}{\{s^{\prime},c^{\prime}\}_{m}}\big{|} and \big{|}{\{s,c\}_{m}}\big{\rangle}, so we need an identity of spin$\otimes$color operators,

[TABLE]

The right hand side of Eq. (95) is an operator on the spin$\otimes$color space for $m$ final state partons. On the left hand side, the operators $\sigma_{i}^{(0)}$ form a basis for this space of operators, as do the operators $\sigma_{j}^{(0)\dagger}$, so the operators $\sigma_{j}^{(0)\dagger}\sigma_{i}^{(0)}$ span this space and are, in fact, over-complete. That is, one can always find coefficients \big{(}{1}\big{|}\mathcal{V}^{(n)}_{i,j}(\mu^{2})\ \big{|}{\{p,f\}_{m}}\big{)} so that we match the operator on the right hand side. However, the choice is not unique. At order $n=1$, we have made a simple choice in Deductor. It is beyond our scope here to investigate what choices might be best at NLO, $n=2$.

Since \big{(}{1}\big{|}\mathcal{X}(\mu^{2}) is infrared finite, Eq. (92) tells us that $\mathcal{V}(\mu^{2})$ is infrared finite.

VI.5 A more sophisticated shower evolution operator

Using Eq. (92), Eq. (80) becomes

[TABLE]

Since $\mathcal{V}(\mu_{\textsc{i}}^{2})\,\mathcal{F}(\mu_{\textsc{i}}^{2})$ does not change the number, momenta, or flavors of partons, it commutes with $\mathcal{O}_{J}$. Thus

[TABLE]

Now we can define the shower evolution operator that we need,

[TABLE]

With this definition, the cross section is

[TABLE]

This moves $\mathcal{F}$ next to \big{|}{\rho_{\textsc{h}}}\big{)} so that at the hard interaction we have the proper factors to make a cross section. It also moves $\mathcal{V}$ next to \big{|}{\rho_{\textsc{h}}}\big{)}. We will see later what the consequences of this are.

In Eq. (99), we use a scale $\mu_{\textsc{i}}^{2}$ that was left undefined except that it should be small compared to $Q^{2}[J]$ (which was the scale of $\mathcal{O}_{J}$) and should be large enough to allow the use of perturbation theory with coupling $\alpha_{\mathrm{s}}(\mu_{\textsc{i}}^{2})$. Our cross section is independent of the value of $\mu_{\textsc{i}}^{2}$. Let us now fix on a standard choice near the lower end of this range. We take $\mu_{\textsc{i}}^{2}\to\mu_{\mathrm{f}}^{2}$, where $\mu_{\mathrm{f}}^{2}$ is on the order of $1\ \mathrm{GeV}^{2}$. Then

[TABLE]

We can write $\mathcal{U}(\mu_{2}^{2},\mu_{1}^{2})$ in a simpler form. We note that, using Eqs. (71) and (79),

[TABLE]

Since the operator $\left[\mathcal{F}(\mu^{2})\circ\mathcal{K}(\mu^{2})\circ\mathcal{Z}_{F}(\mu^{2})\right]$ is independent of scale, this is

[TABLE]

VI.6 Probability preservation in $\mathcal{U}(\mu_{2}^{2},\mu_{1}^{2})$

The operator $\mathcal{U}(\mu_{2}^{2},\mu_{1}^{2})$ has an important property, which we now derive. From Eq. (102), we have

[TABLE]

Then using Eq. (92) twice, we have

[TABLE]

Multiplying any statistical state \big{|}{\rho}\big{)} by \big{(}{1}\big{|} gives the total probability associated with that state. Thus Eq. (104) says that shower evolution as represented in $\mathcal{U}(\mu_{2}^{2},\mu_{1}^{2})$ is probability preserving. Current parton shower algorithms are typically constructed to have this property. Here probability preservation is a derived property.

VI.7 Factorization

Let $Q^{2}[J]$ be the smallest scale associated with the measurement operator $\mathcal{O}_{J}$ in Eq. (100), as discussed in Sec. III.1. If $Q^{2}[J]$ us close to the scale $\mu_{\textsc{h}}^{2}$ of the hard process with which we start the shower, then measuring $\sigma[J]$ does not make use of the full power of a parton shower. Suppose now that $Q^{2}[J]\gg 1\ \mathrm{GeV}^{2}$ but that $Q^{2}[J]\ll\mu_{\textsc{h}}^{2}$. Then perturbation theory for $\sigma[J]$ is still applicable, but there may be large logarithms, $\log(\mu_{\textsc{h}}^{2}/Q^{2}[J])$ in the perturbative expansion of $\sigma[J]$. In many cases, a parton shower is useful for summing such logarithms. It takes a dedicated analysis to show that a given parton shower algorithm does sum the logarithms associated with a given operator $\mathcal{O}_{J}$, but there is at least a chance that if we use a parton shower we will do better than if we simply use fixed order perturbation theory. Thus we consider this sort of measurement operator and examine how factorization works when $Q^{2}[J]\ll\mu_{\textsc{h}}^{2}$.

We argued in Sec. III.1 that shower splittings at scale $\mu^{2}$ can change the measurement by a fraction $\mu^{2}/Q^{2}[J]$. We can neglect these modifications as long as $\mu^{2}$ is small enough, say

[TABLE]

We want $\epsilon_{\mathrm{s}}$ to be small enough that we can regard fractional errors of order $\epsilon_{\mathrm{s}}$ as negligible. However, we may want $\log(1/\epsilon_{\mathrm{s}})$ not to be large.

We can use our knowledge of the scale of $\mathcal{O}_{J}$ by writing

[TABLE]

Then writing

[TABLE]

results in a negligible error. With this substitution, Eq. (100) becomes

[TABLE]

Now, factorization for the cross section measured by $\mathcal{O}_{J}$ requires that splittings at scales smaller than $\epsilon_{\mathrm{s}}\,Q^{2}[J]$ not affect the cross section. Thus we need

[TABLE]

This follows by using Eq. (104) to obtain \big{(}{1}\big{|}\mathcal{U}(\mu_{\mathrm{f}}^{2},\epsilon_{\mathrm{s}}\,Q^{2}[J])=\big{(}{1}\big{|}.

We should emphasize that in order to measure the cross section corresponding to the infrared safe operator $\mathcal{O}_{J}$ with scale $Q^{2}[J]$, it is not necessary to cut off the shower at scale $\epsilon_{\mathrm{s}}\,Q^{2}[J]$ as in Eq. (109). Rather, one simply runs the shower down to $\mu_{\mathrm{f}}^{2}$ and measures $\mathcal{O}_{J}$ on the final state produced by the full shower, as in Eq. (100). When we do that, we are setting $\epsilon_{\mathrm{s}}=\mu_{\mathrm{f}}^{2}/Q^{2}[J]$, so the error estimate $\mathcal{O}(\epsilon_{\mathrm{s}})$ becomes $\mathcal{O}(\mu_{\mathrm{f}}^{2}/Q^{2}[J])$.

VI.8 The shower evolution equation

Using its definition Eq. (98), we see that the shower evolution operator $\mathcal{U}(\mu^{2},\mu^{\prime\,2})$ obeys an evolution equation of the form

[TABLE]

Thus

[TABLE]

Since, according to Eq. (104), \big{(}{1}\big{|}\mathcal{U}(\mu^{2},\mu^{\prime\,2})=\big{(}{1}\big{|}, we have

[TABLE]

Using Eqs. (98) and (72), we see that the shower generator $\mathcal{S}$ in Eq. (110) is

[TABLE]

It is convenient to define

[TABLE]

Also, we can use Eq. (65) for the evolution of $\mathcal{F}(\mu^{2})$ and we can note that since $\mathcal{V}(\mu^{2})$ does not change the number of partons or their momenta or flavors, $\mathcal{V}(\mu^{2})$ commutes with $\mathcal{F}(\mu^{2})$. Then

[TABLE]

The operator $\mathcal{V}(\mu^{2})$ here has a perturbative expansion beginning with

[TABLE]

Then also

[TABLE]

A sensible procedure for determining $\mathcal{S}(\mu^{2})$ is to expand it perturbatively to whatever order is known, e.g. order $\alpha_{\mathrm{s}}^{k}$,

[TABLE]

It is of interest to see how this works out at order $\alpha_{\mathrm{s}}$. Since $\mathcal{S}_{\rm pert}(\mu^{2})$ and $\mathcal{S}_{\mathcal{V}}$ are already order $\alpha_{\mathrm{s}}$, we can simply replace $\mathcal{V}(\mu^{2})$ by 1 in the first and third terms of Eq. (115). Then we have

[TABLE]

Let us multiply by \big{(}{1}\big{|} and use \big{(}{1}\big{|}\mathcal{S}^{(1)}(\mu^{2})=0. In Ref. NSThreshold , we found that \big{(}{1}\big{|}\mathcal{S}_{\rm pert}^{(1,0)}(\mu^{2}) has a simple form,

[TABLE]

where $\bar{\mathcal{S}}_{\rm pert}^{(1,0)}(\mu^{2})$ leaves the number of partons, their momenta and their flavors unchanged but has a non-trivial color structure.111111The operator $[\mathcal{F}(\mu^{2})\circ\bar{\mathcal{S}}^{(1,0)}(\mu^{2})]\,\mathcal{F}^{-1}(\mu^{2})$ was called $\mathcal{V}$ in Ref. NSThreshold , but here we are letting $\mathcal{V}$ denote a different operator. We also divide $\mathcal{S}_{\rm pert}^{(0,1)}(\mu^{2})$ into two pieces

[TABLE]

Here $\mathcal{S}_{{\mathrm{i}}\pi}^{(0,1)}(\mu^{2})$ is the contribution from the virtual graphs that is proportional to ${\mathrm{i}}\pi$, while $\mathcal{S}_{\rm Re}^{(0,1)}(\mu^{2})$ is the remaining part. We note that NSThreshold

[TABLE]

This gives us

[TABLE]

All of the operators here leave the parton state unchanged except for being operators on the color and spin space. We define $\mathcal{S}_{\mathcal{V}}^{(1)}(\mu^{2})$ to have the color and spin structure of the right hand side of the equation, so that

[TABLE]

This result also gives us121212In Ref. NSThreshold , we neglected the ${\mathrm{i}}\pi$ term and we averaged over spin. Then the right hand side of Eq. (125) was called $\mathcal{H}_{I}-\mathcal{V}$.

[TABLE]

Compare this to the more general Eq. (219) in Appendix D.

To use $\mathcal{S}^{(1)}(\mu^{2})$, one solves Eq. (110) in the form

[TABLE]

where

[TABLE]

Here the Sudakov factor $\mathcal{N}$ is the exponential of the part of $\mathcal{S}$ that does not change the number of partons or their momenta or flavors. Normally its spin and color structure is simplified and the ${\mathrm{i}}\pi$ contribution is not included. The splitting operator $\mathcal{F}(\mu^{2})\,\mathcal{S}_{\rm pert}^{(1,0)}(\mu^{2})\,\mathcal{F}^{-1}(\mu^{2})$ adds one parton. Its spin and color structure is also normally simplified. Then Eq. (126) is implemented by solving it iteratively, so that there are some number of splittings interleaved with Sudakov factors.

VI.9 The structure of $\mathcal{V}(\mu^{2})$

In Eq. (100), there is a factor $\mathcal{V}(\mu_{\textsc{h}}^{2})$ that multiplies $\mathcal{F}(\mu_{\textsc{h}}^{2})$ and \big{|}{\rho_{\textsc{h}}}\big{)}. We can write this operator as

[TABLE]

where

[TABLE]

The operator $\mathcal{U}_{\mathcal{V}}(\mu^{2},\mu^{\prime\,2})$ obeys the evolution equation

[TABLE]

where $\mathcal{S}_{\mathcal{V}}(\mu^{2})$ was defined in Eq. (114). Thus

[TABLE]

To use Eq. (131), we apply $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ to a statistical state \big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} that contributes to \big{|}{\rho_{\textsc{h}}}\big{)}. We expand $\mathcal{S}_{\mathcal{V}}(\mu^{2})$ in powers of $\alpha_{\mathrm{s}}(\mu^{2})$ with factors involving the running parton distributions $\mathcal{F}(\mu^{2})$. The available scales other than $\mu^{2}$ come from \big{|}{\rho_{\textsc{h}}}\big{)}, so the relevant matrix elements involve $\mu^{2}/\mu_{\textsc{h}}^{2}$. Since $\mathcal{V}(\mu^{2})$ is an infrared finite operator, the perturbative coefficients in $\mathcal{S}_{\mathcal{V}}(\mu^{2})$ will then be proportional to $\mu^{2}/\mu_{\textsc{h}}^{2}$, possibly times logarithms of $\mu^{2}/\mu_{\textsc{h}}^{2}$. Thus the low $\mu^{2}$ end of the integration is power suppressed. However, it is important that there is a lower bound $\mu_{\mathrm{f}}^{2}$ on the integration. That is because the running coupling $\alpha_{\mathrm{s}}(\mu^{2})$ and the running parton distributions $\mathcal{F}(\mu^{2})$ are not well defined for very small $\mu^{2}$.

We are left with a factor

[TABLE]

We can expand this perturbatively in powers of $\alpha_{\mathrm{s}}(\mu_{\mathrm{f}}^{2})$, using parton distributions $\mathcal{F}(\mu_{\mathrm{f}}^{2})$. The coefficients of $\alpha_{\mathrm{s}}^{n}(\mu_{\mathrm{f}}^{2})$ for $n\geq 1$ are then proportional to $\mu_{\mathrm{f}}^{2}/\mu_{\textsc{h}}^{2}$, possibly times logarithms of $\mu_{\mathrm{f}}^{2}/\mu_{\textsc{h}}^{2}$. Since $\mu_{\mathrm{f}}^{2}\ll\mu_{\textsc{h}}^{2}$, we can safely neglect all of the higher order terms and simply replace

[TABLE]

Thus we make the replacement $\mathcal{V}(\mu_{\textsc{h}}^{2})\to\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ in Eq. (100), giving us

[TABLE]

The factor $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ does two things. First, it provides perturbative corrections to the hard scattering state \big{|}{\rho_{\textsc{h}}}\big{)}, which we need in order to calculate the cross section correct to $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$. For this purpose, it would suffice to expand the exponential in $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ to the desired perturbative order. The second function of $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ is to sum threshold logarithms. For this purpose, it is important that $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ is an exponential.

To understand the relation of the operator $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ to threshold logarithms, it is instructive to look at it at order $\alpha_{\mathrm{s}}$ with $\Lambda$ ordering for the shower. It is structurally the same as the operator introduced in Ref. NSThreshold , which concerned the summation of threshold logarithms.131313In Ref. NSThreshold , an approximate version of $\mathcal{U}_{\mathcal{V}}(\mu_{2}^{2},\mu_{1}^{2})$ was used to sum threshold logarithms, but it appeared between each pair of parton splittings at scales $\mu_{1}^{2}$ and $\mu_{2}^{2}$ and also between the initial hard interaction and the first splitting. This caused problems, which were alleviated by inserting an artificial cut in $\mathcal{S}_{\mathcal{V}}(\mu^{2})$. The analysis in Ref. NSThreshold simply averaged over spins, so we leave out spin here. The operator $\mathcal{S}_{\mathcal{V}}(\mu^{2})$ in Ref. NSThreshold contains several terms. Rather than listing them all, we simply recall the most important terms:

[TABLE]

where $C_{a}=C_{\mathrm{F}}$ when $a$ is a quark flavor and $C_{\mathrm{g}}=C_{\mathrm{A}}$. The operator $[1\otimes 1]$ is the unit operator in the color space. Note that when $\mu^{2}\ll\mu_{\textsc{h}}^{2}$, the range of the $z$ integration shrinks to just a tiny region near $z=1$. Thus, this contribution might seem unimportant. However, in Eq. (135) there is a factor $1/(1-z)$. This multiplies a factor involving the parton distribution functions. The result is that $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ is substantially different from 1 when the parton distribution functions are falling quickly as the momentum fraction grows. This gives a “threshold logarithm” effect that we can sum in a leading approximation by using $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$.

The form of $\mathcal{S}_{\mathcal{V}}(\mu^{2})$ depends on the definition of the infrared sensitive operator $\mathcal{D}(\mu^{2})$. In particular, shower evolution uses parton distributions $\mathcal{F}(\mu^{2})$ that are related to the $\overline{\mathrm{MS}}$ parton distributions according to Eq. (61). This relation is discussed in Appendix B. Briefly, $\overline{\mathrm{MS}}$ parton distributions are defined by removing ultraviolet divergences using $\overline{\mathrm{MS}}$ renormalization, while the ultraviolet region of shower splittings is removed with a cutoff at scale $\mu_{\mathrm{s}}^{2}$ as defined in $\mathcal{D}(\mu^{2})$. See Sec. V.3. The evolution kernel $\mathcal{P}(\mu^{2})$ for the parton distributions reflects the definition of $\mathcal{F}(\mu^{2})$. This evolution kernel appears in $\mathcal{S}_{\mathcal{V}}(\mu^{2})$, Eqs. (124) and (218). Thus the choice of shower oriented parton distribution functions affects how the summation of threshold logarithms appears in the overall result (134) for the cross section: part of the summation of threshold logarithms appears as $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ and part appears as the redefinition of the parton distributions, $\mathcal{F}(\mu^{2})$ in Eq. (134). This is described in Sec. (9.4) of Ref. NSThreshold . If one uses $k_{\mathrm{T}}$ ordering instead of $\Lambda$ ordering for the shower, then the shower oriented parton distributions are very close to the $\overline{\mathrm{MS}}$ parton distributions and all of the leading threshold logarithms appear in $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$. On the other hand, in the factorization scheme in ref. KrkNLO2015 ; KrkNLO2016 , a different definition of factorization results in the leading threshold corrections being absorbed into the parton distribution functions.

VII Summary and outlook

We began with an expression (42) for the cross section $\sigma[J]$ for an infrared safe measurement $J$ calculated at $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$. The pieces in this expression are infrared divergent in four dimensions, so that they are defined by working in $4-2\epsilon$ dimensions. Integrating over the phase space that is unresolved by the measurement leads to some cancellations of poles $1/\epsilon$. Other poles cancel after factorization of initial state infrared sensitivity into parton distribution functions. This leaves a result that is finite in four dimensions, even though it consists of pieces that are divergent in four dimensions. We introduced an infrared sensitive operator $\mathcal{D}(\mu^{2})$ and an inclusive infrared finite operator $\mathcal{V}(\mu^{2})$ to help organize the cancellations.

After some analysis, we have represented the cross section as Eq. (134). Here the separate factors are all finite in four dimensions. If we expand this expression to order $\alpha_{\mathrm{s}}^{k}$, we have the same cross section that we started with except for a power suppressed contribution that we have dropped.

In Eq. (134), we have a hard scattering state \big{|}{\rho_{\textsc{h}}}\big{)} and a factor $\mathcal{F}(\mu_{\textsc{h}}^{2})$ that supplies parton distribution functions and a parton luminosity factor so that if we trace over colors and spins, we have a differential cross section in the space of parton number, flavors, and momenta. Then we have a factor $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ that supplies a summation of threshold logarithms associated with the hard state and also part of the $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$ perturbative corrections to the hard scattering cross section. Next, we have a complete parton shower generated by $\mathcal{U}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$. The parton shower operator preserves inclusive probabilities: \big{(}{1}\big{|}\mathcal{U}(\mu_{2}^{2},\mu_{1}^{2})=\big{(}{1}\big{|}. We have ended the shower at a scale $\mu_{\mathrm{f}}^{2}$. After that, the factor \big{(}{1}\big{|}\mathcal{O}_{J} represents the measurement that we want to make. We suppose that this is an infrared safe measurement that is not sensitive to soft or collinear parton splittings at the scale $\mu_{\mathrm{f}}^{2}$ or below. That means that the error in the calculation, estimated by $\mathcal{O}(\mu_{\mathrm{f}}^{2}/Q[J]^{2})$, is small. With such an infrared safe measurement, the result of the measurement is not sensitive to hadronization. If we wanted to use a measurement operator that is sensitive to hadronization, then we would need to include a model of hadronization before the measurement operator. It is then less certain what a good error estimate is.

There is some temptation to imagine Eq. (134) as being simpler than it is. In our Higgs boson example, if we expand \big{(}{1}\big{|}\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})\,\mathcal{F}(\mu_{\textsc{h}}^{2})\big{|}{\rho_{\textsc{h}}}\big{)} to order $\alpha_{\mathrm{s}}^{k}$, it is the $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$ inclusive cross section to make a Higgs boson. The operator $\mathcal{U}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ generates a probability preserving parton shower. Thus it might seem that one takes the hard scattering cross section and then distributes the probability across different final states according to what the shower generates. However, \mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})\,\mathcal{F}(\mu_{\textsc{h}}^{2})\big{|}{\rho_{\textsc{h}}}\big{)} is not a cross section. It is a statistical state, representing different numbers of final state partons, which come with their own quantum color and spin states. The shower operator acts separately on each component of this statistical state. Then if we measure $\sigma[J]$ for an observable that is more complicated than just the inclusive measurement of a Higgs boson, the separate contributions are not sensible by themselves, but they sum to give $\sigma[J]$ correct to order $\alpha_{\mathrm{s}}^{k}$ with only a power suppressed correction.

We have spoken of getting $\sigma[J]$ correct to order $\alpha_{\mathrm{s}}^{k}$, but, of course, that is not the point of a parton shower. In applying Eq. (134), one would evaluate the splitting function $\mathcal{S}(\mu^{2})$ in the exponent of $\mathcal{U}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ to order $\alpha_{\mathrm{s}}^{k}$, then retain $\mathcal{U}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ as an exponential. When the desired measurement operator $\mathcal{O}_{J}$ contains widely different scales, the cross section will contain large logarithms. Then $\mathcal{U}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ has the potential to sum these logarithms. After all, it exponentiates the soft and collinear singularities of QCD at order $\alpha_{\mathrm{s}}^{k}$. Unfortunately, one needs to study the structure of $\mathcal{U}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ as it relates to the structure of $\mathcal{O}_{J}$ in order to check how well the shower does in summing the logarithms.

One can wonder whether the formalism of this paper is of any use for just a LO shower. We suggest that it is. If one averages over spins and makes the leading color approximation, the shower operator $\mathcal{U}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ generates a rather conventional probability preserving dipole shower. With $\Lambda^{2}$ as the ordering variable, it is the leading color version of the shower in Deductor. One can choose other ordering variables. The operator $\mathcal{U}_{\mathcal{V}}(\mu_{\mathrm{f}}^{2},\mu_{\textsc{h}}^{2})$ generates threshold corrections, as described in Ref. NSThreshold . These corrections are numerically important in some cases and could be included in standard parton shower programs.

In fact, Eq. (134) has been useful in improving Deductor. While working on our paper NSThreshold on threshold corrections, we did not have Eq. (134). The result was a structure that had certain undesirable features that needed to be controlled by means of an ad hoc cutoff. The current more general formulation in Eq. (134) removes this problem, although it does not much change the numerical results. We present phenomenological results from the new version of Deductor in a separate paper NSThresholdII .

The formalism is based on an operator $\mathcal{D}(\mu^{2})$ that encodes the infrared structure of QCD starting with a state with any number $m$ of final state partons. If we know $\mathcal{D}(\mu^{2})$ up to order $\alpha_{\mathrm{s}}^{k}$, then we generate, in an automatic way, the subtraction terms for an $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$ perturbative calculation. See the example of the toy model in Appendix A. This appears to us to be simpler than constructing the subtraction terms directly Daleo:2009yj ; Boughezal:2010mc ; Gehrmann:2011wi ; GehrmannDeRidder:2012ja ; Currie:2013vh ; Somogyi:2006da ; Somogyi:2006db ; Somogyi:2013yk . From $\mathcal{D}(\mu^{2})$ at order $\alpha_{\mathrm{s}}^{k}$, we also generate, in an automatic way, the shower splitting kernels at order $\alpha_{\mathrm{s}}^{k}$.

The perturbative contributions to $\mathcal{D}(\mu^{2})$ are not simple in full QCD. Furthermore, their form depends on choices like the momentum mapping scheme and the choice of a hardness ordering variable. At order $\alpha_{\mathrm{s}}^{1}$, we have made the required choices, made suitable approximations, and calculated the corresponding splitting functions $\mathcal{S}(\mu^{2})$ in Refs. NSThreshold ; NSThresholdII . Similarly, from $\mathcal{D}(\mu^{2})$ we generate the inclusive infrared finite operator $\mathcal{V}(\mu^{2})$. In general, there are some choices that one can make in defining $\mathcal{V}(\mu^{2})$. At order $\alpha_{\mathrm{s}}^{1}$, we have made the required choices, made suitable approximations, and calculated $\mathcal{S}_{\mathcal{V}}(\mu^{2})$ in Ref. NSThreshold ; NSThresholdII .

In the first order Deductor shower, we need only $d\mathcal{D}^{(1)}(\mu^{2})/d\log(\mu^{2})$, which is finite in four dimensions. For subtractions and matching to the shower, one needs the full integrated operator $\mathcal{D}^{(1)}(\mu^{2})$. Techniques for this are described in Ref. Czakon and the earlier papers Robens1 ; Robens2 ; Czakon1 .

We leave it to future work to make suitable choices for a momentum mapping scheme, a hardness ordering variable, and definitions away from the strict soft and collinear limits so as to construct order $\alpha_{\mathrm{s}}^{2}$ contributions to $\mathcal{D}(\mu^{2})$. With a choice for color structure, we could then also construct $\mathcal{V}(\mu^{2})$. We thus hope that the formalism presented in this paper might prove useful in developing a parton shower with order $\alpha_{\mathrm{s}}^{2}$ splitting functions.

We also hope that the formalism presented in this paper might provide support for the view that a parton shower is similar to a more straightforward perturbative calculation at $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$. In this view, the parton shower is an approximate way to calculate cross sections, but the approximation is systematically improvable by working at higher perturbative order. In a practical program, there may be further approximations with respect to color and spin. These need a separate justification.

Acknowledgements.

This work was supported in part by the United States Department of Energy under grant DE-SC0011640. This project was begun while the authors were at the Munich Institute for Astro- and Particle Physics program “Challenges, Innovations and Developments in Precision Calculations for the LHC.” It was continued while the authors were at the Kavli Institute for Theoretical Physics at the University of California, Santa Barbara, program “LHC Run II and the Precision Frontier,” which was supported by the U. S. National Science Foundation under Grant No. NSF PHY11-25915. DS thanks the Erwin Schrödinger Institute of the University of Vienna for its hospitality at its program “Challenges and Concepts for Field Theory and Applications in the Era of the LHC Run-2.” We thank the MIAPP, the KITP, and the ESI for providing stimulating research environments. We thank Tim Cohen, Marat Freytsis, and Jonathan Gaunt for helpful discussions.

Appendix A A toy model for parton shower operators

The construction in this paper relies on an operator $\mathcal{D}$ that contains the infrared singularities of statistical states $|\rho)$ that represent the color and spin density matrix elements of QCD. The operators $\mathcal{D}$ can be used to construct both the infrared subtractions needed for a perturbative calculation of the cross section at $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$ and also the splitting kernels needed to construct a parton shower at the corresponding order, $\mathrm{N}^{k-1}\mathrm{L}\mathrm{O}$. This construction has been quite abstract, especially since we lack an example of $\mathcal{D}$ at $\mathrm{N}^{2}\mathrm{L}\mathrm{O}$, which would correspond to an $\mathrm{N}\mathrm{L}\mathrm{O}$ shower.

In this appendix, we illustrate some of the ideas of the paper using a toy model that provides a concrete example of $|\rho)$ and $\mathcal{D}$ at $\mathrm{N}^{2}\mathrm{L}\mathrm{O}$. In this example, we construct the splitting kernel $\mathcal{S}_{\rm pert}$ at $\mathrm{N}\mathrm{L}\mathrm{O}$. The toy model is very simple. There are no parton distributions. The coupling $\alpha_{\mathrm{s}}$ does not run. The momenta are one dimensional. There is no spin. There is quantum color, but the color structure is vastly simplified compared to what one has in real QCD.

A.1 Statistical states in the toy model

We use momentum states $\{p_{1},p_{2},\cdots,p_{m}\}$ for $m$ partons, with each $p_{i}$ being a real number in the range $0<p_{i}<\infty$. We use “color” represented by basis states labeled by an index pair $(c_{\textsc{r}},c_{\textsc{v}})$ with $c_{\textsc{r}}\in\{0,1,\dots,m\}$, $c_{\textsc{v}}\in\{0,1,\dots\}$. The statistical states have the form

[TABLE]

These are defined to be invariant under permutations of the $p_{i}$. The Born level cross section is \big{|}{\{\}_{0},(0,0)}\big{)}.

We will make use of color operators $\mathcal{C}_{\mathrm{R}}$ and $\mathcal{C}_{\mathrm{V}}$ that act on the space of statistical states. These have the form

[TABLE]

where $\mathcal{A}^{\dagger}_{\mathrm{R}}$ and $\mathcal{A}^{\dagger}_{\mathrm{R}}$ are raising operators:

[TABLE]

The inclusive probabilities corresponding to our statistical states are defined by

[TABLE]

Here $N_{\mathrm{c}}=3$ represents the number of colors. Thus keeping only $c_{\textsc{r}}=c_{\textsc{v}}=0$ is analogous to the leading color approximation.

We include an operator representing an observable. The observable depends only on the parton momenta, not on their color state:

[TABLE]

A.2 Perturbative hard scattering states

We take the perturbative hard scattering states in the toy model to have the form

[TABLE]

Here $|\rho^{(n_{\textsc{r}},n_{\textsc{v}})})$ represents a perturbative contribution with $n_{\textsc{r}}$ real partons emitted and $n_{\textsc{v}}$ virtual loops. The contributions are defined to be

[TABLE]

The integrals here are regularized in the infrared by factors $(p^{2}/Q^{2})^{\epsilon}$. The individual contributions $|\rho^{(n_{\textsc{r}},n_{\textsc{v}})})$ contain singularities when parton momenta become small and $1/\epsilon$ poles that arise from integration from virtual parton momenta. Furthermore, the contributions with real parton emissions have different color states than the corresponding contributions with virtual loops. However, when we calculate the cross section $(1|\mathcal{O}_{J}|\rho)$, we can use the fact that

[TABLE]

so that the color contributions from real emissions and virtual loops match. Then, in fact, there are $\mathit{real}-\mathit{virtual}$ cancellations, with the result that $(1|\mathcal{O}_{J}|\rho)$ is infrared finite.

A.3 The infrared sensitive operator

The infrared structure of this is fairly simple and can be represented using the infrared sensitive operator $\mathcal{D}(\mu_{\mathrm{s}}^{2})$ with

[TABLE]

with

[TABLE]

From $\mathcal{D}(\mu_{\mathrm{s}}^{2})$ we can construct $\mathcal{D}^{-1}(\mu_{\mathrm{s}}^{2})$,

[TABLE]

A.4 Subtractions for the hard scattering states

We can now construct the subtracted statistical state including the measurement operator,

[TABLE]

This has the expansion

[TABLE]

At first order in $\alpha_{\mathrm{s}}$, there are two terms. The first is

[TABLE]

This is

[TABLE]

The first term is singular when $p_{1}^{2}\to 0$, but the subtraction from $\mathcal{D}^{(1,0)}(\mu_{\mathrm{s}}^{2})$ eliminates the singularity. Then $(1|\hat{\rho}^{(1,0)})$ is finite at $\epsilon=0$. The coefficient of $\alpha_{\mathrm{s}}$ that corresponds to virtual graphs is

[TABLE]

This is

[TABLE]

The first term has a $1/\epsilon$ pole from $k_{1}^{2}\to 0$, but the subtraction from $\mathcal{D}^{(0,1)}(\mu_{\mathrm{s}}^{2})$ eliminates the $k_{1}^{2}\to 0$ singularity. Then $|\hat{\rho}^{(0,1)})$ is finite at $\epsilon=0$.

At order $\alpha_{\mathrm{s}}^{2}$, there are three terms in $|\hat{\rho})$. Let us look at the contribution from two real emissions:

[TABLE]

This is

[TABLE]

Each term here exhibits infrared singularities, but the singularities cancel. Then $(1|\hat{\rho}^{(2,0)})$ can be evaluated at $\epsilon=0$. Specifically, for $p_{2}^{2}\to 0$ at fixed $p_{1}^{2}$, the first term cancels the third term and the second term cancels the fourth term. For $p_{1}^{2}\to 0$ with fixed $p_{2}^{2}$, the first two terms are nonsingular, while the third term cancels the fourth term. When $p_{2}^{2}\to 0$ and $p_{1}^{2}\to 0$, the first term cancels the second term and the third term cancels the fourth term.

The contribution to $|\hat{\rho})$ from two virtual emissions is similar. We have

[TABLE]

This is

[TABLE]

Each term here exhibits $1/\epsilon$ poles, but the poles cancel. The pattern of cancellations is the same as for $|\hat{\rho}^{(2,0)})$.

The contribution to $|\hat{\rho})$ from one real emission and one virtual emission is a little more complicated. We have

[TABLE]

We obtain

[TABLE]

Each term here exhibits $1/\epsilon$ poles and singularities but the poles and singularities cancel. Then $|\hat{\rho}^{(1,1)})$ can be evaluated at $\epsilon=0$. Specifically, expanding the square brackets gives us seven terms, $T_{i}$, $i\in\{1,\dots,7\}$. For the singularity in the integrand at $k_{1}^{2}\to 0$ with fixed $p_{1}^{2}$, $T_{2}$ cancels $T_{7}$ and $T_{6}$ cancels $T_{4}+T_{5}$. For the singularity at $p_{1}^{2}\to 0$ with fixed $k_{1}^{2}$, $T_{1}$ cancels $T_{5}$ and $T_{6}$ cancels $T_{3}+T_{7}$. For the singularity when $p_{1}^{2}\to 0$ and $k_{1}^{2}\to 0$, $T_{1}$ cancels $T_{3}$, $T_{2}$ cancels $T_{4}$, and $T_{6}$ cancels $T_{5}$ + $T_{7}$.

We thus see that $\mathcal{D}^{-1}(\mu_{s}^{2})$ provides subtraction terms that cancel all of the singularities of $|\rho)$ at order $\alpha_{\mathrm{s}}$ and $\alpha_{\mathrm{s}}^{2}$. If we had wanted to construct all of the subtraction terms directly, it would have been somewhat difficult. However, constructing $\mathcal{D}(\mu_{s}^{2})$ for our toy model was quite simple. Then obtaining the subtraction terms was automatic.

We also note that the perturbative states $|\rho^{(n_{\textsc{r}},n_{\textsc{v}})})$ have non-trivial color structures in our toy model. Then the operators $\mathcal{D}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu_{\mathrm{s}}^{2})$ must reflect this color structure. If they do not, the subtractions will not work.

A.5 The perturbative shower evolution operator

We are now in a position to construct the order $\alpha_{\mathrm{s}}$ and $\alpha_{\mathrm{s}}^{2}$ terms in the perturbative shower generating operator, Eq. (73),

[TABLE]

The operator $\mathcal{S}_{\rm pert}$ will have the expansion

[TABLE]

where

[TABLE]

To use this, we first need the derivatives of $\mathcal{D}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu_{\mathrm{s}}^{2})$.

The derivatives of $\mathcal{D}^{(1,0)}(\mu_{\mathrm{s}}^{2})$ and $\mathcal{D}^{(1,0)}(\mu_{\mathrm{s}}^{2})$ are simple

[TABLE]

Here in the state $|\{p_{1},\dots,p_{m},p_{m+1}\},(c_{\textsc{r}},c_{\textsc{v}}))$, we have substituted $\mu_{\mathrm{s}}$ for the real number $p_{m+1}$.

The derivative of $\mathcal{D}^{(2,0)}(\mu_{\mathrm{s}}^{2})$ is a little more complicated. There are two terms. If we rename the integration variable and permute the arguments of the statistical state to match in the two terms, we obtain

[TABLE]

The derivative of $\mathcal{D}^{(0,2)}(\mu_{\mathrm{s}}^{2})$ is similar to the derivative of $\mathcal{D}^{(2,0)}(\mu_{\mathrm{s}}^{2})$:

[TABLE]

The singularity at $k_{1}^{2}\to 0$ gives us a $1/\epsilon$ pole after integration.

The derivative of $\mathcal{D}^{(1,1)}(\mu_{\mathrm{s}}^{2})$ is more complicated. After combining terms, we obtain

[TABLE]

In the first term, there is a pole from the integration region $k_{1}^{2}\to 0$, while in the second term there is a singularity at $p_{m+1}^{2}\to 0$.

With the derivatives of $\mathcal{D}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu_{\mathrm{s}}^{2})$ at hand, it is straightforward to use Eq. (161) to construct the shower kernel $\mathcal{S}_{\rm pert}(\mu_{\mathrm{s}}^{2})$.

The shower kernel at order $\alpha_{\mathrm{s}}$ is simple:

[TABLE]

For the shower kernel for two real emissions, we find

[TABLE]

The subtraction removes the infrared singularity at $p_{m+1}^{2}\to 0$, allowing us to set $\epsilon\to 0$. This gives

[TABLE]

For the shower kernel for two virtual loops, we find a result that is analogous to what we found above. The subtraction removes the infrared singularity at $k_{1}^{2}\to 0$ so that we can set $\epsilon\to 0$, giving

[TABLE]

For the shower kernel with one real emission and one virtual loop, we obtain

[TABLE]

Again, the subtraction removes the infrared singularities. Then we can set $\epsilon\to 0$, giving

[TABLE]

Thus Eq. (161) gives us a completely straightforward way to calculate the operators $\mathcal{S}_{\rm pert}^{(n_{\textsc{r}},n_{\textsc{v}})}(\mu_{\mathrm{s}}^{2})$. All of these operators at order $\alpha_{\mathrm{s}}^{1}$ and $\alpha_{\mathrm{s}}^{2}$ are infrared finite.

A.6 The operator $\mathcal{V}(\mu_{\mathrm{s}}^{2})$

We also construct an operator $\mathcal{X}(\mu_{\mathrm{s}}^{2})$ in Eq. (79). In our toy model, which has no parton distributions, we have

[TABLE]

Then we define an operator $\mathcal{V}(\mu_{\mathrm{s}}^{2})$ using Eq. (92),

[TABLE]

In our toy model, at least up to the order that we have defined it, the relation $(1|\mathcal{C}_{\mathrm{R}}=(1|\mathcal{C}_{\mathrm{V}}$ gives us

[TABLE]

Thus

[TABLE]

Then according to the definition Eq. (98) (with $\mathcal{F}=1$) we have

[TABLE]

The full shower operator $\mathcal{U}$ is the same as $\mathcal{U}_{\rm pert}$.

Appendix B Transformation to shower oriented parton distribution functions

In Sec. VI.2, we introduced shower oriented parton distribution functions that are adapted to the choice of the cutoff $\mu_{\mathrm{s}}^{2}$ used in the shower. The shower oriented parton distributions $f_{a/A}(\eta_{\mathrm{a}},\mu^{2})$ are related to the five-flavor $\overline{\mathrm{MS}}$ parton distribution functions f_{a/A}^{{\overline{\phantom{\raisebox{-0.6pt}{\scriptstyle{\textsc{ms}}}}}\hskip-9.5pt\scriptstyle{\textsc{ms}}}}(\eta_{\mathrm{a}},\mu^{2}) by means of a kernel $K$. Following the notation of Sec. VI.2, the transformation of the parton distribution for hadron A is

[TABLE]

The kernel $K$ depends on the flavor indices $a^{\prime}$ and $a$, a momentum fraction variable $z$ and the scale $\mu^{2}$. We also allow it to depend on the momenta and flavors $\{p,f\}_{m}$ of the partonic statistical state. Then $K$ has a perturbative expansion beginning with

[TABLE]

In order to keep the presentation in this paper as simple as possible, we set the masses of all five quarks d,u,s,c,b to zero. In a more complete picture, one uses a variable flavor number scheme. Then, even with $\overline{\mathrm{MS}}$ evolution, we change the renormalization scheme when the scale $\mu$ decreases past the bottom quark mass and then the charm quark mass. Thus the $\beta$ function for $\alpha_{\mathrm{s}}$ evolution and the parton evolution kernel change. For the shower oriented parton distributions $f_{a/A}(\eta_{\mathrm{a}},\mu^{2})$ and $f_{b/B}(\eta_{\mathrm{b}},\mu^{2})$, there is a different dependence on the c- and b-quark masses in the evolution kernels compared to what one uses in the variable flavor number version of $\overline{\mathrm{MS}}$ evolution.

There are various possible versions of $K^{(\mathrm{a})}_{aa^{\prime}}(z,\mu^{2},\{p,f\}_{m})$. The simplest is

[TABLE]

The first order $\overline{\mathrm{MS}}$ DGLAP kernel is

[TABLE]

Here $C_{a}=C_{\mathrm{F}}$ and $\gamma_{a}=3C_{\mathrm{F}}/2$ when $a$ is a quark flavor and $C_{\mathrm{g}}=C_{\mathrm{A}}$, $\gamma_{\mathrm{g}}=11C_{\mathrm{A}}/6-2T_{\mathrm{R}}n_{\mathrm{f}}/3$. The functions $P_{aa^{\prime}}^{\rm reg}(z)$ and $P_{aa^{\prime}}^{(\epsilon)}(z)$ are

[TABLE]

The kernel $K^{(\mathrm{a},1)}_{aa^{\prime}}(z,\mu^{2},\{p,f\}_{m})$ is a distribution in $z$ with a singularity at $z\to 1$. The singularity is represented by the + prescription in the first term of Eq. (179) and by the term proportional to $\delta(1-z)$. The coefficient of $\delta(1-z)$ is associated with how the virtual loop function $D^{(0,1)}(\mu^{2})$ is treated. We have here taken a simple choice based on what is in Ref. NSThreshold , but other choices are possible.

The logarithms of $Q_{\textsc{H}}^{2}/[(1-z)2p_{\mathrm{a}}\cdot Q_{\textsc{H}}]$ in Eq. (179) come about as follows. We attempt to calculate \big{(}{1}\big{|}\mathcal{X}(\mu^{2}) using the definition (79) of $\mathcal{X}(\mu^{2})$. We look at emissions of a parton in the initial state. Call the virtuality associated with this splitting $|k^{2}|=2\hat{p}_{\mathrm{a}}\cdot\hat{p}_{m+1}$. We integrate over $|k^{2}|$ and over the momentum fraction $z$. We use the hardness variable $\Lambda^{2}$, Eq. (55), used in Deductor. This means that there is an upper bound for the integration over $|k^{2}|$,

[TABLE]

There is an infrared divergence coming from the $|k^{2}|\to 0$ limit of the integration. This divergence is regularized by integrating in $4-2\epsilon$ dimensions. Now, dimensional regularization effectively acts as an infrared cutoff on the transverse momentum

[TABLE]

Thus we integrate over $|k_{\mathrm{T}}^{2}|$ with an upper bound

[TABLE]

The $1/\epsilon$ pole produced by the integration over $|k_{\mathrm{T}}^{2}|$ is removed by the factor $\mathcal{Z}_{\mathrm{F}}(\mu^{2})$ in $\mathcal{X}(\mu^{2})$. This leaves us with a $\log(Q_{\textsc{H}}^{2}/[(1-z)2p_{\mathrm{a}}\cdot Q_{\textsc{H}}])$, which multiplies the DGLAP splitting kernel $P_{a\hat{a}}(z)$. This remaining contribution does not have a $1/\epsilon$ pole. It has a term $\log(\mu^{2}/\mu_{\mathrm{s}}^{2})$, where $\mu_{\mathrm{s}}^{2}$ is the scale that defines the upper cutoff in the momentum integration and $\mu^{2}$ is the renormalization scale. We do not see this logarithm because we set these scales equal to each other. The Feynman rules for the splitting functions have some explicit $\epsilon$ dependence, giving a function of the form $f(\epsilon)/\epsilon=f(0)/\epsilon+f^{\prime}(0)+\mathcal{O}(\epsilon)$. The term $f^{\prime}(0)$ gives us the contributions $P_{aa^{\prime}}^{(\epsilon)}(z)$.

This calculation leaves us with an order $\alpha_{\mathrm{s}}$ contribution to \big{(}{1}\big{|}\mathcal{X}(\mu^{2})=\big{(}{1}\big{|}\mathcal{V}(\mu^{2}) that we need to eliminate because it does not vanish in the limit $\mu^{2}\to 0$. If we did not eliminate this term, we would lose Eq. (133). The offending contribution can be removed by the factor $\mathcal{K}(\mu^{2})$ in $\mathcal{X}(\mu^{2})$ if we choose the definition (179).

We also note that if we use $k_{\mathrm{T}}^{2}$ instead of $\Lambda^{2}$, Eq. (55), as the shower hardness variable, then there is no $\log(1-z)$. Then by redefining $\mu_{\mathrm{s}}^{2}$ by a factor $Q_{\textsc{H}}^{2}/[2p_{\mathrm{a}}\cdot Q_{\textsc{H}}]$, we obtain

[TABLE]

That is, the shower oriented parton distribution functions are close to the $\overline{\mathrm{MS}}$ parton distribution functions. However, they are not quite equal. That is because for the shower oriented parton distribution functions we are imposing an ultraviolet cutoff with a theta function, while with the $\overline{\mathrm{MS}}$ prescription we subtract a pole $1/\epsilon$.

We now turn to a construction that puts a different cut on parton splitting, leading to a less simple kernel $K^{(\mathrm{a})}_{a,a^{\prime}}(z,\mu^{2},\{p,f\}_{m})$. We retain the ultraviolet cut (182) when $\mu_{\mathrm{s}}^{2}$ is not too small. However, when $\mu_{\mathrm{s}}^{2}$ is small, the upper bound Eq. (184) can be very small indeed. We can relax this cut to

[TABLE]

where $m_{\perp}^{2}$ is, say, $1\ \mathrm{GeV}^{2}$. Then the matching kernel that defines the shower adapted parton distributions is

[TABLE]

When $m_{\perp}^{2}\ll\mu_{\mathrm{s}}^{2}$, this reduces to our previous definition. When $\mu_{\mathrm{s}}^{2}<m_{\perp}^{2}\,Q_{\textsc{H}}^{2}/2p_{\mathrm{a}}\cdot Q_{\textsc{H}}$, this becomes, after using Eq. (180),

[TABLE]

If we write Eq. (177) to first order as

[TABLE]

and use Eq. (188), recognizing that the $\overline{\mathrm{MS}}$ kernel generates scale changes in f_{a/A}^{{\overline{\phantom{\raisebox{-0.6pt}{\scriptstyle{\textsc{ms}}}}}\hskip-9.5pt\scriptstyle{\textsc{ms}}}}(\eta_{\mathrm{a}},\mu^{2}), we find

[TABLE]

That is, with this definition, for small values of the scale $\mu^{2}$, the shower oriented parton distribution functions approximately equal the $\overline{\mathrm{MS}}$ parton distribution functions at scale $m_{\perp}^{2}$. However, at larger scales, the shower oriented parton distribution functions evolve differently from the $\overline{\mathrm{MS}}$ ones. Deductor uses a definition similar to this, except with non-zero c and b quark masses and a corresponding variable flavor number scheme.141414The contribution from $P^{(\epsilon)}_{a\hat{a}}(z)$ is ignored in Deductor.

Appendix C Choosing $Q_{\textsc{H}}$ and $\mu^{2}$ dynamically

In Sec. IV.3, we introduced a vector $Q_{\textsc{H}}$ that is used to set scales and to help define one measure of hardness, $\Lambda^{2}$, that can be used in the shower. (See Eq. (55).) We stated that one can use the intended measurement operator to set $Q_{\textsc{H}}$ globally. There is another possibility: we can use the statistical state \big{|}{\rho_{\textsc{h}}}\big{)} that represents the hard scattering, Eq. (58). This state has the expansion

[TABLE]

Our example process in this paper is Higgs boson production. With this example, the Higgs boson momentum $p_{\textsc{h}}$ is part of $\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}$. We can set $Q_{\textsc{H}}^{2}=m_{\textsc{h}}^{2}$ and set the rapidity of $Q_{\textsc{H}}$ to be the rapidity of $p_{\textsc{h}}$ while letting the transverse part of $Q_{\textsc{H}}$ be zero. Another example would be jet production, for which the Born process is the production of two jets. In this case, we could use whatever infrared safe jet algorithm we like to find the two highest $P_{\mathrm{T}}$ jets in $\{p\}_{m}$. Letting the momenta of these jets be $P_{1}$ and $P_{2}$, we could set $Q_{\textsc{H}}^{2}=(\vec{P}_{1,\mathrm{T}}^{2}+\vec{P}_{2,\mathrm{T}}^{2})/2$, set the rapidity of $Q_{\textsc{H}}$ to the rapidity of $P_{1}+P_{2}$, and set the transverse part of $Q_{\textsc{H}}$ to zero.

This procedure requires some additional definitions since we have defined \big{|}{\rho_{\textsc{h}}}\big{)}=\mathcal{D}^{-1}(\mu_{\textsc{h}}^{2})\big{|}{\rho(\mu_{\textsc{h}}^{2})}\big{)}, so that we need $Q_{\textsc{H}}$ and $\mu_{\textsc{h}}^{2}=Q_{\textsc{H}}^{2}$ to create \big{|}{\rho_{\textsc{h}}}\big{)}, but we cannot use these variables before they have been defined.

Let us see what is needed in the case of jet production, for which the number of final state partons at the Born level is $m=2$. We content ourselves with what happens with an NLO hard cross section. We have

[TABLE]

We need to implement setting $\mu^{2}$ to $\mu_{\textsc{h}}^{2}$. In the case of the argument of $\alpha_{\mathrm{s}}$, this is a simple replacement. For the Born statistical state \big{|}{\rho^{(2,0)}}\big{)}, there is nothing further to do. In the terms \big{|}{\rho^{(2,1)}(\mu^{2})}\big{)} and \mathcal{D}^{(0,1)}(\mu^{2})\big{|}{\rho^{(2,0)}(\mu^{2})}\big{)}, there are cancelling poles. Then from $\mathcal{D}^{(0,1)}(\mu^{2})$ there is a left over $\log{\mu^{2}}$ from the cutoff $\mu_{\mathrm{s}}^{2}$ in the integration over the virtual loop graph. We simply have to set $\mu^{2}=\mu_{\textsc{h}}^{2}$ here. In the term \big{|}{\rho^{(3,0)}}\big{)}, there is nothing further to do. In the term \mathcal{D}^{(1,0)}(\mu^{2})\big{|}{\rho^{(2,0)}}\big{)}, the momenta $\{p\}_{3}$ determine the momenta $\{p\}_{2}$ in \big{|}{\rho^{(2,0)}}\big{)} and also the splitting variables $k^{2},z,\phi$ in $\mathcal{D}^{(1,0)}(\mu^{2})$. Increasing the ultraviolet cutoff $\mu_{\mathrm{s}}^{2}$ provides an increased range for $k^{2}$. As a result, the range in $\{p\}_{3}$ that is covered increases. Setting $\mu^{2}=\mu_{\textsc{h}}^{2}$, we check whether the point $\{p\}_{m}$ is inside the allowed range. If it is, we multiply by 1. If it is not, we multiply by zero.

It seems clear that this procedure works beyond NLO, although it becomes more complicated. We work term by term in the expansion of \big{|}{\rho_{\textsc{h}}}\big{)}. Once we have set $Q_{\textsc{H}}$ and $\mu_{\textsc{h}}^{2}$, we set $\mu^{2}=\mu_{\textsc{h}}^{2}$ in the argument of $\alpha_{\mathrm{s}}$ and in all explicit logs that come from virtual loop integrals. We check whether the point $\{p\}_{m}$ is generated by a subtraction term with an ultraviolet cutoff $\mu_{\textsc{h}}^{2}$. If not, we omit this term. Our description here has been algorithmic. To formulate this in terms of operators on the statistical space requires additional notation, which we omit.

Appendix D Renormalization

In this paper we use $\overline{\mathrm{MS}}$ renormalization. In particular, this defines $\alpha_{\mathrm{s}}$ and the $\overline{\mathrm{MS}}$ parton distributions that we start with. In the parton shower algorithm that we obtain, each element of the calculation is infrared finite and ultraviolet finite in four dimensions. However, part of the cancellation of infrared divergences is tied to the removal of ultraviolet divergences by renormalization. For this reason, the details of the ultraviolet renormalization scheme are significant. In this appendix, we gather the most important formulas that we use, mostly following Ref. JCCbook .

D.1 Renormalization of the QCD coupling

In the $\overline{\mathrm{MS}}$ scheme, the renormalization of the coupling is151515We find this definition useful for our purposes. Ref. JCCbook uses a different strong coupling, $\tilde{\alpha}_{\mathrm{s}}$, with $\alpha_{\mathrm{s}}=S_{\epsilon}\tilde{\alpha}_{\mathrm{s}}$.

[TABLE]

where $\mu$ is the renormalization scale, $\alpha_{\mathrm{s}}^{\mathrm{bare}}$ has mass dimension $2\epsilon$, and

[TABLE]

The renormalization constant of the strong coupling is given as a sum,

[TABLE]

The scale independent coefficients of the singularities, $Z_{\alpha}^{[n,k]}$, can be given in terms of the expansion parameters of the $\beta(\alpha_{\mathrm{s}})$ function by a recursion relation:

[TABLE]

The running coupling obeys the evolution equation

[TABLE]

where the first two $\beta_{i}$ coefficients are

[TABLE]

D.2 Renormalization of the parton distribution functions

The $\overline{\mathrm{MS}}$ parton distribution functions enter the cross section formula (42) in the form

[TABLE]

The renormalization factor $\mathcal{Z}_{F}$ is a product, so that each of the two parton distributions is transformed separately. The renormalization kernel $Z_{F}$ relates the renormalized parton distribution to the bare parton distribution:

[TABLE]

The kernel has a perturbative expansion

[TABLE]

It follows from the requirement that $f^{\rm bare}_{a/A}(\eta)$ is independent of $\mu^{2}$ and f^{\overline{\phantom{\raisebox{-0.6pt}{\scriptstyle{\textsc{ms}}}}}\hskip-9.5pt\scriptstyle{\textsc{ms}}}_{a/A}(\eta,\mu^{2}) has no poles that the renormalized parton distribution function obeys the DGLAP evolution equation,

[TABLE]

where

[TABLE]

with

[TABLE]

The coefficients of $1/\epsilon$ to higher powers are then determined by the recursion relation

[TABLE]

We can also write the evolution kernel as

[TABLE]

The factors contain poles $1/\epsilon^{k}$, but the poles cancel.

D.3 Renormalization scale dependence

The physical states \big{|}{\rho(\mu^{2})}\big{)} defined in Eq. (43) represent the quantum density operator of the partonic scattering. It is constructed from amplitudes \big{|}{M(\{p,f\}_{m})}\big{\rangle} and conjugate amplitudes \big{\langle}{M(\{p,f\}_{m})}\big{|}. We include the proper Lehmann-Symanzik-Zimmermann (LSZ) factors for the incoming and outgoing partons so that the amplitudes are S-matrix elements in the renormalized theory. With massless partons and dimensional regularization in Feynman gauge, this amounts to multiplying by field strength renormalization factors $Z_{\psi}^{-1/2}$ or $Z_{A}^{-1/2}$ for each external leg of an amputated graph. Then \big{|}{\rho(\mu^{2})}\big{)} is independent of the renormalization scale if calculated at all perturbative orders. When we calculate it up to order $\alpha_{\mathrm{s}}^{k}$ as in Eq. (43), we have

[TABLE]

Let us consider next the infrared sensitive operator $\mathcal{D}(\mu^{2})$. This operator depends on two scales, the renormalization scale $\mu^{2}$ and shower scale $\mu_{\textsc{s}}^{2}$. We have set these scales to equal each other, but in this appendix we highlight their separate roles. Thus we write $\mathcal{D}(\mu^{2},\mu_{\textsc{s}}^{2})$ with two arguments and let the operator with only one argument denote

[TABLE]

For any basis state \big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)}, the state \mathcal{D}(\mu^{2},\mu_{\textsc{s}}^{2})\big{|}{\{p,f,s,s^{\prime},c,c^{\prime}\}_{m}}\big{)} is to be defined so as to have the properties of a physical statistical state \big{|}{\rho(\mu^{2})}\big{)}. It is constructed in the renormalized theory and has the proper field strength renormalization factors for its external partons so that $\mathcal{D}(\mu^{2},\mu_{\textsc{s}}^{2})$ is independent of the renormalization scale up to the order that we calculate. That is, if $\mathcal{D}(\mu^{2},\mu_{\textsc{s}}^{2})$ is calculated up to order $\alpha_{\mathrm{s}}^{k}$, we have

[TABLE]

In order to calculate the parton shower generators $\mathcal{S}_{\rm pert}(\mu^{2})$ and $\mathcal{S}(\mu^{2})$, we need the total derivative of $\mathcal{D}(\mu^{2})$ with respect to the common scale. Since $\mathcal{D}(\mu^{2},\mu_{\textsc{s}}^{2})$ is independent of $\mu^{2}$, this is actually the derivative with respect to $\mu_{\textsc{s}}^{2}$:

[TABLE]

The perturbative expansion of this is

[TABLE]

Using this, the generators of the perturbative shower operator can be obtained from Eqs. (73) and (49) as

[TABLE]

This simplifies the shower generators. The operators $\mathcal{D}^{(n)}(\mu^{2},\mu_{\textsc{s}}^{2})$ always contains a theta function like $\theta(\Lambda^{2}<\mu_{\textsc{s}}^{2})$ or the equivalent constraint for the loop contributions. The partial derivative turns one of these theta functions into a Dirac delta function.

The inclusive infrared finite operator $\mathcal{V}(\mu^{2})$ is defined by the condition (92). This operator is derived from the infrared sensitive operator $\mathcal{D}$ and, just like $\mathcal{D}$, depends on both the renormalization and shower scales. Thus we write $\mathcal{V}(\mu^{2},\mu_{\textsc{s}}^{2})$ with two arguments and let the operator with only one argument denote

[TABLE]

Using Eqs. (79), (92), we have

[TABLE]

Recall that $\left[\mathcal{F}(\mu^{2})\circ\mathcal{K}(\mu^{2})\circ\mathcal{Z}_{F}(\mu^{2})\right]$ is independent of $\mu^{2}$. We have just seen that $\mathcal{D}(\mu^{2}\!,\mu_{\textsc{s}}^{2})$ is independent of $\mu^{2}$ at fixed $\mu_{\textsc{s}}^{2}$. Recall from Sec. VI.4 that we fix the color and spin content of $\mathcal{V}(\mu^{2},\mu_{\textsc{s}}^{2})$ to make Eq. (214) work for \big{(}{1}\big{|}\mathcal{V}(\mu^{2},\mu_{\textsc{s}}^{2}). Once we have done that for one choice of $\mu^{2}$ at a given $\mu_{\textsc{s}}^{2}$, we can keep the same operator $\mathcal{V}(\mu^{2},\mu_{\textsc{s}}^{2})$ for other choices of $\mu^{2}$ and Eq. (214) will continue to hold. That is, we can define $\mathcal{V}(\mu^{2},\mu_{\textsc{s}}^{2})\mathcal{F}(\mu^{2})$ so that it is independent of $\mu^{2}$ up to whatever order we calculate:

[TABLE]

This gives us

[TABLE]

The generator $\mathcal{S}_{\mathcal{V}}(\mu^{2})$ of $\mathcal{U}_{\mathcal{V}}(\mu_{2}^{2},\mu_{1}^{2})$ is defined in Eq. (114) as

[TABLE]

This gives us

[TABLE]

Here the first term represents the evolution of the perturbative part of $\mathcal{V}$ and the second term gives the evolution of the parton distribution functions. This generalizes Eq. (124) for the first order contribution to $\mathcal{S}_{\mathcal{V}}$.

Finally the generator of the probability preserving shower evolution operator is given in Eqs. (115) and (73). We can simplify this by using Eqs. (210) and (218)

[TABLE]