Rapidity distribution at soft-virtual and beyond for n-colorless particles to N$^4$LO in QCD

We present a systematic framework to study the threshold contributions of the differential rapidity distribution for the production of any number of colorless particles in the hadronic colliders. This has been achieved based on the universality structure of the soft enhancements associated with the real emissions, along with the factorization property of the differential cross section and the renormalization group invariance. In this formalism, we present a universal soft-collinear operator to compute the soft virtual differential cross section for a generic $2\rightarrow n$ scattering process up to next-to-next-to-next-to-next-to-leading order (N$^4$LO) in perturbative QCD. We also provide a universal operator to perform the threshold resummation to next-to-next-to-next-to-leading logarithmic (N$^3$LL) accuracy. We explicitly present the approximate analytical results of the rapidity distributions at N$^4$LO and N$^3$LL for the Higgs boson production through gluon fusion and bottom quark annihilation, and also for the Drell-Yan production at the hadronic collider. \textcolor{black}{We extend our framework to include the next to threshold contributions for the diagonal partonic channels.


Introduction
The upcoming era of high energy physics will confront a huge boom of data driven by the upgraded run of the Large Hadron Collider (LHC). The High-Luminosity LHC will also come into effect in a few years of time. In order to fully exploit the increased quantity of data which will not only help to detect a rare phenomena but also improve the precision, an enhanced amount of effort by the theoretical physia e-mail: taushif@mpp.mpg.de b e-mail: ajjathah@imsc.res.in c e-mail:poojamukherjee@imsc.res.in d e-mail:ravindra@imsc.res.in e e-mail: aparnas@imsc.res.in cists will be called for. In improving theoretical precision, the higher order quantum chromodynamics (QCD) and electroweak (EW) corrections play an important role. Among different observables, since the differential cross-section allows a wider range of comparisons with the experimental data, over the past few decades several attempts have been made to incorporate the higher order QCD and EW radiative corrections to this observable. The topic of this article is concerning the differential cross-section with respect to rapidity, in particular, we address the question of computing the higher order QCD corrections to this observable for any generic process at a hadron collider with all the final state particles as colorless.
Despite its high importance, unlike the inclusive crosssection, the differential rapidity distribution and its radiative corrections are computed only for a limited number of scattering processes. The rapidity distributions in Drell-Yan and of the scalar Higgs boson were computed to next-to-next-toleading order (NNLO) QCD in refs. [1,2] and [3], respectively. In case of the scalar Higgs boson produced through gluon fusion, the next-to-NNLO (N 3 LO) QCD correction was incorporated in ref. [4]. Shortly before, it was approximated in ref. [5] in the formalism of q T -subtraction. For the Higgs boson production through bottom quark annihilation, it was computed to NNLO in ref. [6].
Needless to say, achieving a full QCD correction to any order is not easy and with increasing perturbative order, the complexity level increase substantially which often prevents us from achieving it. In absence of the full QCD correction, it is often desirable to find an alternative method averting the full complexity to capture the dominant contribution. In this article, we discuss such a method, called soft-gluon or soft-virtual (SV) approximation. In the soft limit, the momenta of all the real emission diagrams are assumed to be infinitesimally small which leads to an all order exponentiation of this contribution. In ref. [7,8], a formalism to in-corporate the soft-gluon contribution to the rapidity distribution for the production of a colorless final state in hadron collider is presented which, in this article, is extended to the case of any number of final state colorless particles. The formalism is based on QCD factorization, which dictates that the soft part of the real emission diagrams factorizes from the hard contribution, and renormalization group (RG) invariance. The factorized soft part is conjectured to fulfil a Sudakov type differential equation with respect to the final state invariant mass square and as a consequence, it is found to get exponentiated which not only provides us with the fixed order result under soft limit but also enables us to perform a resummation in the soft limit. For the production of arbitrary number of colorless particles in hadronic collision, the soft part essentially remains identical to the case of Sudakov type process since the real emission can only takes place from the initial state partons. Using this idea, in this article, we extend the formalism [7] to the case of 2 → n scattering, where n denotes the number of final state colorless particles. We show that combining the virtual matrix element that captures the process dependence, the universal soft part and mass-factorization kernel in an elegant way, we can calculate the SV differential rapidity distribution for any generic 2 → n process. In addition, we also show how naturally it leads to the threshold resummation for the same observable. Needless to say, in some kinematic regions such as when partonic center-of-mass energy becomes very close to threshold, the logarithms often become so large that it endangers the validity of the perturbation theory and we are left with no other choice but to perform an all order resummation in order to get a reliable prediction.
In the literature, several results for the rapidity resummation employing different methods are available. In ref. [9], following the conjecture given in [10], the resummation of rapidity of W ± gauge boson and in ref. [11] of Drell-Yan are computed in Mellin-Fourier (M-F) space. A detailed theoretical underpinnings and phenomenological implications of threshold resummation of rapidity are examined in ref. [12] emphasizing the role of prescriptions that take care of diverging series at a given logarithmic accuracy. Our method belongs to a category, so called direct QCD approach [13], which is based on ref. [7,8,14] resums the soft gluons in two dimensional Mellin space (M-M). In [15], the merits of different approaches are discussed in details.
Over the past decade, the formalism [7] for the 2 → 1 scattering has been tested extensively through the computations of the rapidity distribution of several Sudakov type process. In ref. [7], the observable was computed partially at N 3 LO under the SV approximation for the Higgs boson production through gluon fusion and for di-lepton in Drell-Yan. Later, some of us completed the full SV correction at N 3 LO in ref. [16] by incorporating the missing part from real emission diagrams. One of the salient features of our formalism is that the soft part that enters into the rapidity distribution is shown to be connected to the respective part of inclusive cross-section through a very simple relation involving gamma function of the dimensional regulator. This relation was used to extract the required soft part from the respective quantity of the SV cross-section [17,18,76]. The SV rapidity distribution at N 3 LO is also computed for the Higgs boson production through bottom quark annihilation [19]. The threshold resummed results at next-to-next-to-leading log (NNLL) for the Higgs boson [14] and Drell-Yan [20] are also obtained. In this article, we extend the formalism for computing the SV differential rapidity distribution and its resummation to any 2 → n scattering process at the hadron collider and we provide the results to next-to-N 3 LO (N 4 LO) and next-to-NNLL (N 3 LL) in terms of generic quantities in Online Resource. We also provide the approximate results of the rapidity distributions at N 4 LO and N 3 LL explicitly for the Higgs boson production through gluon fusion and bottom quark annihilation, and for the Drell-Yan production. The approximation comes from the unavailability of the full four loop virtual matrix elements and the soft-collinear distributions. In recent times, there has been a surge of interest to understand the next-to SV (NSV) contributions to inclusive as well as differential distributions, see [21][22][23][24][25][26][27][28][29][30][31][32][33] for more details. In [34,35,77], some of us have demonstrated how these contributions can be systematically included for single colorless particle productions. This was done for the diagonal partonic channels. In the context of rapidity distributions, it was shown that like soft virtual terms, NSV terms can also be resummed to all orders both in z and double Mellin {N 1 , N 2 } spaces [35]. The goal of this article is to present the formal methodology of computing threshold rapidity corrections including NSV terms for any generic process of 2 → n kind at hadron collider.
The paper is organized as follows: in section 2, we introduce the notion of soft-virtual correction in the context of differential rapidity distribution and then describe our formalism in details in section 2.1. The universality of soft part leads us to define a quantity called the differential softcollinear operator that essentially captures the process independent part is also introduced in section 2.1. In the next section 3, we extend our formalism to incorporate the threshold resummation of the rapidity distribution.In section 4, we discuss how next-to soft virtual corrections can be included in the rapidity distribution in particular for diagonal channels. All the results to N 4 LO and N 3 LL in terms of generic quantities are presented in Online Resource in Mathematica format which can be used for any generic 2 → n process once the corresponding virtual matrix element becomes available provided the soft distribution function for Sudakov type process is known.

Soft-virtual rapidity distribution
We begin by introducing the regime of soft gluon contribution to differential rapidity distribution for the production of n-number of colorless particles in hadron collisions within the framework of perturbative QCD. Our prescription that we will develop subsequently to capture this contribution is within the scope of QCD improved parton model where the collinear divergences factorize to all orders in strong coupling constant α s = g 2 s /4π. We consider a generic hadronic collision between two hadrons H 1,(2) having momentum P 1(2) that produces a final state consisting of nnumber of colorless particles, denoted as F i (q i ) Through the quantity X, we represent an inclusive hadronic state. q i stands for the momentum of corresponding colorless particle F i . We denote the invariant mass square of the final state by q 2 which is related to the momenta {q i } through q 2 = (∑ i q i ) 2 . Without loss of generality, the rapidity of the final state invariant mass system is defined as The differential rapidity distribution at the hadronic level can be written as where σ B is the leading order inclusive cross section. The dimensionless variables, τ and z, are respectively defined as the ratios of invariant mass square of the final state to the square of hadronic (S) and partonic (ŝ) center-of-mass energies i.e.
We denote the fraction of the initial state hadronic momentum carried by the partons (a, b) that take part in the scattering at the partonic level as x 1 (2) , and these are constrained through the relation τ = zx 1 x 2 as reflected by the presence of the respective δ -function in the definition of W . The remaining delta functions reflects the momentum conservation and also the rapidity of final state invariant mass system as defined by (2).f a(b) are the bare parton distribution functions (PDF). The spin and color averaged square of the scattering matrix element is denoted by |M ab | 2 . The corresponding m + n-particle phase space is given by [dPS m+n ] where the integer n indicates the number of colorless particles at the final state. Note that the numerical value of the integer m depends on the number of radiated partons which is solely controlled by the perturbative order we are interested in. The scattering matrix element or equivalently the partonic differential distribution can be perturbatively expanded in powers of strong coupling constant a s (µ 2 R ) ≡ α s (µ 2 R )/4π with µ R being the renormalization scale. In this article, we confine ourselves to the regime where the leading order (LO) processes can only be initiated through color neutral quark or gluon channels, with the corresponding momenta p 1 (2) . Moreover, we are interested in computing the differential rapidity distribution only in the soft limit which constrained all the partonic radiation to be only soft i.e. those can only have nearly vanishing momenta. The goal of this article is to present a prescription to compute the leading contribution of the differential rapidity distribution in the soft limit which is commonly referred as soft-virtual (SV) contribution. In addition, we also discuss how the current formalism can be extended to study the next-to-soft virtual (NSV) contribution. In order to define the soft limit for the rapidity distribution, we choose to work with a set of symmetric scaling variables x 0 1(2) instead of y and τ which are related through Note that unlike the inclusive cross-section, the choice of variables which one needs to take in order to define the soft limit is not unique and as it turns out, our choice of these new set of variables is crucial for our prescription. In terms of these variables, the partonic contributions arising from the subprocesses are found to depend on the ratios which play the role of scaling variables at the partonic level. After performing the mass-factorization and also evaluating the δ -function integration over z, the W (τ, q 2 , y) in (3) can be rewritten as where µ F is the factorization scale. In (8) , ∆ d,ab is the partonic coefficient function which can be expressed as where δ 1 and δ 2 are the shorthand notations for the two δ -functions appearing in the last line of (8). In the above equation, M ab is the UV renormalized pure virtual corrections for the 2 → n scattering while the quantity inside the square brackets encapsulates all the real emission contributions as it is already normalised by the pure virtual contributions after the phase space integration. The Γ 's are the mass-factorization kernels which are introduced to cancel the initial state collinear singularity present in both virtual and real emission contributions in (9). In the next section, we will study the IR factorisation property of M ab . Being a scattering process containing n-number of final state colorless particles, the partonic coefficient function does, in fact, depend on the Mandelstam variables constructed out of all the independent external momenta which is concisely denoted through {p j · q k }. In order to find the definition of soft limit in terms of the new partonic scaling variables, we take the double Mellin moment of W with respect to the variables N 1(2) which turns out to be All the quantities with functional dependence of N 1(2) are in Mellin space where the soft limit is defined by the simultaneous limit of N 1(2) → ∞. In terms of partonic scaling variables this condition gets translated to z 1(2) → 1. Note that we normalize the coefficient function ∆ d,ab in such a way that at the leading order W satisfies In the following section, we will present the prescription to calculate the infrared safe SV differential rapidity distribution to N 4 LO in QCD for any 2 → n scattering process which can be computed order by order in perturbation theory: λ is the order of strong coupling constant at the leading order partonic process. Here the arguments in ∆ sv d,ab are suppressed.

Universal soft-collinear operator for SV rapidity distribution
In this section, we setup a framework to compute the softvirtual corrections to the rapidity distribution to all orders in strong coupling constant. The infrared safe SV rapidity distribution can be obtained by combining the ultraviolet (UV) renormalized virtual matrix elements with the soft gluon contribution and performing appropriate mass factorization to get rid of initial state collinear singularities.It is wellknown that the combined soft and collinear divergences, conveniently denoted as infrared (IR), in virtual matrix elements factorize from the corresponding UV renormalized part to all orders in perturbation theory and thereby in dimensional regularization we can write with the space-time dimensions D = 4 + ε. Without loss of generality, we choose the renormalization scale to be equal to the scale of aforementioned factorization which, of course, in general can be different. Upon multiplying the renormalization constant Z ab,IR , the IR divergent part of the UV renormalized matrix element M ab gets compensated and we end up with the finite part of the matrix element M ab,fin . The renormalization constant is a universal quantity as it is independent of the details of the process, it only depends on the nature of external color particles. It is fully independent of the number and nature of external colorless particles. The exponentiation of the Sudakov form factor in dimensional regularization was first demonstrated in refs. [36][37][38]. For the multiparton amplitudes, the universal structure of the IR divergences in dimensional regularization, up to the single pole in ε at two-loop, was first depicted in ref. [39] and subsequently analyzed in ref. [40]. For scattering involving only two external colored particles, the universality of the single pole at two-loop was established in terms of soft and collinear anomalous dimensions in ref. [41] which was verified at three-loop in ref. [42]. An all order conjecture on the form of IR divergences in QCD for any generic scattering process is depicted in terms of soft anomalous dimension matrix in ref. [43,44]. An elaborated discussion on the universality of subleading infrared poles in form factors is given in ref. [45]. We provide the results of the renormalization constant to four loops in Online Resource. Before moving further, we wish to iterate a well-known fact. For processes involving only conserved operator, such as Drell-Yan, the coupling constant renormalization is sufficient to get rid of all the UV divergences. However, for other processes, such as the Higgs boson production in heavy quark effective theory, an additional renormalization which often called the operator renormalization is needed. This is a property inherent to the operator itself.
In order to get the infrared safe and finite differential rapidity distribution, we need to combine the UV renormalized virtual matrix element to the real emission contributions in the soft limit and perform mass factorization which ensures the removal of collinear singularities arising from the initial state colored particles. Therefore, the universal nature of IR divergences in virtual matrix element implies that the combined contribution from the real emission diagrams and mass-factorization kernels must exhibit the same universality. By employing the criteria of universal IR structure and imposing the finiteness property of the rapidity distribution, we develop the prescription to compute the rapidity distribution under SV approximation for any generic 2 → n scattering process and present the result in terms of universal quantities to N 4 LO QCD. Once the pure virtual matrix element for any process of the kind under consideration becomes available, our expression can immediately be employed to calculate the SV rapidity distribution at that order in QCD.
In this article, we extend the prescription, which was introduced for the Sudakov type process in ref. [7], to the case of 2 → n scattering. We propose that the coefficient function for the rapidity distribution in (12) can be written as a Mellin convolution of the pure virtual contribution F , softcollinear distribution Φ d and mass-factorization kernel Γ , which read as where δ (z l ) = δ (1 − z l ) for l = 1, 2. The pure virtual contribution is captured through the form factor F aā that is defined as where M (k) aa represents the k-th order UV renormalized matrix element of the underlying partonic level process a(p 1 )+ a(p 2 ) → ∑ n i=1 F i (q i ). The symbol "C " stands for the convolution whose actions on a distribution g(z 1 , z 2 ) is defined as where ⊗ denotes Mellin convolution. In this article in the context of SV corrections, we encounter only δ (1 − z i ) and The contribution from the real emission diagrams is contained in Φ d,aa which we call as soft-collinear distribution. The soft divergences arising from the real emission and virtual diagrams, which are respectively encapsulated in Φ d and F , get cancelled. The final state collinear singularity is guaranteed to go away, as dictated by Kinoshita-Lee-Nauenberg (KLN) theorem, once the sum over all the final states is performed. The mass factorization kernel takes care of the initial state collinear singularities. As a result, the coefficient function ∆ sv d,aā in (14) becomes finite. By demanding the finiteness of this quantity we can put a constraint on the soft-collinear distribution which turns out to be a Sudakov type renormalization group (RG) equation. This has profound implications which not only reveals a significant amount of insights about the IR world but also it enables us to perform threshold resummation as we will see in the next section. To be more precise, the solution of the RG equation results an all order exponentiation of the soft-collinear distribution. So, the whole job of computing the SV correction depends on our ability to determine and explore the unknown distribution Φ d,aa to which we now turn to 1 .
As we have discussed, the soft-collinear distribution essentially captures the contribution arising from real emission diagrams which only can occur from colored partons. Naturally, Φ d,aa for Sudakov form factor i.e. 2 → 1 and 2 → n scattering should essentially be identical. The presence of more Mandelstam variables in the latter process just makes it more involved in its kinematic dependence when it is expressed in terms of {p j · q k }. However, in terms of the total invariant mass square of the final state colorless particles i.e. q 2 , it has to be exactly same as that of Sudakov process. In ref. [7], it was conjectured to satisfy a integrodifferential RG equation to all orders in QCD coupling constant. The underlying reason behind this all order conjecture is inspired by the akin integro-differential Sudakov equation [36,37,47] fulfilled by the form factor whose solution is present explicitly to five loops order in massless QCD in refs. [42,[48][49][50]. By integrating the differential rapidity distribution, we get the inclusive cross-section. Upon taking the Mellin moment with respect to the same Mellin variable N of this relation we get By taking the limit N → ∞ on both sides of this relation, we can relate the soft-collinear distributions in rapidity and that of inclusive cross-section. This is remarkable in a sense that given the soft-collinear distribution for inclusive crosssection, we can automatically calculate it for the rapidity.
Since this is the only quantity that is unknown in comparison to the ingredients for the computation of SV crosssection, we can immediately calculate the SV rapidity distribution. The Φ d,aa for the Sudakov form factor is deter-mined to NNLO in refs. [7] and in [16] at N 3 LO in QCD.
In the current article, for the first time we present the general analytical form of Φ d,aa in terms of universal quantities at N 4 LO for any generic 2 → n scattering. One of the most notable features of this quantity is it satisfies the maximally non-Abelian property i.e.
where the quadratic Casimirs in Adjoint and fundamental representations of SU(n c ) are denoted by C A = n c and C f = (n 2 c − 1)/2n c , respectively. The property says that the softcollinear distribution for quark and gluon initiated processes which are respectively denoted by qq and gg are related by simple scaling of quadratic Casimirs. This essentially signifies the universality of the real emission in the soft limit i.e. it is independent of the details of the process, it solely depends on the nature of the external partons. Needless to say, it is also quark flavour blind. The relation in (19) was explicitly verified to NNLO in refs. [7] and at N 3 LO in [16]. The flavour dependence of the Φ d,aā was exploited in ref. [19] in order to calculate the SV rapidity distribution at N 3 LO for the Higgs boson production in bottom quark annihilation. We expect the Casimir scaling to hold true to all orders in perturbation theory since it originates entirely from the softcollinear part of the differential cross-section, and therefore it would indeed be interesting to see whether truly it holds beyond N 3 LO with generalised Casimir scaling [51].
We decompose all the quantities into its singular (sing) and finite (fin) parts as where Through the decomposition of the quantities into singular and finite parts in (20), we put together all the singular components of the rapidity distribution into I d,aā which must be unit distribution δ (z 1 )δ (z 2 ) in order to get a finite ∆ sv d,aā . In (21), the form factor and the leading order matrix element are the only process dependent quantity. The remaining part which comprises of the finite segments of the soft-collinear distribution and mass factorization kernel is a process independent universal quantity which we call as differential soft-collinear operator The expression of S d,aa being process independent can be used for any generic 2 → n scattering process. Since we are confining our discussion to only those scattering processes with initial state quark-antiquark pair of same flavours or a pair of gluon, we conveniently rewrite the (21) as where I = qq, gg. We can calculate the SV coefficient function for the rapidity distribution order by order in perturbation theory by expanding it in powers of a s according to (12). We provide the results of S d,I and ∆ sv d,I for any generic process, which can be expressed in terms of universal lightlike cusp (A I ), eikonal or soft ( f I ) and virtual or collinear (B I ) anomalous dimensions along with the process dependent form factors, in Appendix A and Appendix B up to N 2 LO QCD, respectively. The results beyond N 2 LO can be found in the Online Resource provided with this article. Though the scale dependence can be restored employing the renormalization group evolution, nonetheless for users' convenience we provide the results of ∆ sv d,I by keeping the explicit scale dependence in Online Resource. Once the virtual matrix element becomes available, one can directly use our generic results to calculate the SV rapidity distribution for any generic process of the kind under consideration. The anomalous dimensions are expanded in powers of a s (µ 2 R ) as where X = A, B, f . Thanks to recent calculations, the lightlike cusp anomalous dimensions are available to four loops [52][53][54][55][56] in QCD. The soft and collinear anomalous dimensions can be extracted [41,42] from the quark and gluon collinear anomalous dimensions [57,58] through the conjecture [41] γ I = 2B I + f I to three loops. At four loop, only partial results are available in refs. [56,[59][60][61]. We present the new results of ∆ sv d,I for 2 → 1 processes such as the Drell-Yan and the Higgs boson productions through gluon fusion as well as bottom quark annihilation at N 4 LO in Appendix D. This has been achieved by using the explicit results of the recently computed four loop cusp anomalous dimension [52][53][54][55][56] along with the form factors for Drell-Yan and the Higgs boson productions which are approximately available to fourth order in QCD [56,[62][63][64][65][66]. The previously missing coefficients of δ (z 1 )δ (z 2 ) in ∆ sv d,I at N 4 LO were due to the missing O(ε 0 ) results of form factors and the soft gluon contributions at four loop. The full explicit results of the quark and gluon form factors corresponding to Drell-Yan and the Higgs boson production in gluon fusion are available at three loop to O(ε 2 ) in ref. [67,68] which are also required. The corresponding partial four loop form form factors are computed in several articles over the past few years [56,[62][63][64][65][66]. For the case of Higgs production through bottom quark annihilation, the three loop [65,69] and partial four loop results are available in ref.
[65]. Moreover, the one-and two-loop results are needed to expand to O(ε 5 ) and O(ε 3 ), respectively. Similarly for the coefficient which contributes to the O(ε 0 ) part of four loop form factor, the one-, two-, and three-loop results are needed to order O(ε 6 ), O(ε 4 ) and O(ε 2 ), respectively. The four loop explicit results which we present in this article are still incomplete due to the unavailability of the full explicit results for form factors as well as soft contributions resulting from the real emission processes at four loop, so far these are the state-of-the-art available results in the literature.
In Online Resource, the explicit expressions of all the anomalous dimensions including the QCD β -functions to three loops can be found. In the following section, we describe how our formulation of SV rapidity distribution naturally leads us to the soft gluon resummation for any generic 2 → n scattering process.

Threshold resummation and its universal soft-collinear operator
Here in this section, we develop the resummation formalism for the differential distribution with respect to the rapidity variable y, for the production of n-colorless particles which also paves the way for a wider range of comparisons with the experiments. Earlier we have seen that there exists an operator, S d,I , referred as the differential soft-collinear operator, which embeds the universality of all the soft enhancements associated with the soft gluon emissions in the production of n-colorless particles in the hadronic collision. The universality lies in the fact that the operator, S d,I in (23), for the SV differential cross-section depends only on the initial state partons and is completely independent of the hard process under study. Besides being the process independent operator, interestingly it also exhibits an exponential behaviour. Recall that the threshold resummation [46] relies on the fact that the soft contribution exponentiates to all orders in perturbation theory, owing to the Sudakov differential equation and the renormalization group invariance. Following the same argument we proceed towards the resummation formalism for differential cross-section as well.
The relevance of resummation of differential cross-section arises from the fact that, in the limit z 1(2) → 1, the logarithms of type a n s log , give rise to large contributions which could potentially spoil the reliability of the perturbative series. Hence a systematic way of exponentiating these large logarithms and resumming them to all orders in perturbation theory becomes indispensable. In ref. [13] it was shown, in the context of differential distribution with respect to the Feynman variable x F , that the potential logarithms which give dominant contributions in certain kinematic regions can be resummed to all orders in perturbation theory in Mellin-Mellin (M-M) space approach. This approach was also extended to rapidity distributions in the earlier works by one of the authors of this paper (See [14,20] for details). Note that this approach is different from the Mellin-Fourier (M-F) approach [10] proposed by Laenen & Sterman which was earlier discussed in the introduction. In M-F formalism partonic cross-section is expressed in terms of scaling variable z and rapidity variable y and then the threshold limit is taken only for z → 1 which resums delta (δ (1 − z)) and distributions ([ log(1−z) 1−z ] + ) in z, but for rapidity variable y only delta (δ (y)) piece is resummed. In ref. [20], one of the authors of this paper has made a detailed numerical comparison of M-M approach against the M-F approach and found that both the approaches converges to a few percent correction to the fixed order prediction at NNLL level. In the following we further extend the M-M approach and derive the resummation formalism for the production of n-colorless particles in a partonic collision.
Within the framework of M-M approach, both the partonic scaling variables z 1 (2) are simultaneously taken to the threshold limit 1 and the corresponding delta, δ (1 − z i ), and , are resummed to all orders in perturbation theory. Due to the involvement of convolutions in the z 1(2) space, the resummation is performed in two dimensional Mellin space where the differential crosssection is expressed in terms of simple normal products along with the aforementioned exponential structure in form. In the following we derive the generic formalism in terms of the Mellin variables N 1 and N 2 corresponding to the z 1 and z 2 variables, respectively. Hence the threshold limit z 1(2) → 1 translate to N 1(2) → ∞ in Mellin space and the large log-arithms proportional to log N 1 (2) are resummed to all orders in perturbation theory.
To derive the all order behaviour of the SV differential cross-section, ∆ sv d,I (z 1 , z 2 ), in the two dimensional Mellin space withN i = N i e γ E , we begin with the Mellin moment of the same, which takes the following form: γ E is the Euler-Mascheroni constant. In the previous section in (14), ∆ sv d,I is decomposed into constituents corresponding to the virtual as well as the soft-collinear real emission contributions. Now in this section, we further decompose those contributions into a process dependent and a process independent quantities. We denote the process dependent coefficient C I d,0 in the context of 2 → n scattering process as, Here C I d,0 accounts for all the finite contributions coming from the virtual corrections and the coefficients proportional to δ (z 1 )δ (z 2 ) of the real emission contributions. Besides, it also contains the finite part of the mass factorized kernel Γ I,fin in terms of log(µ 2 F /µ 2 R ) which results from the coupling constant renormalization. The quantity S d,I res,δ which we name as the differential soft-collinear operator for threshold resummation, embeds the δ (1 − z i ) contributions from the soft distribution function Φ d,I and from Γ I,fin in the following way: The subscript δ indicates δ (z 1 )δ (z 2 ) coefficients of the aforementioned quantities. In a similar way, we denote the process independent contributions to ∆ sv d,I as Φ res d,I which comprises of the terms proportional to plus distributions from Φ d,I and Γ I,fin . Mathematically it can be written as, where the subscript D indicates the terms proportional to plus distribution which includes, D i (z 1 )δ (z 2 ) , D i (z 2 )δ (z 1 ) and D i (z 1 )D j (z 2 ). Now following the approach given in ref. [7], the above equation can be written in an integral form which is given as, here the subscript + indicates the standard plus distribution and the other constants are defined asz i = (1 − z i ) and z 12 = q 2 z 1 z 2 . The finite functions, , are related to the threshold exponent D I i of inclusive cross section owing to the relation given in (18) ( See [7,14] for more details). For completeness, we provide the coefficients A I i and D I d,i in the Online Resource. Consequently the SV differential cross-section decomposes into a process dependent and a process independent way and can be re-written in the following form: Substituting (32) in (27) and after doing the two dimensional Mellin transformation systematically, we obtaiñ with ω = β 0 a s (µ 2 R ) log(N 1N2 ). The first coefficient of QCD β -function is denoted by β 0 ≡ (11C A − 2n f )/3, n f is the number of active light quark flavours. Here, the decomposition in the exponent is done in such a way that the coefficient G I d,N contains N 1(2) dependent terms, and the remaining ones are embedded in G I d,0 . Besides this, G I d,N (q 2 , µ 2 F , ω) also vanishes in the limit ω → 1. Needless to say that both of these coefficients has a universal structure in terms of the anomalous dimensions A I and process independent coefficients D I d and thus are dependent only on the incoming partons. Further we combine the N 1(2) independent coefficients G I d,0 with C I d,0 from (33) and define, which can be expanded in terms of a s (µ 2 R ) as, From (34) it can be seen, that the coefficient g I d,0 contains finite contribution from virtual corrections, differential soft-collinear operator for threshold resummation and N independent terms coming from Mellin transformation of plus distribution. Consequently, (33) gets modified as, where the exponent G I d,N can be organized as a resummed perturbation series in Mellin space as, The explicit form in (47) when expanded till k-th order in powers of a s (µ 2 R ), gives the logarithmically enhanced contributions to the fixed order results∆ sv d,I (N 1 ,N 2 ) up to the same order. The successive terms in the above series given in (37) along with the corresponding terms in (35) define the resummed accuracy as LL, NLL, NNLL, N 3 LL and so on. In general for N k LL accuracy, terms up to g I d,k+1 must be included along with g I d,0 up to order a k s (µ 2 R ). The general expression for the coefficients g I,i d,0 and g I d,i up to N 3 LL are provided in the Online Resource.
The coefficients G I d,N remains unaltered even for 2 → n scattering process owing to its universality. However, the process dependent coefficient function g I d,0 changes for the production of n-colorless particles due to the inclusion of process specific form factor via (28) and (34). The results of these coefficients appear as a product of N 1 and N 2 in the Mellin space, and all those terms which are only function of N 1 or N 2 cancel internally.
We have also observed that the coefficients g I d,0 and G I d,N coincides with their inclusive counterparts g I 0 and G IN respectively in the limit N 1 → N 2 → N, provided the coefficients D I d in (31) is expressed in terms of D I of inclusive soft distribution function using the relation (18). Hence we infer that all the above observations which hold true for 2 → 1 scattering processes are further extended and verified for any generic system of n-colorless particles in the final state.

Beyond soft virtual rapidity distribution
Having obtained the results for the SV part of coefficient functions, namely ∆ sv d,I , both in z i and N i spaces, we would like to extend our approach to include NSV terms. These are logarithms of the form log k (1 − z i ), k ≥ 0 present in coefficient functions that are often comparable or even larger than SV contributions. Recently, some of us have developed a formalism to study the NSV structure in inclusive reactions, namely, the production of leptons pairs in Drell-Yan process, Higgs boson production in gluon fusion or in bottom quark annihilation [34] and deep inelastic scattering and semi-inclusive e + e − annihilation processes [77]. Later the formalism has been extended to rapidity distributions in [35]. We considered the observables such as invariant mass (total cross section) and/or rapidity distributions of a lepton pair (a Higgs boson) in Drell-Yan (gluon fusion or bottom quark annihilation). We restricted ourselves to only diagonal channels where either quark and anti quark annihilate (light quarks in Drell-Yan, bottom quarks for Higgs boson) or gluon fusion in Higgs boson production. We used the principles, namely mass factorization, RG invariance that we had used for studying SV logarithms. Remarkably, the IR structure of the NSV terms in the coefficient functions can be understood in terms of the corresponding contributions in real emission processes and mass factorization kernels which are the building blocks. In the case of production of a single colorless state, we found that a part of the NSV terms is controlled by certain universal anomalous dimensions similar to the case of SV terms and the remaining terms in NSV depend on the process under study beyond second order in a s [34,35]. In particular, the differential soft-collinear operator for the NSV terms turns out to be process dependent. To the end, we have shown that both SV and NSV logarithms can be systematically summed up to all orders both in z i and N i space and this allowed us to predict certain SV and NSV logarithms for all i > n at every order in a i s , using results known up to order a n s . In the following, we apply the same formalism to study the NSV logarithms present in the rapidity of a state of ncolorless particles. We restrict ourselves to the diagonal coefficient functions (dCFs). For the diagonal channels, it is straightforward to show that the NSV terms in dCFs arise only from diagonal partonic sub processes and diagonal part of mass-factorization kernels. The pure virtual part of the partonic sub-processes does not contain any NSV terms and hence can be factored out from them. Hence, in the mass factorization formula, we need to keep both SV as well as NSV terms in Φ d,aa and in Γ aa . The latter is process independent and is known to third order [53,54]. Recall that the challenging task was to determine Φ d,aa . It contains singular and finite terms and the singularities are from IR region .The mass factorization demands that the soft part of singularities must cancel with those from pure virtual part of the partonic channels and the initial state collinear singularities should cancel against those from mass factorization kernels. The cancellation of singularities and the knowledge of dCFs allowed us to parametrise the SV part of Φ d,aa in terms of δ (1 − z l ), plus distributions containing logarithms and we can use the same method to determine the NSV logarithms. We find that the logarithmic structure of finite part in Φ d,aa does not depend on the process under consideration, while their numerical coefficients will depend on the process. As it was done for the SV case, these numerical coefficients can be determined using the inclusive results by working in the Mellin N space. To end, we obtain with and which again reduces to identity in z 1 , z 2 space, namely δ (1− z 1 )δ (1−z 2 ). We provide the results of NSV part of the dCFs in (38) denoted by ∆ nsv d,I andS d,I for any generic process with I = (gg, qq) in Appendix F and Appendix G up to N 2 LO QCD, respectively. The results beyond N 2 LO can be found in the Online Resource provided with this article. In summary, the coefficient functions for the diagonal terms taking into account the NSV contributions reduce to (24) with the replacement of S d,aa byS d,aa . Note that the modified differential soft-collinear operatorS d,aa contain additional NSV terms inΦ d,aa andΓ aa and we need to keep both SV and NSV terms whenC operator acts on the functions like exponential and logarithms in (16), (39) and (40). We find that the NSV terms being proportional to D k (z i ) log l (1 − z j ), k, l ≥ 0, i, j = 1, 2 and log k (1−z i )δ (1−z j ), k ≥ 0, i, j = 1, 2, do not alter the δ (1 − z 1 )δ (1 − z 2 ), D k (z 1 )D l (z 2 ) and D k (z i )δ (1 − z j ), i, j = 1, 2, k ≥ 0 terms in the coefficient functions. This results in with Φ res d,I = δ (z 1 ) 2 q 2 z 2 µ 2 F dλ 2 λ 2 P I a s (λ 2 ), z 2 + Q I d a s (q 2 2 ), z 2 + + 1 4 1 z 1 P I (a s (q 2 12 ), z 2 ) + 2L I (a s (q 2 12 ), z 2 ) + q 2 d dq 2 Q I d (a s (q 2 12 ), z 2 ) + 2ϕ f d,I (a s (q 2 12 ), z 2 ) where q 2 l = q 2 (1 − z l ) and q 2 12 = q 2 z 1 z 2 . The function P I in (42) is given as with A I is the cusp anomalous dimension and L I (a s , z l ) ≡ C I (a s ) log(z l ) + D I (a s ). The subscript + denotes plus distribution. The function Q I d in (42) is given as The constants C I and D I can be obtained from the the splitting functions given in [53,54] and are known to three loops in QCD. The finite function ϕ f d,I depends on a s and z l , l = 1, 2 and we use explicit fixed order results to parametrize in the following way, We determine the upper limit on the sum over k by studying the dimensionally regularised Feynman integrals that contribute partonic cross sections in fixed order perturbation theory. We know that z space result in the case of SV part [7,16,19] can predict all the distributions to dCFs to all orders starting from a n s provided the lower order results up to a n−1 s are known. In the present case, the inclusion of NSV part can predict terms of the form δ (z l ) log k z j , n + 1 ≤ k ≤ 2n − 1 and D i (z l ) log k z j for i, k = 0, 1, · · · , n; i + k < 2n − 1 at every order a n s provided the form factors, differential softcollinear and mass factorization kernels are known to order a n−1 s . Using the second order inclusive results, some of the authors obtained the third order NSV contributions to dCFs, for the first time in [34], that are contributing to Drell-Yan process and Higgs boson productions. Further the rapidity NSV coefficients ϕ I,(k) d,i and dCFs to third order are obtained in [35] using third order inclusive results in [34] and also with the use of (18). The complete third order results for the Higgs production in gluon fusion are already known in [4,15]. For the DY, the third order prediction of [35] using the same formalism for 2 → 1, is in complete agreement with the [15] for terms of the type D i (z l ) log j (z m ), i, j ≥ 0, l, m = 1, 2. The remaining δ (z l ) log j (z m ) terms in DY and the complete NSV predictions for Higgs production in bottom quark annihilation channel at third order are new and can be found in [35]. The generic expression for the NSV rapidity coefficients ϕ where β i are the coefficients of QCD-β function which are known to five loops [70][71][72][73]. The SV coefficientsG I,(k) j from SV part of soft-collinear distribution as well as the NSV anomalous dimensions C I and D I for I = g, q up to third or-der are provided in the Online Resource supplied with this article. The NSV inclusive coefficients ϕ (k) I, j for I = g, q up to third order can be found in [34]. For completeness, we also provide the explicit results of ϕ I,(k) d, j for Higgs production and Drell-Yan process up to third order along with the partial fourth order results in the Online Resource.
The z space result for coefficient functions expressed in the integral representation in (42) is suitable for studying large N behaviour of the rapidity distribution and hence we use the modified differential soft-collinear operator to obtain the resummed result in N space taking into account the NSV terms. The Mellin moment of dCFs is found to bẽ where the exponentG I,sv+nsv d,N can be organized as a resummed perturbation series in Mellin space as, with h I d,0 (ω, N l ) = h I d,00 (ω) + h I d,01 (ω) log N l , The NSV resummation coefficients are g I d,i and h I d,i . The coefficient g I d,1 is found to be zero. The coefficients g I d,i+2 depend on universal cusp anomalous dimension A I and D I d , while h I d,i s are determined by the NSV coefficients ϕ f d,I as well as by C I , D I from P I (a s , z l ) as given in (50). The resummation coefficients g I,i d,0 , g I d,i (ω), g I d,i (ω) and h I d,i (ω) contain leading, next-to-leading, · · · , SV and NSV logarithms in the dCFs.
Rescaling the constants by β 0 as A we present below the results of the NSV exponents g I d,i (ω) and h I d,i j (ω) after setting µ 2 R = µ 2 F = q 2 . The full list of these exponents with explicit dependence on µ 2 R and µ 2 F are provided in the Online Resource supplied with this article.
In summary, thanks to the remarkable simplification of mass factorised formula for the rapidity distribution for diagonal channels and the knowledge of logarithmic structure of differential soft-collinear distribution and the mass factorization kernels from fixed order results, we can systematically resum both SV and NSV terms in z as well as in double Mellin space. Note that we could extend the framework that was used for the production of a single colorless state to the state of n-colorless particles. This was possible simply because of the fact that IR structures of the building blocks, namely, the pure virtual contribution, soft-collinear distributions and the mass factorization kernels for 2 → n is identical to those for 2 → 1. In addition, the IR singularity structure of these quantities in dimensional regularisation plays an important role to understand the leading SV distributions and sub-leading NSV logarithms in the real emission processes. In the case of SV, there exist two universal or process independent soft-collinear functions, namely for quark-anti quark (qq) and gluon-gluon (gg) initiated processes. For the NSV, we find that a part of the soft-collinear functions is partially determined by universal anomalous dimensions and the remaining part depends on the underlying hard process. However, for a given process, the resummed results either in z or {N 1 , N 2 } can be used to predict SV and NSV terms at higher orders provided the lower order results are known. With this in mind, we have derived a general formula for coefficient functions up to fourth order and resummed results to N 3 LL accuracy which will be useful when the complete four loop form factor, soft-collinear functions and lower order process dependent terms become available.

Conclusions and Outlook
In today's era reducing the theoretical uncertainties remain one of the main motivations for higher order radiative corrections. It is also particularly relevant for constraining beyond SM scenarios and validation of the SM itself. Among several observables, differential cross-sections allow a wider range of comparisons with the experiment and hence several attempts were made in the past for better theoretical understanding of the same. In this article, we restrict ourselves to the discussion of differential rapidity distribution for the production of n-number of colorless particles in the hadronic collision within the realm of perturbative QCD.
Through this publication, we intend to present a systematic framework for the study of soft-plus-virtual corrections to the differential distribution with respect to the rapidity variable y, for the production of n-colorless particles in the hadron collider. The infrared structure of rapidity distribution which was earlier studied in ref. [7] for Sudakov type processes is further extended to the case of 2 → n scattering. We employ the universality of the soft enhancements associated with the real emission diagrams. The main deviation from the Sudakov type formalism comes from the virtual corrections where the kinematic dependence is much more involved and hence these are now expressed in terms of scalar products of the kind {p j · q k }. The rest of the formalism relies on the collinear factorization of the differential cross section, the renormalization group invariance, universality of perturbative infrared structure of the scattering amplitudes, and the process independence of the soft-collinear distribution. Besides this, we also use an additional fact that the N-th Mellin moment of the differential distribution has a relation with its inclusive counterpart in the limit N → ∞, as depicted through (18). The mere use of this fact enables us to to get an all order relation between the soft-collinear distribution of inclusive cross-section and that of rapidity. Hence from the given quantity in inclusive part, we can determine it for the rapidity and thereby avoid performing the explicit computation of the real emission processes for rapidity distribution. The goal of this current article is to present the general structure for the SV differential rapidity distribution up to N 4 LO and also the resummed predictions till N 3 LL level in QCD, which can be expressed in terms of universal anomalous dimensions along with the process dependent virtual matrix elements. The former, which comprises of process independent finite segments of soft-collinear distribution and the mass factorized kernels, remains unaltered irrespective of the number of colorless particles in the final states. Furthermore, the soft-collinear distributions for the quark and gluon initiated processes are found to be related to each other through simple quadratic Casimir scaling, known as the maximally non-Abelian property.This is explicitly verified up to N 3 LO. However, whether the validity will remain intact beyond this order with generalized Casimir scaling, that will be an interesting thing to look into. Often, one finds that in certain kinematic regions, the subleading logarithms, namely NSV terms can not be ignored in phenomenological studies. Our investigation on these terms for diagonal partonic channels reveals that there are similarities in the structure of IR terms with those of 2 → 1 process allowing us to propose resummed predictions for NSV terms within the same framework.
In summary, in order to obtain the fixed order as well as resummed prediction for the differential rapidity distributions of a generic n-colorless final states, one merely requires the form factor corresponding to the hard process under study provided the soft-collinear distribution for Sudakov type process is known. We present the analytical re-sults for the fixed order up to N 4 LO and the resummed predictions up to N 3 LL level in the appendix for the scale choice of µ 2 R = µ 2 F = q 2 and the same with the explicit scale dependence are provided in the Online Resource supplied with this article.

Acknowledgements
The algebraic computations have been done with the latest version of the symbolic manipulation system FORM [74,75]. We would like to thank S. Catani

Appendix A: Soft-collinear distribution for rapidity distribution
In this section, we present soft-collinear distribution S d,I , as defined in (23), in powers of a s (µ 2 R ) up to N 2 LO. Expanding the quantity in powers of a s as we present the results for µ 2 R = µ 2 F = q 2 . The results up to N 4 LO with explicit scale dependence can be found from the Online Resource supplied with this article.
The symbolsG ,I, j d,k that appear in the aforementioned softcollinear distributions are provided explicitly in the Online Resource with this article.

Appendix B: Soft-virtual partonic rapidity distribution
We expand the ∆ sv d,I , as defined in (24), in powers of a s (µ 2 R ) through Here, we present ∆ sv d,I to N 2 LO for the specific scale choice The results up to N 4 LO with explicit dependence on µ R and µ F are provided in the Online Resource supplied with this article.
In the aforementioned equations, we define M (m,n) where |M (n) I is the UV renormalized pure virtual amplitude at n-th order in a s as introduced in (13).

Appendix C: Soft-collinear distribution for threshold resummation
The universal soft-collinear operator that is required to obtain the resummed cross-section in z-space, defined in (29), is expanded in powers of a s (µ 2 R ) as We present the results for q 2 = µ 2 R = µ 2 F below to fourth order in coupling constant. The result with explicit scale dependence can be obtained from the Online Resource provided with this article.  In this section, we present the explicit results of ∆ sv d,I , as defined in (24) and (F.16), for the Drell-Yan (I = q), and the Higgs boson productions through gluon fusion (I = g) as well as bottom quark annihilation (I = b) at fourth order in coupling constant. Setting µ 2 F = µ 2 R = q 2 , in the following, we provide only the new results, and the old results for Drell-Yan and Higgs boson productions can be found in [7,16,19]. The results with explicit dependence on µ R and µ F are provided up to N 4 LO in the Online Resource files supplied with this article.  ,C 2 F C 2 A ,C 3 F C A ,C 4 F are the unknown coefficients of the color factors in four loop soft and collinear anomalous dimensions. In the aforementioned equations, n f v is proportional to the charge weighted sum of the quark flavours and N 4 = (n 2 c − 4)/n c [67]. Following [55], we have
(G. 19) The results up to N 4 LO with explicit scale dependence are provided in the Online Resource supplied with this article.