Exponentiating virtual imaginary contributions in a parton shower

The operator in a parton shower algorithm that represents the imaginary part of virtual Feynman graphs has a non-trivial color structure and is large because it is proportional to a factor of $4\pi$. In order to improve the treatment of color in a parton shower, it may help to exponentiate this phase operator. We show that it is possible to do so by exponentiating matrices that are no larger than $14\times14$. Using the example of the probability to have a gap in the rapidity interval between two high transverse momentum jets, we test this exponentiation algorithm by comparing to the result of treating the phase operator perturbatively. We find that the exponentiation works, but that the net effect of the exponentiated phase operator is quite small for this problem, so that one can as well use the perturbative approach.


I. INTRODUCTION
There have been efforts to improve the accuracy of leading order parton shower programs by moving away from the leading color (LC) approximation in the treatment of QCD color [1][2][3][4][5][6][7][8][9]. For this effort to be consistent with quantum mechanics, one needs to work with bra and ket color amplitudes and to implement the color content of virtual graphs acting on these amplitudes [1,3,4,[7][8][9].
In the paper [8], we explored how to include QCD color in a parton shower in a systematically improvable way. Then in Ref. [9] we applied this algorithm, as implemented in our program Deductor, to the problem of gaps [10][11][12][13][14][15][16][17][18][19][20][21][22][23] in the rapidity interval between two high transverse momentum jets. In this paper, we explore the possibility of exponentiating the contributions from the imaginary part of virtual graphs, which are proportional to iπ times a non-trivial color operator.
The basic idea in Ref. [8] is simple. In an ideal, but not practical, leading order parton shower with exact color, the shower is represented as an operator U(t, t 0 ). This operator takes the shower from an initial stage at a hard momentum scale, corresponding to a shower time t 0 , to a later stage with a softer momentum scale, corresponding to a shower time t [24,25]. The operator U(t, t 0 ) acts on a linear space of states ρ , the statistical space. Ignoring spin but accounting for color, this space consists of functions of partonic momentum and flavor variables together with a density matrix in the space of the colors of the partons. The partonic momentum and flavor variables, Here H I (t) is an operator that generates parton splittings (at first order in α s ), increasing the number of partons by one each time it acts, and V(t) is an operator that approximates virtual graphs and leaves the number of partons unchanged. The operators H I (t) and V(t) have complicated color structures, so that it is not practical to implement Eq. (2) directly. There is an approximation, the LC+ approximation, that generalizes the leading color approximation [1]. In the LC+ approximation, we replace H I (t) and V(t) by approximate operators H lc+ (t) and V lc+ (t). Then we solve The solution takes the form U lc+ (t, t 0 ) = N lc+ (t, t 0 ) Here N lc+ (t 2 , t 1 ) is the no-splitting operator, The operator V lc+ (τ ) is diagonal in the color basis that we use, the trace basis, so that it is practical to calculate its exponential. Then the diagonal matrix elements of N lc+ (t 2 , t 1 ) constitute the Sudakov factor that is needed to choose the next splitting time in the parton shower algorithm. Now, the LC+ approximation leaves out operators ∆H(t) and ∆V(t), The operator ∆V(t) contains two terms: The first term, ∆V Re (t), approximates the terms in the real part of virtual graphs that are not included in V lc+ (t). The second term, V iπ (t), corresponds to the imaginary part of virtual graphs. In order to set up the general expansion in powers of ∆H(t), ∆V Re (t), and V iπ (t), we first define an operator X (t, t 0 ) by This equation is to be solved iteratively to give a series in powers of ∆V: Then the evolution equation for U(t, t 0 ) becomes This is also to be solved iteratively. The iterative expansions of Eqs. (9) and (10) give U(t, t 0 ) as series in powers of ∆H(t), ∆V Re (t), and V iπ (t).
Then we can specify non-negative integers N Re and N iπ and retain contributions proportional to In the present paper, we ask whether one could exponentiate V iπ in Eq. (7), so that we effectively set N iπ = ∞. If one could do this, one would have a parton shower that is well adapted for observables in which V iπ is effectively a large operator.
To include all powers of V iπ , we would like to replace N lc+ (t 2 , t 1 ) in Eq. (5) by Here T indicates time ordering, with the latest shower times to the left. Then we would use not the full ∆V(t) in Eq. (8), but only ∆V Re (t): With V iπ exponentiated, the evolution equation (10) becomes

II. EXPONENTIATING Viπ
In order to examine whether it is practically possible to include V iπ along with V lc+ in the no-splitting operator, as in Eq. (11), we need to understand the color structure of V iπ . The structure of V iπ (assuming massless partons) is given in Eq. (10.14) of Ref. [1], Here T a represents the insertion of a color matrix T c on incoming parton line "a", T b represents the insertion of a color matrix T c on incoming parton line "b", and the dot in (T a · T b ) represents a summation over the octet color index c. In [(T a · T b ) ⊗ 1] the color matrices act on the ket state, while in [1 ⊗ (T a · T b )], they act on the bra state. When we take the color trace, we get 1 V iπ (t) = 0 because the two terms in Eq. (14) have opposite signs. At first sight, the goal of exponentiating V iπ in a practical algorithm seems unreachable because the color operators in Eq. (14) act on a color space that has a very large dimensionality. In fact, an n gluon statistical state vector is expanded in [(n − 1)!/2] 2 basis vectors if we use the trace basis. For n = 10 this is 3 × 10 10 basis vectors.
However, let us look at this more closely. When the operator N e (t 2 , t 1 ) operates on a statistical space basis vector {p, f, c ′ , c} m , it leaves the number of partons, their momenta, and their flavors unchanged. It acts as a matrix on the colors. However, this matrix is simpler than it would seem because it is a product of two matrices, one for the ket colors and one for the bra colors: In order to simplify the notation, we can write these matrices as operators, an operator n on the ket color space, (16) and an operator n † on the bra color space, The operator n has the form Here T indicates ordering according to the shower time t. In n(t 2 , t 1 , {p, f } m ), which operates on the ket color space, the latest shower times are on the left. In n † (t 2 , t 1 , {p, f } m ), which operates on the bra color space, the latest shower times are on the right and the sign of the imaginary part of the exponent is reversed. The operator a in the exponent of Eq. (18) has the decomposition [8]: First, there is a lc+ coll . This operator comes from "direct" graphs, in which a parton is emitted from parton line l and absorbed on the same parton line, There is a sum over the index of the emitting parton l. We integrate over momentum splitting variables Λ, z, φ and flavor splitting variables. Here Λ is a hardness variable, z is a momentum fraction, and φ is an azimuthal angle. The splitting variables are collectively denoted by ζ. The function T l (ζ, {p} m ) defines the corresponding shower time t. There is a function λ ll ({p, f } m , ζ) that contains a factor α s and depends on the splitting variables and the momenta and flavors of the partons before the splitting [8]. This function contains a color factor C F , C A , or T R that is determined by the flavor of parton l. However, this factor is independent of the color state {c} m on which it acts. That is, v lc+ coll is proportional to the unit operator I c on the ket color space.
The operator a lc+ soft corresponds to graphs in which a gluon is exchanged between two partons, l and k. The operator a lc+ soft (t, {p, f } m ) has the form, when acting on a color basis vector, There is a function λ lk ({p, f } m , ζ) that contains a factor α s and depends on the splitting variables and the momenta and flavors of the partons before the splitting. Finally, the color information is provided by a function χ(k, l, {c} m ), which is 1 if partons l and k are next to each other along a color string in the basis state {c} m and is 0 otherwise [8]: Notice that the color basis states {c} m are eigenvectors of a lc+ soft (t, {p, f } m ). The color basis states {c} m are not eigenvectors of T a ·T b , so it would seem difficult to include the imaginary term proportional to T a · T b in the exponent in Eq. (18). To see what is involved, let us start by considering the case that both incoming partons are antiquarks. Then incoming parton "a" carries an outgoing quark color index that we can call α and incoming parton "b" carries an outgoing quark color index that we can call β. The outgoing partons carry color indices that we can collectively denote by {r}. This gives us a (possible unnormalized) color basis state G with color wave function G(α, β, {r}). With this notation, we can define what T a · T b does to G : where, using the definition that T a inserts a color generator on line "a" and T b inserts a color generator on line "a," There are implied sums over repeated color indices. We will want to consider two basis states G(n) , n ∈ {1, 2}. We take G(1) to be the color basis state with which we start: G(1) = {c} m . Then we define a second basis state G(2) . For the color wave functions G(n, α, β, {r}), we choose where R is a fixed color state vector that carries indices α ′ , β ′ for two outgoing quark lines together with the indices {r} for the rest of the partons. The two color wave functions C(n; α, β) α ′ β ′ are Then after applying T a · T b , we get new basis states where We can evaluate this. We can use the Fierz identity, This gives us This calculation gives us where For our complete basis states G(n) , this is That is, for two antiquark incoming states, the operator T a · T b acts as a matrix on just a two dimensional subspace of the color ket space or the color bra space. We do not need to deal with a color space of very high dimensionality. Now we can return to the operator n in Eq. (18). The time ordered exponential in Eq. (18) can be implemented as an ordinary exponential of a matrix by using the Magnus expansion [26]. We have The first term in the exponent is simply [α s (τ )/(2π)] a(τ, {p, f } m ) without time ordering: The second term in the exponent is The higher order contributions can be obtained by recursion.
The commutator in the second order term simplifies when we use the decomposition (19) of a. The operator a lc+ coll (t, {p, f } m ) contains the soft×collinear singularities of a. Thus at fixed t, this operator contains a factor of t, which comes from integrating over the momentum fraction in a splitting. However, we saw in Eq. (20) that a lc+ coll (t, {p, f } m ) is proportional to the unit matrix, so that it commutes with the other terms in a. We also saw in Eq. (21) that the soft part of a is a diagonal operator. However, T a · T b does not commute with a lc+ soft . Thus the commutator is The operator a lc+ soft (t, {p, f } m ) contains a soft singularity only, so it does not have a contribution proportional to t. From its functional form [8], one sees that this operator is independent of the shower time except for shower times close to 0. It has the expansion With this, the commutator in Eq. (37) is We see that the second term of the exponent, ω 2 , in Eq. (34) is not only second order in α s but is also suppressed by powers of e −t so that it does not contribute in the soft and collinear enhanced regions. The higher order terms in Eq. (34) are similarly suppressed. Thus these contributions can reasonably be neglected and we simply write n in Eq. (18) as the ordinary exponential of a.
We have seen that the operator n, including its iπ terms, can be written as an exponential of a finite dimensional matrix in the case of an incomingq,q state. In the following subsections, we examine the other possible choices for the flavors of the initial state partons. In each case, we obtain a finite dimensional matrix that needs to be exponentiated numerically. The largest dimensionality is 14 × 14.
A. Incoming q q For an incoming q q state, the incoming quarks carry outgoing antiquark indices. We then consider basis states of form with When we exchange a gluon with octet color index c between the two outgoing antiquark lines, we have Here the SU (3) generator for the first outgoing antiquark line is −t(c) α ′′ α and the generator for the second antiquark line is −t(c) β ′′ β . The two −1 factors cancel. After using the Fierz identity, we obtain with the same matrix as in theq,q case: B. Incoming qq Here we have one outgoing quark index α and one outgoing antiquark index β. We then consider basis states of the form When we exchange a gluon with octet color index c between the two external lines, we have Here the SU(3) generator for the outgoing quark line is t(c) α α ′′ and the generator for the antiquark line is −t(c) β ′′ β . After using the Fierz identity (29), we obtain with the a new matrix: C. Incomingq g Here we have one outgoing quark index α and one outgoing gluon index a. We then consider basis states G(n) , n ∈ {1, . . . , 4}, with G(1) = {c} n and with color wave functions of the form The matrix R representing the rest of the color vector is given a quark index β ′ and an antiquark index β instead of a gluon index because in the trace basis the external gluon couples to a quark line. We define We use four basis states because there are four ways to attach the indices. Suppose now that we exchange a gluon with octet color index c between the outgoing quark line and the outgoing gluon line, giving us Let us start with C(1; a, α) β α ′ β ′ . We use the commutator rule, (53) Then we use the first term of the Fierz identity (29), noting that the second term of the Fierz identity gives canceling contributions. This gives This is The analogous calculation for C(2; a, α) β α ′ β ′ gives This is Now we can turn to C(3; a, α) β α ′ β ′ : This is The analogous calculation for C(4; a, α) β α ′ β ′ gives This is This calculation gives us where In this case, we have one outgoing antiquark index α and one outgoing gluon index a. We then consider basis states of form We define Suppose now that we exchange a gluon with octet color index c between the outgoing antiquark line and the outgoing gluon line, giving us (66) There is a minus sign here because we attach the color matrix to an antiquark line.
Then a calculation using the identities (53) and (29) gives us where E. Incoming g g In this case, we have two outgoing gluon indices a and b. We consider basis states of form Each vector C(n; a, b) carries two (outgoing) quark indices α and β, and two (outgoing) antiquark indices α ′ and β ′ . There are 14 ways to couple the two gluon indices to the quark and antiquark indices: Suppose now that we exchange a gluon with octet color index c between the outgoing gluon lines, giving us Then a calculation similar to the calculations that we have done previously gives where (73)

III. ASSEMBLING THE SHOWER
It is now clear what to do to complete the shower evolution equations. Starting with states G(1) = {c} m and G ′ (1) = {c ′ } m , we define color states G(n) and G ′ (n) . This gives us an operator n(t 2 , t 1 , {p, f } m ) defined as a matrix that mixes the states G(n) and an operator n † (t 2 , t 1 , {p, f } m ) that mixes the states G ′ (n) according to Eqs. (18), (19), and (33). Then the combined no-splitting operator N e (t 2 , t 1 ) is defined by the direct product of these matrices, according to Eqs. (15), (16), and (17). Now, in a complete algorithm, one would calculate X e (t 2 , t 2 ) as a power series in ∆V Re (τ ) by using Eq. (12). Then one would produce the complete shower using Eq. (13). One would keep any number of splittings generated by H lc+ (τ ) and ∆H(τ ). This would generate contributions proportional to powers of ∆H(τ ) and ∆V Re (τ ). One would retain contributions proportional to [∆H] A [∆V Re ] B with A + B ≤ N Re for a chosen value of N Re .
For this paper, we have not implemented the complete algorithm. Rather, we have set N Re = 0. This allows us to test whether the exponentiation of the iπ contributions as outlined above is computationally practical. It also allows us to test whether exponentiation of the iπ contributions makes a substantial numerical difference compared to simply keeping just a small number of factors of V iπ . For this purpose, we will use the gap fraction, for which we have found previously that the inclusion of color beyond the LC+ approximation can be numerically important [9].

IV. RESULTS
In this section, we explore the exponentiation of V iπ by examining its effect on the gap survival probability for jet production in proton-proton collisions at √ s = 13 TeV. Using the anti-k T jet algorithm [27] with a radius parameter R = 0.4, we find jets with transverse momenta P T and rapidities y in the region −Y cut < y < Y cut . We will use Y cut = 4.4. Label the two highest P T jets 1 and 2, with y 1 < y 2 . Definē p T = 1 2 (P T,1 + P T,2 ) , Now, define a cut parameter p cut T . We will take p cut T = 20 GeV. Look at those jets with P T > p cut T in the rapidity region y 1 < y < y 2 between the two leading jets. We will say that the event has a rapidity gap if there are no such jets in this rapidity region. We denote the fraction of events with a rapidity gap by f (p T , y 12 ).
As discussed in the previous section, we have not implemented the inclusion of the operators ∆H(τ ) and ∆V Re (τ ) in the shower along with the exponentiation of V iπ (τ ). For that reason, we do not include ∆H(τ ) and ∆V Re (τ ) in the shower operator U(t, t 0 ) in any of the calculations in this section.
Just after the hard scattering, Deductor inserts an operator U V that produces a summation of threshold logarithms [28]. The operator U V gives results as an expansion in powers of ∆V Re . Deductor retains only those terms with no more than N thr ∆ factors of ∆V Re . We choose N thr ∆ = 0 or N thr ∆ = 1, as we will specify, in the investigations in this section.
We begin by examining results for f (p T , y 12 ) with several values of N iπ and with the V iπ contributions exponentiated. Our aim is to see if the perturbative results with increasing values of N iπ are consistent with the exponentiated result. In this investigation, we examine f in the range 300 GeV <p T < 400 GeV and 4 < y 12 < 5. We let the parton shower run over the evolution range µ s > Λ > 30 GeV and stop evolution at Λ = 30 GeV.
Later in this section, we will choose the maximum color suppression index I max = 4. However, this choice allows a systematic error of about 1/N 4 c ≈ 0.01 and we want to avoid that here, so we choose I max = 6, even though this gives us larger statistical errors.
In Fig. 1, we exhibit the results for f with N iπ =  Fig. 1, we choose N thr ∆ = 0. We see that there is a jump of 0.014 between N iπ = 0 and N iπ = 2. For larger values of N iπ , the results appear to oscillate, but the statistical errors increase. All that we can really say is that the result for large N iπ seems to be between 0.200 and 0.215. This is consistent with the exponentiated result, f = 0.204 ± 0.04. That is, the perturbative results and the exponentiated result agree that the change from the N iπ = 0 result is ∆f ≈ 0.003±0.004. It is remarkable that the effect of V iπ is so small. If the effect were of order ∆f ≈ 0.1, it would be easy to see, but an effect of order ∆f ≈ 0.01 is difficult to see in the face of systematic and statistical errors of the program.
In Ref. [9], we examined the effect of color on the rapidity gap survival fraction for a range of y 12 andp T . We found that the effect of having up to two powers of V iπ was always quite small. It is of interest to see if the effect of including V iπ in exponentiated form is also small for all y 12 andp T .
For this study, results were obtained with the stated values N iπ or with exponentiated V iπ in the shower between Λ = µ s and Λ = Λ min = 30 GeV. For the rest of the shower, down to Λ = 1 GeV, we used the LC+ approximation. Throughout, we set the maximum color suppression index to I max = 4 and we set N thr ∆ = 1. We choose five different bins for y 12 and examine f (p T ) as a function ofp T . The results are shown in Fig. 2. In each y 12 bin, we show three curves. The first, in blue, is obtained with N iπ = 0. Then, in green, we show results obtained with N iπ = 2. Finally, in red, we show the result with V iπ in exponentiated form, N iπ = ∞. The result with up to two powers of V iπ (N iπ = 2) is generally somewhat larger than that with exponentiated V iπ , but the difference is typically not larger than 0.02.
It is remarkable that f (p T , y 12 ) with V iπ in exponentiated form is always quite close to its value with N iπ = 0. This is consistent with what we found in Ref. [9] and in Fig. 2: allowing two powers of V iπ perturbatively makes a rather small difference in the gap fraction. This does not mean that color beyond the LC+ approximation is unimportant for the gap fraction. We found in Ref. [9] that allowing up to two powers of ∆H and ∆V Re can make a substantial difference when log(p T /p cut T ) and y 12 are large.

V. CONCLUSIONS
The accuracy of leading order parton shower programs could be enhanced by improving the treatment of QCD color [2,[4][5][6][7]. We have worked to make our program, Deductor, more accurate with respect to color [1,3,8,9]. In the organization of Deductor, there are three operators that have a color structure that cannot be incorporated into the LC+ approximation that Deductor uses: ∆H, ∆V Re , and V iπ .
We found in Ref. [8] that it is possible to include the operators ∆H, ∆V Re , and V iπ perturbatively. That is, one can work with the result expanded in powers of these operators and retain the contributions up to some number of powers. There are practical limitations on the number of powers that one can retain, but two powers turns out to be quite practical.
The operator V iπ , which corresponds to the imaginary part of one-loop virtual graphs, is of special interest because of the role that it plays in analytical summation of large logarithms and because it contains a factor of 4π, which makes it a numerically large operator. One wonders if one could include V iπ to all orders. That is, one wonders if one could use V iπ as part of the exponentiated Sudakov factor that occurs between real emissions.
In this paper, we have found that one can include V iπ in the Sudakov exponent and we show how to do that. One has to calculate numerically the exponential of matrices, but the dimensionality of the matrices is at most 14 × 14.
In principle, one could include V iπ in the Sudakov exponent while including powers of ∆H and ∆V Re perturbatively. This would, however, require substantial effort to realize in computer code, so we have not attempted this further step.
We have tested this summation numerically, using the rapidity gap survival probability, which we have found is generally sensitive to color [9]. For the numerical test, we compare the results from the perturbative approach to the result with V iπ included in the Sudakov exponent. The result of the two calculations agree that the effect of V iπ is small, near the limit of what we can measure numerically. We have also examined the effect of V iπ , both in exponentiated form and perturbatively, for a wide range of the variablesp T and y 12 that appear in the rapidity gap problem. We find that the effect of V iπ is always quite small.
It is perhaps not a surprise that the net V iπ contribution is quite small. We note that if we integrate Eq. (14) for V iπ over a range of scales from p cut T = 20 GeV tō p T = 300 GeV and leave out the color factors, we have a net phase φ ∼ 2α s log((300/20) 2 ). With α s ∼ 0.1, this is φ ∼ 1. Thus V iπ is not a small operator in the context of the rapidity gap problem. However the contributions from V iπ have opposite signs when acting on ket color states compared to the sign when acting on bra color states. We note that for real numbers a and φ, exp(a + iφ) exp(a − iφ) is just exp(2a), with no contribution from the phase φ. Thus the effect of V iπ has to arise because the relevant color matrices don't commute. In our numerical tests of the perturbative approach to color in Ref. [8], we found that the effects of V iπ separately on the bra and ket states were of order 1, but that these effects cancelled, as they should, when we examined the unitarity property of shower evolution. Only when we have a color sensitive observable can the contributions fail to cancel. Apparently, the gap fraction is not sufficiently color sensitive to prevent a high degree of cancellation.
Because the effect of V iπ is small for the gap fraction observable, we believe that it suffices to examine observables that might be sensitive to color by treating V iπ perturbatively, along with ∆H and ∆V Re . If one finds an observable in which V iπ has a numerically substantial effect, then we could consider exponentiating this effect in Deductor. Of course, other workers may want to use the method of this paper to include V iπ as part of the Sudakov exponent in other parton shower programs.