Unitarity methods for Mellin moments of Drell-Yan cross sections

We develop a method for computing Mellin moments of single inclusive cross sections such as Drell-Yan production directly from forward scattering diagrams, by invoking unitarity in the form of cutting equations. We provide a diagram-independent prescription for eliminating contributions from unwanted cuts at the level of an expansion in the reciprocal ω = 1/z variable. The modified sum over powers of ω produces the result from physical cuts only, with the nth coefficient precisely equal to the nth Mellin moment of the cross section. We demonstrate and validate our method for representative one- and two-loop diagrams.


Introduction
The efficient computation for higher order QCD corrections for scattering processes has been a mainstay of research in theoretical particle physics for decades, as this directly impacts the potential for discovery at colliders: the more precisely signal and background are computed, the more significant comparison of theory with data can be.In recent years a notable increase in computational capacity has taken place, spurred largely by the development of unitarity techniques for computing scattering amplitudes.The ability to determine a scattering amplitude from its poles and branch cuts [1][2][3] has been a watershed in these efforts, especially for one-loop high-multiplicity processes.
Unitarity methods for few external legs but at higher loop have proven to be highly valuable as well.Reverse unitarity techniques have been important in relating real emission amplitudes to virtual ones [4] for two and three-loop calculations.An inspiration for the present paper is the body of work on the computation of the 2-and 3-loop splitting functions and deepinelastic scattering (DIS) Wilson coefficients [5][6][7][8][9][10].Computing the DIS structure functions by moments has been a very succesful approach for massless partons [5,8,[11][12][13][14], and also for heavy quarks [15][16][17][18][19][20][21][22][23].For the inclusive DIS structure functions it is possible to use the optical theorem to compute all the cuts of the forward scattering process, and directly extract the Mellin moments of the coefficients.In essence, from unitarity considerations, one may expand the forward scattering amplitude in reciprocal powers of the Bjorken scaling variable x, the coefficient at order n then being the nth Mellin moment of the coefficient function.
In this paper we aim to generalise this method to a single-inclusive cross section, specifically the Drell-Yan cross section, which is then prototypical for other processes such as Higgs production in the large top mass limit, for which remarkable results for high-order corrections to Higgs have recently been achieved [24][25][26][27][28][29].The optical theorem does not directly apply to the Drell-Yan cross section, not being a fully inclusive observable.We show in this paper that it is nevertheless possible to compute the Mellin moments of the Drell-Yan cross section directly from forward diagrams, using unitarity and an expansion in reciprocal powers of z.The key aspect of our method is the efficient subtraction of unwanted cuts, through complex-valued shifts of the moment variable n and through a replacement rule for harmonic sums.
The paper is organised as follows.In section 2 we review the relevant aspects of unitarity and the optical theorem, and their role in DIS.In section 3 we treat the one-loop corrections in "scalar" Drell-Yan correction with our new method in some detail, highlighting key features.In section 4 we test our method at two loops, showing how our methods work for representative two-loop forward scattering scalar diagrams.Here we show explicitly how remove contributions from unphysical cuts, such that those from physical cuts are unaltered.We conclude with a summary and some remarks towards further development.

Forward amplitudes and unitarity
In this section we outline the general ideas of the paper, postponing technical details to the following sections, where one and two loop examples will be discussed.
We start in section 2.1 reviewing the essential points that make the optical theorem successful for DIS.Then, in section 2.2 we move to the Drell-Yan case, stressing the differences that make a generalization of the DIS method highly non-trivial.Among these, the most problematic one is the presence of unphysical cuts, absent in DIS, that need to be removed from the discontinuity of the forward amplitude.Therefore in section 2.3 we classify all unphysical cuts, showing that most of them either vanish or are easily treated.For the remaining unphysical cuts, we outline a solution in section 2.4, referring the reader to section 3 and section 4 for more technical explanations.

DIS and the optical theorem
Let us first review the role of the unitarity, in the form of the optical theorem, in deep-inelastic scattering through off-shell photon exchange.Our exposition follows largely that of ref. [30].It is well-known that the fully inclusive cross section for this process, e(l) + P (p) −→ e(l ) + X, can be written in the form with q = l − l , and with W µν (L µν ) the hadronic (leptonic) tensor.They are defined as where implicit spin quantum numbers in the external states are summed over.Note that the sum over final states |n is fully inclusive in terms of QCD as both explicit momenta p and q are incoming.Current conversation and parity invariance in both indices then imply the structure W µν (p, q) = − η µν − q µ q ν q 2 W 1 (x, Q 2 ) + p µ − q µ p • q q 2 p ν − q ν p • q q 2 W 2 (x, Q 2 ) , (2.3) so that the structure of the proton is encoded into two scalar functions that depend on the variables (2.4) The optical theorem applies to the hadronic tensor W µν , since the sum in eq.(2.2) is fully inclusive so that we can write W µν (p, q) = 2 Im T µν (p, q) , ( where T µν is the forward Compton amplitude γ * (q) + P (p) −→ γ * (q) + P (p), having the same tensor structure as in eq. ( 2.3) but now in terms of the scalar functions Note that we have chosen to indicate the functional dependence in terms of the reciprocal x variable, for reasons we discuss below.The functions W i and T i both have cuts starting at branch points at x = ±1, corresponding to the kinematical conditions for normal thresholds, (p ± q) 2 > 0. For the W i (x, Q 2 ) functions, the cut then runs from x = −1 to x = 1.For the T i (ω, Q 2 ) functions, consequently, the cuts lie along the ω-intervals (−∞, −1] and [1, ∞).These are also the only cuts, and we have in general for the T i Figure 1: Branch cut structure of T i (ω, Q 2 ) with the two contours used for T (n) i (Q 2 ).On the left the contour C 0 wraps around the origin, while on the right the contour C 1 encloses the two branch cuts.Note that the combination ω −n−1 T (ω) of eq.(2.8) has an additional pole at the origin.
One can now compute Mellin moments of the W i (x, Q 2 ) functions by expanding the T i (ω, Q 2 ) amplitudes.The nth derivative of T i at ω = 0 may be rewritten by Cauchy's theorem in terms of the contour in fig.1a T The contour C 0 may be deformed into the contour C 1 shown in fig.1b.Then, using eq.(2.7), we have where the discontinuity of a function in the variable x is defined in general by (2.10) The presence of the factor (1 + (−1) n ) in eq.(2.9) implies that odd series coefficients vanish.
For n even, instead, using the optical theorem in the form of eq.(2.5), and changing integration variables to x = 1/ω, we get where the second equality defines the Mellin transform M n .Thus, indeed, the expansion of the forward scattering amplitude in ω yields the Mellin moments of the cross section.
A few remarks are in order.This way of using the optical theorem, computing directly the Mellin moments of the DIS structure functions by expansion in 1/x, has been marvellously successful for 2-and 3-loop calculations for DIS [5][6][7].Translating back to momentum space is readily done, and produces known combinations of functions (Harmonic Polylogarithms (HPL's) [31]).The presence of the branch points at ω = ±1 and the analytical behaviour of the T i (ω, Q 2 ) near ω = 0 is very helpful towards the consistency and also practicality of the method.
The question whether the DIS method can be generalized to semi-inclusive cross sections such as Drell-Yan and Higgs production will be addressed in the next sections, and is indeed the central issue of the present paper.

Analytical structure of one-particle inclusive processes
We consider the inclusive production of an electroweak boson V of invariant mass Q 2 by quarkantiquark annihilation where X represents any partonic contribution to the final state and the vector boson V may be an off-shell photon γ * , or an on-shell W ± or Z boson.As such, this process is described by two scales only, the mass Q 2 and the squared partonic center-of-mass energy s, from which it is possible to define the dimensionless variables where z is the variable analogous to the Bjorken variable x for DIS.In the following we will focus on the case of an off-shell photon, referring to this process as Drell-Yan.Apart from interesting in its own right, this process is prototypical for many other partonic processes relevant at hadron colliders, especially Higgs boson production via gluon fusion.Differences with the Drell-Yan case resides only in numerator factors which go along for the ride in our method.Before discussing how the generalisation of the DIS formalism to the Drell-Yan case can be set up, let us review those similarities and differences between the two processes that are relevant for our purposes.
At face value, the differences seem not very large.Focusing on the partonic part, the set of diagrams for the Drell-Yan process can be obtained from the DIS ones, by crossing the exchanged off-shell photon to the final state, and the outgoing quark to the initial state.However, this crossing has significant consequences.First, the off-shellness of the photon becomes time-like and can be effectively regarded as a mass.Therefore the forward amplitude q(p) + q(p) −→ q(p) + q(p) will contain a massive propagator.Then, most importantly, the vector boson must be present in the final state.Hence, the process is not fully inclusive like DIS, but only singleparticle inclusive, so that the optical theorem, as the simplest realization of unitarity, cannot be used and the cross section is not given by the imaginary part of the full forward amplitude.
A further complication arises when moving to Mellin space.Looking at the analytical structure of the DIS forward amplitude, branch cuts in the ω-plane are at (−∞, −1] and [1, ∞).Due to the symmetry in eq.(2.7) of the forward amplitude, one may consider only the cut along the positive real axis, which can be eventually converted to a Mellin transform, as discussed in the previous subsection.In the Drell-Yan case instead the forward amplitude generally will have more branch cuts, in particular also along (−∞, 0] and [0, ∞) and no symmetry relates the forward amplitudes with opposite value of ω.A new strategy is needed if we want to extract the series coefficients of the forward amplitude expanding it around the branch point ω = 0.
These considerations suggest that extending the DIS techniques for directly computing Mellin moments to the Drell-Yan case is not straightforward.However, we shall see that it is possible when using unitarity cuts.Let us discuss the key aspects of this in somewhat more detail.
The optical theorem relates a cross section to the imaginary part of the relevant forward amplitude.At the same time, this imaginary part is, by the Cutkosky rules [32], equal to the sum over all cuts of the amplitude.For a fully inclusive processes like DIS, these cuts correspond to the phase space integration of the squared matrix elements of the process.For Drell-Yan instead this is not the case, as cuts that do not cut the massive photon are not to be included for the cross-section.However, the use of unitarity cuts is considerably more general, and holds on diagram-by-diagram basis.In general, branch cut discontinuities in different channels correspond to different cuts of the diagram [33,34] for any Feynman diagram F. Our approach exploits this fact fully when F is a forward scattering diagram.
Our goal is to compute (moments of) the cross section from the forward amplitude.The cross section can be reconstructed from the discontinuity of the forward amplitude across the physical branch cut.In general an amplitude has discontinuities around unphysical branch cuts as well, and these must be subtracted.This does not seem a very efficient procedure, as it apparently requires one to compute unphysical-cut diagrams nonetheless.Moreover, the unphysical cuts may be even more complicated than the physical cuts.However, as we will see in the following sections, one can modify the analytic structure of the forward amplitude such that its discontinuity is given by the sum of physical cuts only.In particular, we will see that, after moving to Mellin space, it is possible to automatically select the physical cut without the need to subtract (and compute) the unphysical cuts.
Let us first review the classification of the cuts appearing in the forward amplitude of the Drell-Yan process.

Classification of Drell-Yan cuts
The set of diagrams we would like to discuss are those required in a NNLO calculation, though many of the features will be valid also at higher order.To set up our classification, we regard forward diagrams as amplitudes that may depend on different channels, and therefore can be cut in all possible ways, as long as the diagram is cut into two (connected) subgraphs.In this regard, we distinguish four different classes of cuts, depicted in fig. 2 and denote them as vertex, s-channel, t-channel and u-channel cuts.Of course, for a forward amplitude and with on-shell external lines the only possible invariant is s, but the nomenclature will be useful, and is based on the case when final momenta are different from the initial ones, such that also the t and u channels would be open.
The s-channel cuts are the only ones that can be interpreted as a phase space integration of squared matrix elements.Among these, physical cuts, i.e. those that contribute to the crosssection, are only those s-channel cuts that pass through the massive photon, and we thus call massive s-channel cuts.More generally, we call massive (massless) every cut that does (does not) cut the massive boson propagator.At first it seems that the number of unphysical cuts might grow dramatically with the order of the computation, making it difficult to control them.However, a number of simplifications are possible, making some of these cuts give a vanishing contribution.
The vertex cut in fig.2a vanishes, because it measures the discontinuity of the forward amplitude in the p 2 -channel.But this discontinuity is zero, because the forward amplitude does not actually depend on this variable due to the on-shell condition p 2 = 0.The same holds when any other of the four vertices is cut.Furthermore, by the same token, the t-channel cut in fig.2c vanishes as well.This leaves only the s-and u-channel cuts to be considered: the massless s-channel cuts and the (massive and massless) u-channel cuts.In section 4.1 we will treat these unphysical cut diagrams in greater detail and in specific examples.In the following subsection we first review the general ideas how to deal with these cuts.

Extracting the physical cuts from the forward amplitude
At this point we make an important observation.The forward amplitude carries more information than needed; indeed we are only interested in its imaginary part.We have the freedom to modify the amplitude, as long as the branch cut structure remains the same.For instance, adding a constant or even an analytic function will not affect the physical information one wishes to extract from its cuts.This consideration leads us to disregard lower order Mellin moments.Indeed, assuming that the forward amplitude can be expanded around ω = 0 as its series coefficients c n will be defined only for n ≥ n 0 .However, any shift in n 0 making the sum start from a new positive integer is equivalent to adding to the original f (ω) simply polynomials in ω, which does not affect the branch cut structure.Therefore, we conclude that no physical information is carried in the lower bound of the sum, and henceforth we shall omit it in series expansions except where necessary.We can even take a further step in this line of reasoning.Since we are interested in extracting the discontinuity of the forward amplitude only across the physical branch cut, we have the freedom to redefine the forward amplitude, modifying also its branch cut structure, as long as this leaves the physical branch cut unaltered.Also, poles in ω can be removed from (2.15) because such poles do not correspond to physical cuts.These steps form the essence of the strategy we shall implement to deal with the unphysical cuts.We start with the first two types of cuts presented in section 2.3: massless s-channel cuts and massless u-channel cuts.These classes of cut diagrams contain factors ω , where = 4−d 2 is the dimensional regulator.Hence, they belong to branch cuts starting at ω = 0.As such, those are the cuts that prevent the forward amplitude to be written as a power series in ω around ω = 0. We introduce a shifting procedure, through which it will be possible to define a new function f (ω) with no branch point at ω = 0.This will be necessary already at one loop and will be further discussed in section 3.4.
The most difficult kind of unphysical cut appears only from two loops: the massive u-channel cut, because it corresponds to a branch cut in the forward amplitude starting at ω = −1.We shall go beyond the simple shifting procedure and subtract the contributions from this cut directly in Mellin moment space, following the extending reasoning above.We are able to compose a dictionary of replacements for harmonic sums, which may be applied to any diagram.This procedure will be first discussed in section 4.2 and is applied to a two-loop crossed box in section 4.3.2.
After all unphysical cuts are removed, one can repeat the procedure carried out for DIS in fig. 1.This time only the branch cut [1, ∞) is present and the generalization of eq.(2.11) (for both even and odd Mellin moments) reads where F (n) phys are the series coefficients of the forward amplitude modified such that it contains only the physical cut.Let us now turn to a more explicit illustration of these methods at one-loop.

Drell-Yan at one loop
Here we develop the concepts from section 2 for the one-loop Drell-Yan cross section.The standard approach requires the evaluation of the matrix elements for real and virtual corrections, and the subsequent integration over the phase space.This requires the computation of three phase space integrals, represented in fig.3, which in the language of section 2 are massive schannel cuts.In order to minimise technical complications we will henceforth omit numerator factors in those diagrams, i.e. we restrict to the scalar case, as these are anyway irrelevant for illustrating our method.This is motivated by the fact that the analytical structure of Feynman integrals in QCD arises in essence from denominators in diagrams. 1n the next three subsections the calculation of the scalar equivalent of the diagrams in fig. 3 is presented, verifying the cutting equation in Mellin space in the form of eq.(2.16).In this process we encounter two relatively harmless types of unphysical cuts and we show how to deal with them as discussed in section 2. In particular, we devote one subsection to the notion of the shifting procedure, which is needed to remove both kinds of cuts.Given the simplicity of those one-loop calculations, the computation is actually performed to all orders in .However, extracting the series coefficients of the ω expansion exactly in is not feasible at higher loops.Therefore, in the last subsection section we illustrate how to compute the Mellin moments using IBP identities, by the example of the one-loop box diagram.

Triangle diagram
The simplest of the three diagrams that contribute to the one-loop Drell-Yan forward amplitude involves a triangle loop.The cutting equation for this graph is schematically depicted in fig.4, showing both a physical and an unphysical cut.It appears that in order to compute the physical Cut 1 T contribution one would need to 'subtract' the unphysical Cut 2 T contribution from the full discontinuity Disc T .
To verify explicitly the presence of an unphysical cut, let us compute both the forward amplitude and the physical cut diagram.The forward amplitude T reads where the prefactor coming from the loop integration is given by and we have set µ2 = Q 2 .We rescale forward diagrams F by their mass dimension s dim [F ] , such that they become dimensionless.Applying the Cutkosky cutting rule, the physical cut reads where we expressed the result as function of z = 1/ω.Clearly, the discontinuity of the forward amplitude is not given by the physical cut alone, as can be seen by the presence of (−ω) − in the forward amplitude. 2Indeed, upon expanding this factor in , it is evident that T has a branch cut for ω < 0, whereas Cut 1 T is different from zero only at ω = 1.Therefore, to recover the physical cut diagram, we should compute also the unphysical Cut 2 T .As we shall see, this can be bypassed in Mellin space.We start computing the Mellin moments of the physical cut where the phase e −iπ is due to the minus sign in (−z) and can be fixed by keeping track of the Feynman iη in the propagators.Then, we would like to compare eq.(3.4) with the series coefficients c n of the forward amplitude expanded in powers of ω.However, as can be seen in eq.(3.1), this expansion cannot be perfomed, since T contains a non-integer power of ω.This is of course expected: the forward amplitude has a branch point at ω = 0 and therefore cannot be expanded around that point.However, the structure of this branch cut starting from the origin is simply given by (−ω) − .Therefore, we apply the following procedure.We write a series representation for the other factors, to get Then, we shift n → n + in the summand, leaving the lower bound of the sum unaltered.This defines a new function T with no branch cut from the origin, but with the physical pole in ω = 1 still present.This procedure will be also used in section 3.3 for the crossed box diagram, where a non-trivial ndependence of the summand will make the shift less trivial.Thus, writing T = n c n ω n , we find Comparing eq.(3.4) and eq.(3.7) we get which verifies the cutting equation in Mellin space in the form of eq.(2.16).The triangle diagram is the easiest example that exhibits an unphysical cut and where it is possible to test our Mellin space approach.In the next subsection we discuss an example with non-trivial n-dependence.As we did for the triangle diagram, let us compute explicitly the forward amplitude B 1 and the cut diagram Cut 1 B 1 .From a direct computation we have

Box diagram
where we expressed the result as a function of ω = s/Q 2 .The cut diagram Cut 1 B 1 is easily computed as well.It is defined as As for all massive s-channel cuts, this integral is non-vanishing when s > Q 2 , and we can perform the computation in the centre-of-mass frame, where Computing the integral in this frame yields where we expressed the result as a function of z = 1/ω.Note that both results (3.9) and (3.12) are valid to all orders in .We can now verify the cutting equation in Mellin space.Specifically, the forward amplitude may be written as a series representation B 1 = n c n ω n , where while the Mellin transform of the cut diagram can be trivially computed and reads Comparing eq.(3.13) and eq.(3.14) we find which verifies the cutting equation in Mellin space in the form of eq.(2.16).This example has shown the cutting equation in Mellin space with a non-trivial n-dependence.To increase further the complexity, in the next section we will compute the crossed-box diagram, which has both non-vanishing unphysical cuts and non-constant moments.

Crossed-box diagram
The crossed box B 2 is the last contribution to the one-loop Drell-Yan forward amplitude.It results from the normal box B 1 by interchanging the two final state momenta.This has the consequence that the u-channel cut does not vanish so that the discontinuity of the forward consists of two cuts, a physical and an unphysical cut, as shown in fig.6 As we did for the other two examples, we compute the forward amplitude and the cut diagrams.The former it is defined as Again, applying standard techniques with Feynman parameters, it is possible to work out a result for this integral exact in , which is The computation of the physical cut Cut 1 B 2 (see fig. 6) is essentially the same as the cut diagram of the normal box B 1 , since it again has a two-particle phase space non-vanishing for s > Q 2 .It reads In order to clarify the structure of the cutting equation, we explicitly compute also the unphysical u-channel cut diagrams (the second diagram on the right-hand side of fig.6).This is defined as In contrast with the s-channel cut, this is non-vanishing when s < 0, so we perform the computation in the following frame: Putting the right-most vertical propagator on-shell fixes as shown in fig.6, proving that in z-space both cuts are needed to reproduce the discontinuity.Hence, in order to work out the physical cut contribution, one would need to subtract the unphysical cut from the discontinuity of the forward diagram.As we shall see, this can be bypassed in Mellin space by using the shifting procedure that we introduced in section 3.1 for the triangle diagram.We start writing eq.(3.17) as a series representation for the hypergeometric functions and moving the overall ω inside the sums.This gives We note that the last term contains a non-integer power of ω, which prevents us from constructing a formula for the series coefficients of the forward amplitude.Therefore, we apply to this term the same trick that we used for the triangle diagram.We change n → n + in the summand but not in the range of the sum (i.e.we sum from n = 1).This gives for the last term We have now defined a new function B 2 with no branch cut starting at the origin, from which we can extract the series coefficients.Applying this prescription to eq. (3.23), we see that last term cancels against the first sum and we are left with the second term.In conclusion, writing B 2 = n c n ω n , we have Now we move to the cut diagrams.As for the normal box B 1 , they can be computed after writing them as function of z = Q 2 /s, and then performing the standard Mellin transform.Using the results in eqs.(3.18) and (3.21) this yields The latter moments are zero due to the step function θ(−z).Comparing the series coefficients in eq.(3.25) to the moments of the physical cut in the first line of eq.(3.26), we see that which verifies eq.(2.16).We conclude that the Mellin moments of the physical cut are indeed reproduced by the series coefficients of the modified forward amplitude B 2 and that the unphysical cut has been removed by the shifting procedure.

Shifting procedure
For the triangle and the crossed-box diagrams, we introduced a prescription to deal with a forward amplitude f (ω) that cannot be expanded around ω = 0. We also saw that for the crossed box diagram this is due the presence of a non-vanishing massless u-channel cut, which encodes the part of the forward amplitude having a branch cut along the negative real axis in the ω-plane.For the triangle diagram, instead, this is due to a non-vanishing massless s-channel cut, with branch cut in the positive real axis in the ω-plane.
In order to clarify the shifting procedure let us review its general features.Assume that (a non-analytic piece of) the forward amplitude f (ω) can be written as where k is some integer and g(ω) is analytic in ω = 0.The following discussion will trivially generalise to the case with more terms, which for instance might have different values of k or different signs like (±ω) k .The functions f (ω) and g(ω) may also depend on , which is left implicit for brevity.Writing g(ω) as an expansion around ω = 0, we have The shifting procedure is defined by replacing n with n − k in the summand, but not in the lower bound of the sum. 3 This produces a new function This function f (ω) is a modified version of the forward diagram that is precisely what is needed for our purpose, if the following two criteria are met: (i) The unphysical cut must be absent in Disc ω f (ω); (ii) The discontinuity around the physical cut must be unaffected.
We discuss the validity of these two conditions in turn.
The first criterium goes back to the assumption made in eq.(3.28), namely that the nonanalyticity of f (ω) around ω = 0 is captured by an overall factor ω k .This can be argued with dimensional analysis, when looking at the physical complex s-plane (remembering that ω = s/Q 2 ).Branch cuts starting at s = 0 are described by the single dimension-full quantity s, irrespective of the value of Q 2 .Since Feynman diagrams have a fixed integer mass dimension, the only way in which s can occur is as an overall power of s and not as the argument of some other (elementary) function.Fractional powers are excluded by dimensional analysis.The only deviation from integer s powers allowed is due to the d-dimensional integration measure.Feynman integrals yield results proportional to s k which, combined with the dimensional 3 An alternative definition of the shifting procedure is the following.First rewrite the forward amplitude as f (ω) = ∞ n=n 0 +k c n−k ω n , where sums starting at non-integer lower bound α ∈ C are to be interpreted as The shifting procedure may then alternatively be defined as setting = 0 in the lower bound of the sum, so that f (ω regularisation mass scale set to Q 2 , produces overall factors ω k .We thus conclude that any unphysical cut starting at ω = s = 0 is captured by functions of the form in eq.(3.29).The modified forward amplitude in eq.(3.30) is analytic at ω = 0 by construction and therefore the unphysical cut is indeed completely removed from Disc f (ω) by the shifting procedure.
The second criterium ensures that altering the forward amplitude, does not destroy the connection between the discontinuity around the physical branch of the forward amplitude and the sum of physical cut diagrams.This issue can be clarified through some toy examples.Let us first consider a simple case where the "forward amplitude" is given by Applying the shifting procedure, we arrive at the new function While the original function has its branch cut along the negative real axis removed, both functions have the same pole structure in the region ω ≥ 1, namely Therefore in this example the discontinuity in the physical region is unaltered by the shifting procedure.Another example is given by a function with a branch cut, rather than a simple pole, in the physical region.This mimics more closely the cases we encountered at one loop.Consider The shifting procedure produces Again, these functions have the same discontinuity in the physical region.This is most easily verified upon writing both functions as an expansion in , Using the identities Disc one readily confirms that the discontinuity of f 2 (ω) in the physical region is equal to that of f 2 (ω) and is thus unaltered by the shifting procedure.

Direct extraction of series coefficients from IBP's
The shifting procedure from the previous subsection requires the result of an integral to be given in terms of ω − in unexpanded form.Indeed, if the integral were expanded in , then the logarithmic divergence at ω = 0 could no longer be removed by shifting n.In the one-loop examples of the previous section, there is no problem since the forward amplitude diagrams are exact in .Moreover, for each of the one-loop diagrams a simple series representation is known, which allows us to extract their series coefficients exact in as well.At higher loops it is not realistic to expect exact results in for all the forward diagrams.However, it is in fact sufficient that the divergent part of a forward diagram f (ω, ) around ω = 0 is written as ω − g(ω, ), where g(ω, ) is regular at ω = 0, and may also be given in expanded form g(ω, ) = k≥k 0 k g k (ω).Such a hybrid expression can be obtained by making a series ansatz for the forward diagram.Specifically, for one-loop diagrams one writes where dim[f ] denotes the integer mass dimension of the forward amplitude f .This structure for a forward diagram is not surprising, since in general the function is non-analytic at the origin ω = 0. Without loss of generality, one can decompose such a function into a sum of analytic and non-analytic pieces.As discussed in the previous subsection, the non-analyticity can always be captured by a factor ω − multiplied by another function, which is regular at the origin and thus admits a series representation.The coefficients c n ( ) and d n ( ) may be given exact in or as an expansion in ; in either case the shifting procedure works.
A further benefit of making the series ansatz in eq.(3.38) is that the series coefficients can be extracted more directly by deriving equations for them and subsequently solving the equations.One way to proceed along these lines is to generate a differential equation for f (ω, ) from integration-by-parts (IBP) identities.Inserting the series ansatz into such differential equation yields in turn a difference equation for the series coefficients, which takes the form where the A i,n and F n are rational functions of , whose form depends on the particular differential equation for f (ω, ) under consideration.If r = 1, then eq.(3.39) is simply a recursion, in which case the series coefficients can be found exact in .For diagrams with multiple loops the order of the difference equation becomes typically quite high (we find up to r = 8 for two-loop diagrams).In that case it will be advantageous to expand the difference equation eq.(3.39) in and solve for the coefficients c n ( ) = k≥k 0 k c k,n , order-by-order in .The task of solving the resulting difference equations for the set {c k,n } may be even further simplified by following the approach in ref. [37], which exploits the expectation that these coefficients are given in terms of harmonic numbers.
For the computation of the two-loop diagram in this paper we indeed adopt the approach in ref. [37] and seek solutions to eq. (3.39) of the form for reasonable choices of k, , m.The functions S (n) are harmonic sums, whose properties are well-known [38].The unknown coefficients A k, ,m contain both rational and transcendental numbers, and may be determined as follows.The simplest approach is to evaluate the difference equation for c k,n at suitably many values of n, so as to obtain a system of equations which may be solved for the unknown A k, ,m .In a more sophisticated method [37] each term in the difference equation for c k,n is projected onto a basis of synchronised harmonic numbers, after which the coefficients of each harmonic number is equated to zero.This also yields a system of equations for the unknown A k, ,m .We have implemented both techniques and successfully applied them to the two-loop examples in section 4.3.
Before closing this section, let us present an example of the methods in this subsection for obtaining the series coefficients from IBP's.To this end, consider the one-loop box B 1 from section 3.2.After shifting the loop momentum k → q = k + p in eq.(3.9), the scalar box integral becomes a special case B(2, 1, 1) of the topology The derivative of B(a, b, c) with respect to Q 2 produces another integral in the topology: Performing IBP reduction on the right-hand side of eq.(3.42) thus yields a differential equation for B(a, b, c) in terms of (typically) simper integrals.In the case of the one-loop box, the differential equation for B(2, 1, 1) reads4 The inhomogeneous terms on the right-hand side are simpler integrals because they have fewer propagators.Inserting the known bubbles B(a, 0, 1) on the right-hand side and the definition of B(2, 1, 1) in terms of the forward box B 1 on the left-hand side, leads to a differential equation for B 1 : (3.44) One way to solve this differential equation is by inserting an ansatz for B 1 , like in eq.(3.38), in terms of unknown coefficients c n ( ) and d n ( ).Upon doing so, one finds that d n ( ) = 0, while the other coefficients c n ( ) satisfy the equation Equating terms with equal powers of ω on both sides leads to a recurrence relation in n, This recursion is solved by This result for the series coefficients of B 1 fully agrees with eq.(3.13).
In this simple example it has been possible to solve the recursion exact in .For higher loop diagrams this is typically not expected to be possible to do exactly but rather order-by-order in .In the next section we extend our investigations to two-loop diagrams, where it is shown how to solve a higher-order difference by making an ansatz for the series coefficients in terms of harmonic numbers.

Drell-Yan at two loop
The previous section discussed for Drell-Yan production at the one-loop level, how Mellin moments of cut diagrams can be computed as series coefficients of forward diagrams.A feature in our approach is that unphysical cuts in forward amplitude diagrams can be removed by a shifting procedure.At higher loops this shifting procedure is no longer sufficient.Indeed, in this section we extend our investigations to two loops, for which we develop an additional prescription to subtract unphysical cuts.Two-loop diagrams serve furthermore as non-trivial applications of our method for direct extraction of Mellin moments from integration-by-parts identities, as discussed in section 3.5.In the next subsection we start by listing all possible types of unphysical cuts, placing particular emphasis on the new type of unphysical cut appearing at two loops.We then describe our methods to remove them, working out two examples in detail.

Unphysical cuts of two-loop diagrams
Let us analyse the possible unphysical cuts of two-loop forward diagrams.At one-loop level there are two types: unphysical s-channel cuts and massless u-channel cuts (which do not cut the massive photon).At two loops, there is the possibility for a new type: massive u-channel cuts (which do cut the massive photon).All types of unphysical cuts are illustrated in fig. 7. We briefly discuss the differences between these types of unphysical cuts.
Massless s-channel cuts A cut of this type appears already in the case of the one-loop triangle in fig. 4. A two-loop example is given in fig.7a, which features a three-particle massless phase-space integral.In general, diagrams in this category are always massless phase-space integrals, because the massive line is not cut by definition.Such massless phase-space integrals always come with a step function θ(s), which indicates that it arises from the discontinuity around a logarithmic branch cut starting at the origin s = 0 (or ω = 0).This situation is reminiscent of the unphysical branch cut corresponding to massless u-channel cut diagrams, which can be removed by the shifting procedure from section 3.4.Indeed, we find that the shifting procedure is sufficient to deal with these unphysical s-channel cuts, which is supported by the example to be treated in section 4.3.1.
Massless u-channel cuts Cuts in the u-channel are for our case unphysical by definition, as they do not occur in the cut-diagrammatic expansion of the Drell-Yan cross section.The simplest class of unphysical u-channel cuts are the ones where only massless lines are cut.
Examples of such cuts are depicted in fig.6 and fig.7b, in the case of one-and two-loop crossed-box diagrams, respectively.These cuts correspond to branch cuts of forward diagrams with the branch point at the origin, so they can be removed by the shifting procedure from section 3.4.In section 4.3.2 a two-loop example is discussed where this procedure is applied.
Massive u-channel cuts This is a new type of u-channel cut which first appears at the two-loop level.The presence of the massive line has the effect of shifting the branch point to ω = −1, as compared to situation of the massless u-channel cuts.In this case the shifting procedure cannot be applied, so new method must be introduced to remove this type of unphysical cut.In the next section we focus on this problem and propose a solution in the form of an extra prescription.A non-trivial test of that procedure is then presented in the context of the two-loop crossed box in section 4.3.2.
The various types of unphysical cuts correspond to branch cut discontinuities of forward amplitude diagrams, where the branch cut does not extend from ω = 1 to ω = ∞.The connection between the above-mentioned cut diagrams and branch cut discontinuities is summarized in fig.8.

Extracting series coefficients for two-loop forward amplitude diagrams
The analysis of the types of unphysical cuts of two-loop diagrams in the previous subsection calls for an extension of our method for extracting the series coefficients of forward amplitudes, as described in section 3 for one-loop diagrams.In particular, the massive u-channel cut requires a new prescription besides the shifting procedure in section 3.4.The shifting procedure itself works also at higher loops, but the series ansatz for a forward diagram in eq.(3.38) is particular to the one-loop case and needs to be generalised.Furthermore, when dealing with higher-loop diagrams, the discussion in section 3.5 on how to obtain the series coefficients from IBP's should be combined with the notion of master integrals.We start by discussing the latter.
In many calculations of scattering amplitudes at higher orders in the QCD coupling constant, the use of master integrals has proven to be extremely useful.There can be many diagrams contributing to a cross-section, which produce even more Feynman integrals upon working out tensor reduction.Typically, all those integrals can be written as special cases of a handful of topologies: integrals with as many linearly independent propagators as Lorentz invariants formed out of at least one loop momentum, raised to arbitrary powers.It has been shown that all integrals in a topology can be written in terms of a finite set of master integrals [39].The computation of Feynman integrals for cross-sections thus boils down to computing master integrals.For this reason we focus in the remainder of this section on applying our method to master integrals.
Let M denote a vector of n such master integrals which depend on ω and .Assume that the first k master integrals can be computed by applying known analytical formulae for oneloop bubbles successively.We indicate these with a superscript B. Then the vector of master integrals is written as M = (M B 1 , . . ., M B k , M k+1 , . . ., M n ).Given the fact that the bubbletype integrals M B i are known exactly in , they can serve as inhomogeneous terms for the differential equations for the remaining unknown M i .This works as follows.Gathering the unknown integrals in the vector M = (M k+1 , . . ., M n ), taking its derivative with respect to ω and reducing the result to master integrals yields a system of first-order coupled differential equations d dω M = A • M. Notice here that the right-hand side generally depends on all master integrals.This situation can be avoided by decoupling the differential equations, at the expense of raising the order of the differential equations [40].As a result, the differential equation for a given M i will then be of order r, which takes some value 1 ≤ r ≤ k depending on the particular system, and has the form Here, the free index i is bound by k + 1 ≤ i ≤ n and the coefficients a m and b j,m are rational functions of ω and .We emphasise that the right-hand side is known exactly in , since it consists of the known bubble-type integrals and derivatives thereof.The number of unknown integrals k can be large, in practice.This means that the order of the differential equation r ∈ [1, k] could be equally large, making it challenging to solve.Moreover, the rational functions a m and b j,m also grow in size as r increases.Such situations may be ameliorated by decoupling differential equations to subsets of master integrals in M.After each iteration the solutions that can be found exact in may be used as inhomogeneous terms as well, thus lowering the order of the differential equations for the next integrals to be calculated.
The next task is to solve the differential equations in eq. ( 4.1) for i = k + 1, . . ., n.The approach we shall take in this paper is that of inserting an ansatz for the M i in terms of series expansions in ω.The one-loop ansatz in eq.(3.38) already displays the key feature, namely of decomposing the function into analytic and non-analytic pieces.The non-analyticity at ω = 0 is captured by powers ω − .In the case of two-loop diagrams the ansatz will take the following form where the series coefficients c n , e n depend on but not on ω.Substituting this expression for M i into the differential equation eq. ( 4.1) and equating equal powers of ω produces difference equations for the series coefficients.The difference equations have the general form of eq.(3.39).Non-zero coefficients of ω k on the right-hand side of the differential equation supply boundary conditions to the difference equations.This means that additional computations to ascertain such boundary conditions, e.g. using expansion-by-regions, are (typically) not necessary.In simple cases, when the order of the difference equation is relatively low, the series coefficients might be solved exactly in , involving ratios of Gamma-functions.Otherwise, the series coefficients can always be solved order-by-order in in terms of harmonic sums, using an ansatz of the form given in eq.(3.40).Once one has the series coefficients c n in hand, we have essentially computed the forward diagram.
Then we are in the position to start modifying the forward diagram in such a way that its discontinuity no longer has any contribution from unphysical cuts.The first step in this process is the shifting procedure.Conceptually, this procedure prescribes exactly what was discussed for one-loop diagrams: terms in the series are replaced according to c n ω n−k → c n+k ω n .After shifting, a two-loop forward master integral thus becomes As discussed in section 3.4, this procedure removes the unphysical branch cut discontinuity that arise from massless s-channel cuts and/or massless u-channel cuts.Indeed, the function M i is expanded as a series around ω = 0.At a technical level, this shifting procedure can be a bit more involved than the one-loop case.Namely, if the coefficients c n , e n were expressed order-by-order in in terms of harmonic numbers, then the series in eq. ( 4.3) features harmonic numbers evaluated at non-integer values: S (n + k ).These functions must be expanded in to match the form of the rest of the expression, which boils down to taking derivatives of harmonic numbers with respect to their argument.To this end one uses the known analytic continuation of harmonic numbers from the integers to the real line.In practice, we make use of the package HarmonicSums [31,38,[41][42][43][44][45] to expand the harmonic numbers evaluated at non-integer values.
So far we have discussed the shifting procedure and the use of IBP's to extract series coefficients in the context of two-loop diagrams.This is sufficient to deal with all types of unphysical cut of two-loop diagrams, except for massive u-channel cuts.The latter type of cut diagrams correspond to a branch cut of the forward along (−∞, −1] in the complex ω-plane.Since the branch point is not at the origin ω = 0, it is not removed by the shifting procedure.In order to remove discontinuities around such branch cuts, we extend our method further. Figure 9: The analytic structure of the example functions f (ω) and f (ω), given in eq.(4.4) and eq.(4.11), respectively.The function f (ω) represents a forward diagram, whose discontinuity contains contributions from both physical and un-physical cut diagrams.The second function, f (ω), is a modified version of the forward, such that only the physical branch cut is present.
We shall replace the forward diagram by a new function, whose branch cut along ω ∈ (−∞, −1] is removed while its branch cut discontinuity around the physical region ω ∈ [1, ∞) remains unchanged.Our technique for obtaining a function that satisfies these requirements is perhaps best explained with the help of an example.Consider the following product of logarithms, denoted f (ω) in view of the absence of a branch point at ω = 0, In the complex ω-plane this function f (ω) has two branch cuts, which are located along the disconnected intervals (−∞, −1] and [1, ∞).This situation is shown in fig.9(a).The discontinuity of f (ω) is simply the sum of the discontinuities around the individual branch cuts: In eq.(4.5), the first term on the right-hand side may be interpreted as a contribution coming from unphysical cut diagrams.Removing the unphysical cuts thus amounts to removing the branch cut along (−∞, −1] from the function f (ω), leaving a new function, f (ω), such that The corresponding analytic structure is displayed in fig.9b.The question is how to find f (ω).Note that in there is no unique answer to this question.Indeed, any constant (or smooth function, for that matter) may be added without changing the discontinuity.This ambiguity is lifted by imposing the constraint f (0) = 0, which reflects the physical property of scattering cross-sections that they vanish in the limit of zero centre-of-mass energy.This constraint, together with the analyticity of f (ω) around the origin, allows us to write down a series representation The coefficients c n can be obtained from the Cauchy integral formula, taking a small contour enclosing the origin.Inflating the contour such that it wraps around the branch cut, the contour integral becomes the integration of the discontinuity along the real line, analogous to the discussion in section 2.1.Subsequently changing variables to the reciprocal z = 1/ω leads to the following Mellin-transform integral This standard integral transform may be performed (for more complicated cases one may use the MT package [46]) and the result is In the analogy with perturbative computations, these coefficients correspond to the Mellin moments of the sum of cut diagrams obtained from the forward f (ω) by taking all possible physical cuts.Considering the aim of this paper, these moments therefore provide a satisfactory answer.
For completeness, we will also determine the full function f (ω).Obviously, in a small neighbourhood around the origin, f (ω) is given by the series eq.(4.7) with coefficients in eq.(4.9).Its analytic continuation to the complex ω-plane is given in terms of polylogarithms.This continuation may actually be constructed by first rewriting the series coefficients as linear combination of harmonic sums with multiple indices [37], which essentially projects c n onto a convenient basis of the function space: where . With this expression in hand, the sum in eq.(4.7) may be evaluated in closed form, using the fact that the series coefficients of harmonic polylogarithms are harmonic numbers [31], and one finds One can check explicitly that this expression has the correct branch cut discontinuity, as required by eq.(4.6).This completes the example.The same method for removing the unphysical branch cut can be applied to two-loop forward diagrams.Apart from branch cuts, one then also deals with poles, typically at ω = 1.One simple way to implement the removal of the wrong branch cut is by deriving replacement rules for the individual harmonic numbers, which appear in the result of the shifting procedure eq.(4.3).For a given harmonic number S (n), one first evaluates the corresponding sum n S (n) ω n in closed form.Based on similar analysis as in the previous example, one then constructs a function which has the unphysical branch cut removed and whose series coefficients define the replacement of S (n).For example, in the case of which is equivalent to the effective replacement rule Following these steps with all harmonic sums produces a 'dictionary' of replacement rules, which may then be applied to any diagram.We shall use such replacement rules in section 4.

Two-loop examples
This subsection provides examples that serve to illustrate two main lessons from our studies of two-loop diagrams, namely: how to remove unphysical cuts from a forward diagram, and how to compute the series coefficients of forward diagrams at higher loops from differential equations.The first example below will illustrate how to deal with massless s-channel cuts.We demonstrate that the shifting procedure not only removes massless u-channel cuts, but also removes any unphysical s-channel cut.The second example shows the power of the method by applying it to a rather difficult forward amplitude diagram, the two-loop crossed box.The latter admits massive u-channel cuts, which can be treated along the lines of section 4.2.

Two-loop self-energy diagram
In our first two-loop example we study a forward self-energy diagram, whose cutting equation is depicted in fig.10.As illustrated in the figure, the forward diagram admits two cuts: a two-and a three-particle cut.The two-particle cut is physical, but the three-particle cut is an unphysical s-channel cut which needs to be removed from the forward.In this subsection we show how to compute the physical cut from the forward diagram and point out the differences with a direct calculation the cut diagram.
Let us start by computing the forward diagram, before proceeding to remove the unphysical cut in order to extract the moments of the physical cut.The forward two-loop self-energy diagram is given by where C( ) = (iπ 2− r Γ ) −1 and G 1,1,1,1,0 is embedded in the integral topology with denominators D i given by the following expressions in terms of P = p + p, We proceed to compute G 1,1,1,1,0 by establishing an appropriate differential equation.To this end, notice that raising the power of the massive propagator may be achieved by differentiation with respect to the mass Q 2 , that is Using IBP reduction, the integral on the left-hand side can be expressed in terms of simpler integrals The first integral on the right-hand side is the self-energy diagram at hand (up to a prefactor), the integral on the left-hand side is its derivative, and the last two integrals on the righthand side are simpler bubble-type integrals.The latter can readily be computed exactly in , producing As we discussed in the previous subsection, we now insert a series ansatz for S into this differential equation, to turn it into a difference equation for the series coefficients.From the inhomogeneous terms in eq.(4.20) one can infer that the forward diagrams will have the structure S = f 1 (ω)+ω −2 f 2 (ω), where f 1 (ω) and f 2 (ω) are regular functions of ω close to the origin.We thus proceed to make the series ansatz5 Equating same powers of ω produces two recursions, complete with boundary conditions: The solutions to these elementary recursions are ratios of gamma functions, As a result, the forward self-energy diagram S is now known as a series expansion around the origin: These series can easily be recognised as representations of 2 F 1 -hypergeometric functions, but for our purposes the current form is actually more useful.Indeed, the aim of the remainder of this section is to extract the Mellin moments of the physical cut in fig. 10 from the forward amplitude diagram in eq.(4.27).
Extracting the Mellin moments of the physical cut from the forward is done in the following way.We construct a new function S, which has the same branch cut discontinuity as S around ω ∈ [1, ∞), but does not possess a branch cut starting at the origin ω = 0.In practice, we find such a function by means of the shifting procedure, as explained in section 3.4.Applied to the series in eq.(4.27) this produces where we dropped an ω-independent term, without affecting the discontinuity.The series coefficients of this new function S (in contrast to S) are well-defined.If we write S = ∞ n=1 c n ω n , then its series coefficients c n are equal to Based on our arguments presented in section 4.2 we claim that these series coefficients provide the Mellin moments of the physical cut on the right-hand side of the cutting equation in fig.10.The coefficients in eq.(4.29) therefore constitute the main result of this example.
Let us now verify our claim.To this end we shall compute the physical cut diagram explicitly.One way to proceed is by applying reverse unitarity [4] to the IBP reduction in eq.(4.17), in order to derive a differential equation for the cut diagram.Alternatively, one may actually perform the phase-space integration directly.In the latter approach one simply integrates a massless sub-bubble over a two-particle (one-mass) phase space.The massless sub-bubble reads Because the massive line is cut, k 2 = Q 2 , the bubble can be pulled out of the phase-space integral.As a result, the cut diagram is given by The Mellin moments can be computed exactly in due to the simple dependence on z.They are given by This relation holds at all orders in , as claimed.

Two-loop crossed-box diagram
We turn to our second example: the two-loop crossed-box diagram, depicted in fig.11.This diagram is distinguished from previous examples in two key aspects.First, it is sufficiently complicated so that it cannot be calculated exactly in , thereby providing a testing ground for the techniques of the previous subsection for computing the series coefficients of forward diagrams order-by-order in from differential equations.Second, the diagram is the first example to admit a massive u-channel cut, for which a new procedure was developed also in the previous subsection.In the example below we focus on these two aspects: first we compute the forward diagram, after which the moments of the physical cut are recovered by means of the shifting procedure and the replacement rules.We finally cross-check our results against the literature.Our first task is to compute the crossed-box diagram in fig.11.It may be written as where G 1,1,1,1,1,1,1 is one of the integrals in the following two-loop double-box topology where the denominators D i are given by The integral of interest is the last (and most complicated) master integral M 14 .Following the notation introduced in section 4.2, the first eight integrals are marked with the superscript "B" to indicate that they can be readily computed as iterated bubble integrals.These integrals are

.38)
Being exact in , these expressions are allowed to appear as inhomogeneous terms in differential equations for the six remaining unknown master integrals.We proceed to derive decoupled differential equations for the master integrals M 9 through M 14 of the form in eq.(4.1), using the Laporta reduction algorithm in FIRE [49,50].Inserting the series ansatz eq.(4.2) for the two-loop forward master integrals, the differential equations then transform into difference equations.It turns out that three of those equations can be solved exact in , producing series coefficients expressed as ratios of Gamma functions.As a result, the ansätze are easily recognised as series representations of hypergeometric functions: These expressions are useful because the corresponding integrals may appear as inhomogeneous terms in differential equations for the remaining unknown integrals: M 10 , M 11 and M 14 .
In the remainder we focus on the computation of M 14 , which is the forward crossed-box diagram.Inspecting the first-order differential equation for this integral reveals that it is coupled to all other master integrals, in particular to the unknown integrals M 10 and M 11 .After decoupling those two, as described in section 4.2, we obtain a third-order differential equation for M 14 .As before we insert the series ansatz eq.(4.2) for the forward integral, which produces an eighth-order difference equation.The latter can be solved order-by-order in in terms of harmonic numbers, c.f. eq.(3.40), using the strategy outlined in section 3.5.If the full diagram then its series coefficients are found to be where S ≡ S (n − 1).We have checked the validity of the representation eq.(4.40) for the forward crossed-box diagram, by reconstructing from the infinite sums the full diagram in terms of harmonic polylogarithms and performing a numerical cross-check using SecDec [51][52][53].These series coefficients now form the starting point for the next phase, which is to extract the Mellin moments of the corresponding physical s-channel cut diagram.
Inspecting the analytical structure of B 2 from its representation in terms of harmonic polylogarithms reveals three branch cuts.They are located along the real axis in the domains ω ∈ (−∞, 0], ω ∈ (−∞, −1] and ω ∈ [1, ∞), which correspond to massless u-channel cuts, massive u-channel cuts and the massive s-channel cut, respectively.The first and second types of branch cuts in the forward diagram are unphysical; they will be removed by performing the shifting procedure and applying replacement rules, respectively.Let us start with the shifting procedure.As we have seen in previous examples, this amount to the transformation More explicitly, the newly defined coefficients c n are given by where again S ≡ S (n − 1).In order to arrive at this form for the series coefficients c n , we have made use of the package HarmonicSums to expand the harmonic numbers S (n + k ), which appear in the coefficients d n+ and e n+2 as a result of shifting, as a Taylor series in .From the formula in eq. ( 4.42) it is clear that B 2 is regular at the origin, so we have successfully removed the branch cut along ω ∈ (−∞, 0] from the forward diagram.Crucially, the discontinuities around the remaining two branch cuts are unchanged.This can be verified by explicitly computing and comparing the discontinuity of both B 2 and B 2 using the HPL package [54, 55] 6 .In terms of cutting equations, this elimination of unphysical branch cut in the forward diagram is to be interpreted as the elimination of cut diagrams on the right-hand side of the cutting equation, as indicated by the first two lines in fig.12. In the next stage we modify the forward diagram even further, in such a way that the second unphysical branch cut is removed as well.At the level of individual harmonic polylogarithms this task is performed simply along the lines of the example in section 4.2.The results translate to replacement rules for the harmonic numbers.In particular, harmonic numbers with only positive indices do not need to be altered: the corresponding "resummed functions" do not contain unphysical branch cuts.The first two orders in of c n therefore do not need to be modified.For the remaining harmonic numbers we apply the following replacement rules.We recall that these rules are derived in a diagram-independent way.At order −2 we need a single replacement rule: These coefficients c n constitute the main result of this example section.We have checked explicitly that resumming these coefficients yields an expression The validity of the above statement can be verified by comparing the coefficients c n against an explicit computation of the Mellin moments of the physical cut of the forward diagram B 2 , depicted in fig.13.An explicit result for this particular cut diagram was given in eq.(B.21) of ref. [4].Correcting for small misprints (see appendix A of ref. [56]) and adopting our normalisation convention, we have that

Conclusions
In this paper we have presented a method for computing Mellin moments of single-particle inclusive cross sections, such as Drell-Yan and Higgs production, directly from forward scattering diagrams by invoking unitarity in the form of cutting equations.Due to the non-inclusive nature of these processes, the cutting equations contain unphysical cuts.The main achievement of this paper is to a provide diagram-independent prescription for "removing" such unphysical cut contributions to the discontinuity of a forward diagram, once these are expressed in the reciprocal ω = 1/z variable.The removal occurs through a complex shift in the summation index, and through a replacement rule dictionary for harmonic sums in the results.After this, the modified sum over powers of ω reproduces precisely the desired sum of physical cuts, and the coefficients are precisely the Mellin moments of the corresponding contribution to the cross section.We have demonstrated our method for various one-and two-loop diagrams.
The approach of this paper is conceptually similar to the computation of three-loop DIS splitting functions [6,7].While DIS is a fully-inclusive process, our method provides a nontrivial extension to semi-inclusive processes.Other methods exist for obtaining cross sections of semi-inclusive processes, but they do not make use of the optical theorem or cutting equations.For example, one very successful approach [4] computes cut diagrams as solutions to differential equations.Technically the latter need to be augmented with boundary conditions coming from a separate calculation (typically expansion-by-regions).In our approach the boundary conditions to difference equations for the Mellin moments are provided by the results for bubble-type loop integrals.In these other approaches calculations are moreover performed in z-space, except in [9].
Our method thus provides a new means of computing semi-inclusive cross sections, at least up to two-loop order.Since the main ingredients to the method are forward loop diagrams, as opposed to cut diagrams, it is particularly useful as an alternative to corrections involving real radiation, but provides no alternative way to compute virtual corrections.Being exclusively made out virtual diagrams, numerical cross-checks may be performed in a uniform way for all contributions (see the two-loop examples in the previous section).
In regards to extending our method beyond two-loop order, we note that the work in this paper is based on an analysis of (un)physical branch cuts and the assumption that the solution space is spanned by harmonic sums.Both aspects will need to be reviewed at higher loops, but we are hopeful that progress can be made towards single-scale cross sections at N 3 LO.

Figure 2 :
Figure 2: Generic cuts of forward aplitudes with two initial massless particles with momenta p and p.The cuts of type (a) and (c) vanish by the general cutting rules.Note that the u-channel cut differs from the t-channel cut by interchange of the two outgoing momenta.

Figure 3 :
Figure 3: Diagrams needed for the one-loop DY cross sections.Diagrams obtained by complex conjugation or exchanging p ↔ p are omitted.Arrows on the lines indicate momentum flow.

Figure 4 :
Figure 4: Cutting equation for the one-loop triangle diagram.The right-hand side features both Cut 1 T , a physical massive s-channel cut, and Cut 2 T , an unphysical massless s-channel cut.

For the box diagram B 1 Figure 5 :
Figure 5: Cutting equation for the one-loop box diagram.

Figure 6 :
Figure 6: Cutting equation for the one-loop crossed-box diagram B 2 , featuring on the right-hand side both the physical s-channel cut (Cut 1 B 2 ) and the unphysical u-channel cut (Cut 2 B 2 ).

Figure 7 :
Figure 7: Types of unphysical cuts appearing at two loops.Example are: (a) massless s-channel cut; (b) massless u-channel cut; (c) massive u-channel cut.

Figure 8 :
Figure8: Branch cut structures that appear in two-loop forward diagrams, listed together with cuts that describe the corresponding branch cut discontinuities according to the cutting equation.Only the first type of cut is physical, i.e. contributes to the Drell-Yan cross section.

Figure 10 :
Figure 10: Cutting equation for a two-loop self-energy diagram.

Figure 13 :
Figure 13: The only physical cut of the two-loop forward crossed-box diagram.
only has a branch cut along ω ∈ [1, ∞) and whose discontinuity along that branch cut is the same as for the original diagram B 2 .This means that these series coefficients c n in eq.(4.47) must be equal to the Mellin moments of the sum of physical cuts of the forward diagramB 2 !We claim that 1 2πi M n [Cut phys B 2 ] = c n .(4.49)