Scattering and gluon emission in a color field: a light-front Hamiltonian approach

We develop a numerical method to nonperturbatively study scattering and gluon emission of a quark from a colored target using a light-front Hamiltonian approach. The target is described as a classical color field, as in the Color Glass Condensate effective theory. The Fock space of the scattering system is restricted to the $\ket{q}+\ket{qg}$ sectors, but the time evolution of this truncated system is solved exactly. This method allows us to study the interplay between coherence and multiple scattering in gluon emission. It could be applied both to studying subeikonal effects in high energy scattering and to understanding jet quenching in a hot plasma.


I. INTRODUCTION
The general picture of a high-energy dilute probe scattering off a color field is a commonly used approach for many different processes in QCD phenomenology. Scattering processes that probe the color glass condensate (CGC) [1] state of small-x gluons inside a high-energy hadron or nucleus are described in terms of infinitely energetic partons passing through an infinitesimally thin color field sheet, using the eikonal approximation. In order to study the phenomenon of jet quenching and radiative energy loss, one studies the situation when a high-energy parton passes through an extended colored medium and loses energy by gluon emission [2][3][4][5][6][7]. In both cases, one often performs analytical calculations in a kinematical approximation where the probe has an infinitely large energy. For realistic collider phenomenology in both physical situations, it is important, however, to be able to relax this approximation. For scattering off a CGC color field, subeikonal effects [8,9] can be important at realistic collider energies, such as at the upcoming Electron-Ion Collider [10]. This is the case, in particular, for the physics of spin at high energies [11][12][13][14][15]. Also for jet quenching, understanding the interplay between the coherence time of the emission and the timescales of the scattering centers of the medium is an area of active study [16][17][18][19][20]. Here, we address this problem using a nonperturbative approach. We consider the scattering of a highly energetic quark off a strong classical background field, and we treat the quark in a Fock space consisting of |q and |qg sectors. We explicitly solve the time evolution of this system with the light-front Hamiltonian formalism, using the time-dependent basis light-front quantization approach (tBLFQ) [21].
The tBLFQ approach is a nonperturbative computational method to investigate time-evolution problems. It is based on light-front quantum field theory and the Hamiltonian formalism. The implementation of the basis function representation allows one to choose a basis with the same symmetries of the system under investigation, and is therefore advanta- * mliy@jyu.fi † tuomas.v.v.lappi@jyu.fi ‡ xbzhao@impcas.ac.cn geous for carrying out efficient numerical calculations. This method has been previously applied to nonlinear Compton scattering [21,22], to the interaction of an electron with intense electromagnetic fields [23], and to quark-nucleus scattering [24].
In the earlier treatment of quark-nucleus scattering with tBLFQ presented in Ref. [24], the Fock space of the quark was truncated to the leading sector |q . In this limit, the subeikonal effect was revealed in the transverse coordinate distribution of the quark. In this work, we extend the Fock space to |q + |qg , thus including gluon emission and absorption in the process. We treat the target nucleus as a classical SU(3) color field given by the McLerran-Venugopalan (MV) model [25][26][27]. In the usual CGC treatment, the scattering only depends on the field integrated over the longitudinal coordinate. The method introduced here, however, can be applied to a more general situation where the process can be sensitive to the structure of the field in the longitudinal direction. We explicitly solve for the time evolution of the quark as a quantum state inside the target color field. The time dependence is sensitive to all three parts of the Hamiltonian of our system: the interaction with the background field, gluon emission and absorption, and phase rotation with the light-front energy of the state. The phase rotation is neglected in the eikonal limit usually used in CGC calculations, and it encodes the physics of the formation time of the radiated gluon. In our full nonperturbative treatment, we can smoothly vary the magnitudes of these three effects separately.
We study the evolution of the quark by looking into its distribution in phase space, including the longitudinal momentum, the transverse momentum, light-front helicity, and color. Our focus in this paper is on presenting and testing the numerical method, and demonstrating it in different physical regimes. For clarity, we use an initial condition of a pure |q state with a specific color and helicity so that the |qg components are generated only by the interactions. The only exception is when studying the sole effect from the interaction with the background field, where we also include a |qg component in the initial state. In the future, we aim to apply this numerical method to different physical situations, such as a high-energy scattering with subeikonal effects, which requires choosing initial conditions and measured observables corresponding to the physical process of interest. The layout of this paper is as The quark is moving along the positive z direction and it scatters on the nucleus which moves along the negative z direction. The dashed line is the worldline of the quark, z = β q t with β q the speed of the quark. The quark state is a superposition of the |q and the |qg states. The quark line is dressed by helical lines representing the gluon in the |qg state. The band represents the worldlines of the target nucleus, bounded by z = −β A t and z = −β A t + d . Here, β A is the speed of the nucleus and d = d 1 − β 2 A with d the width of the nucleus in its rest frame. In the ultrarelativistic limit of β A → 1, the red band in the diagram shrinks to a single line aligned with x + = 0.
follows. We first introduce the formalism of tBLFQ for the case of a quark emitting/absorbing a gluon and scattering on a color field in Sec. II. We then present and discuss numerical results in Sec. III, highlighting the effects of the three different parts of the Hamiltonian separately and together. We conclude the work in Sec. IV.

II. METHODOLOGY: TIME-DEPENDENT BASIS LIGHT-FRONT QUANTIZATION (TBLFQ)
The basic physical situation in our study is a high-energy quark moving in the positive z direction, scattering on a highenergy nucleus moving in the negative z direction, as shown in Fig. 1. The quark has momentum P µ with P + P − , P ⊥ whereas the nucleus has momentum P µ A with P − A P + A , P A,⊥ . We treat the quark state at the amplitude level and the nucleus as an external background field. The quark state is a superposition of the |q and the |qg states. The quark interacts with the nuclear field over a finite distance in light-front time 0 ≤ x + ≤ L η . The light-front quantization formalism is manifestly boost invariant in the z direction. Thus the same physical process can be described in different Lorentz frames with equivalent results. In practice, this means that the change of the P + momentum of the incoming quark can be compensated by a corresponding Lorentz contraction of the x + dependence of the target (both its size and internal structure). The physically genuinely different regimes correspond to different relative timescales of coherence and the background field interactions. For practical simulations, however, we choose specific numerical values, expressed here in GeV for concreteness.

A. The light-front Hamiltonian
The Lagrangian for the process we are considering is the QCD Lagrangian with an external field, where F µν a ≡ ∂ µ C ν a − ∂ ν C µ a − g f abc C µ b C ν c is the field strength tensor, D µ ≡ ∂ µ + igC µ the covariant derivative, and C µ = A µ + A µ is the sum of the quantum gauge field A µ and the background gluon field A µ .
The light-front Hamiltonian is derived from the Lagrangian through the standard Legendre transformation [28] in the light-cone gauge of the quark, i.e.,A + = A + = 0, and we show the detailed derivation in Appendix B 1. Here, we focus on the Hamiltonian in the truncated Fock space that we are actually working with.
The interacting quark state admits an infinite Fock space expansion in terms of the bare states. The dimensionality of this Fock space grows with the number of basis states (color, helicity, and momentum states) to the power of the number of particles. This growth makes it intractable when numerically going beyond higher orders in the Fock state expansion. Here, we truncate this expansion to the leading two sectors, |q and |qg , |q dressed = ψ q |q + ψ qg |qg + · · · , where ψ q (ψ qg ) is the probability amplitude of the |q (|qg ) sector, and "· · · " includes all the other Fock sectors with gluons and sea quarks, such as |qgg and |qqq , which are not considered in this work. In the truncated Fock space, the light-front Hamiltonian consists of two parts, P − (x + ) = P − KE + V(x + ), where P − KE is the kinetic energy and V(x + ) the interaction. Note that we do not consider the kinetic energy of the background field. The kinetic energy part of the Hamiltonian is a sum of single particle energies, The interaction part of the Hamiltonian consists of two terms, V(x + ) = V qg + V A (x + ), and its diagrammatic representation is illustrated in Table I. The first term V qg is the interaction between the quark and the dynamical gluon: It accounts for gluon emission and absorption inside the dressed quark state. The second term V A (x + ) includes the interaction of the background field with the quark and that with It admits an explicit time dependence introduced by the background field. Note that for an infinitesimal time step, the quark and the gluon interacting with the background field are separate interaction terms. The usual CGC picture of both the quark and gluon being rotated by the shockwave field of the target arises after iterating these interactions over several time steps.
The background field A µ accounts for the target, and we describe it using the MV model [25,26,29]. This is a classical field satisfying the reduced Yang-Mills equation, and it has only one nonzero component A − . In the MV model, one assumes that the target field is independent of x − (the light-front time for a left-moving target). This is justified by the probe's large momentum P + , which means that the x − dependence of the probe is larger than the x − scales of the target. The consequence of this approximation is that the longitudinal momentum P + of the probe is preserved in the interaction. The gluon mass m g is introduced to regularize the infrared (IR) divergence in the field, which simulates color neutrality on the source distribution [30]. We take m g = 0.1 GeV in the numerical simulations. The background field can be expressed in terms of Green's function as where The color charges are treated as Gaussian stochastic variables that are uncorrelated between different points in the transverse plane and between different points in light-front time. They satisfy the correlation relation Note that the parameterμ 2 has dimensions of GeV 3 , consisting of GeV 2 for the transverse dimension x ⊥ and GeV for the target's longitudinal dimension x + . This corresponds to the transport coefficientq in jet quenching [16]. For a highenergy scattering process, what matters is the charge density (g 2μ ) 2 integrated over the extension of the field along x + [31,32]. This integrated quantity, corresponding to the typical transverse momentum transferred by the target color field, is known as the saturation scale Q 2 s . For a field with constant charge density, it can be obtained, up to logarithmic corrections, from the product of (g 2μ ) 2 and the duration of the field L η [33]. The conventions regarding factors of π and 2 differ between different sources in the literature. Here, we use the fundamental representation saturation scale, which we take to be given by the relation neglecting the logarithmic corrections. Here, C F = (N 2 c − 1)/(2N c ) = 4/3 is the second-order Casimir invariant in the fundamental representation.

B. Time evolution of the state
The evolution of quantum states is governed by the timeevolution equation on the light front. Since we are interested in how the quark evolves under the interaction, it is natural to use the interaction picture (denoted by the subscript I), In the interaction picture, the interaction Hamiltonian is x + , and the interaction picture state is related to the Schrödinger picture state by |ψ; The solution of Eq. (11) describes the state of the investigated system at any given light-front time x + , where T + denotes light-front time ordering. In perturbative calculations, the time-ordered exponential is written as an expansion in powers of V I (z + ), and is approximated by retaining the leading terms in the series. However, in cases where the external fields are strong, a perturbative treatment may not be sufficient.
One possible nonperturbative treatment is decomposing the time-evolution operator into many small steps of the lightfront time x + , then solving each time step in the sequence numerically, The step size is δx + ≡ x + /n, and the intermediate time is x + k = kδx + (k = 0, 1, 2, . . . , n) with x + 0 = 0 and x + n = x + . This product sequence is equivalent to the time-ordered exponential in the continuum limit n → ∞.
In practice, the calculation is carried out in a finitedimensional basis space, where the state becomes a column vector, and the interaction operator is in matrix form. The choice of the numerical method, to some extent, depends on the basis representation of the system. Here, we consider two typical treatments for general purposes, and we will discuss the numerical method in solving this problem after introducing the basis in the next section.
Knowing that Eq. (11) is an ordinary differential equation, one primary group of numerical methods is the finitedifference method (FDM). FDM approximates the derivatives with finite differences in each small time step. Typical methods of the group include the Euler method, the secondorder difference scheme MSD2 [34], and Runge-Kutta methods [35]. For example, with the most straightforward method, the forward Euler, one would treat the evolution in each time step as This method is, however, not numerically stable since the formula is not invariant under time reversal. For practical use, stable methods such as MSD2 and the fourth-order Runge-Kutta methods are recommended. The Runge-Kutta methods propagate a solution over each step by combining the information from several smaller Euler-style steps and eliminating lower-order errors. Thus it has the advantage of simulating the time dependence even inside each time step. One could adjust the step size δx + to achieve a desired accuracy in the calculation. There are also implicit methods, such as the Crank-Nicholson method, which uses the backward difference in time and is always stable. However, in these cases, one might need to pay the price of inverting the interaction matrix in a large basis space, which is not always an easy task, especially when the interaction matrix is complicated. Another treatment is to compute the exponential directly, which is automatically unitary. When the time step is sufficiently small, the interaction during every single step can be considered as constant in time, and the evolution operator reduces to an ordinary exponential, However, this way, one loses the time dependence of V I (z + ) within each time step. This method would be favorable if the matrix exponential is straightforward to evaluate, which is the case especially when the interaction matrix is diagonal.
These introduced methods all simulate the time evolution by computing the interaction in a sequence of time steps, and they provide a nonperturbative solution. Our algorithm is a combination of the Runge-Kutta method for the gluon emission and absorption and the matrix exponentiation for interaction with the background field, as is explained in detail in Sec. II C 4.

Constructing the basis
We are interested in how the momentum states, i.e., eigenstates of the kinetic energy part of the Hamiltonian P − KE , evolve due to gluon emission/absorption and interactions with a background field. Therefore we choose the basis state |β as the eigenstates of the free Hamiltonian P − i.e., the "bare" 1-and 2-particle Fock states. The quark state is a sum over the basis states where c β (x + ) ≡ β|ψ; x + I are the basis coefficients. The initial state at x + = 0 can be specified by assigning values of c β (0), and the information of a state at x + is encoded in the column vector c(x + ).
In each Fock sector, the many-particle basis states are direct products of single particle states. The basis state in the |qg sector is in the format of β qg = β q ⊗ β g . Each single particle state carries five quantum numbers, The first quantum number, k + l , labels the longitudinal momentum of the particle. For this degree of freedom, we employ the usual plane-wave basis states, i.e., eigenstates of the longitudinal momentum operator P + , with corresponding eigenvalues p + l . In this paper, we compactify x − to a circle of length 2L (i.e., x + to a circle of length L ). We impose (anti-)periodic boundary conditions on (fermions) bosons. As a result, the longitudinal momentum p + l in the basis states takes discrete values as with the dimensionless quantity k + g = 1, 2, 3, . . . for bosons (neglecting the zero mode) and k + q = 1/2, 3/2, 5/2, . . . for fermions.
For each Fock state, let K = l k + l be the total k + of all the l particles in that state. Since the background field that we are considering does not provide extra longitudinal momentum to the state, the total p + of the system and K are conserved. In the |q sector, the quarks in all basis states have k + q = K. In the |qg sector, there are a number of (K − 0.5) K-segments, where in each K-segment, the quark and the gluon have definite values of k + q and k + g . For example, with K = 8.5, the quark in the |q sector has k + = K = 8.5, and the |qg sector comprises eight K-segments, each with {k + q = 0.5, k + g = 8}, {k + q = 1.5, k + g = 7}, . . ., {k + q = 7.5, k + g = 1}, respectively. The next two quantum numbers, k x l and k y l , label the momentum components in the transverse directions. The twodimensional transverse space is a lattice extending from −L ⊥ to L ⊥ in each direction with periodic boundary conditions. The number of transverse lattice sites in each dimension is 2N ⊥ , so the lattice spacing is a ⊥ = L ⊥ /N ⊥ . Thus the transverse coordinate vector r ⊥ = (r x , r y ) is discretized as The corresponding momentum space is also discrete with periodic boundary conditions. The transverse momentum vector p ⊥ = (p x , p y ) on the momentum grid reads where d p ≡ π/L ⊥ is the resolution in the transverse momentum space, which effectively acts as an IR cutoff λ IR = d p . The ultraviolet (UV) cutoff from the transverse momentum grid is λ UV = N ⊥ d p = π/a ⊥ . The transverse coordinate and the transverse momentum spaces are related through the Fourier and the inverse Fourier transformations. For the interaction with the background field, we go to transverse coordinate space, where the basis states are characterized by the quantum numbersβ l = {k + l , n x l , n y l , λ l , c l }, l = q or g .
The basis states are eigenstates of the kinetic energy operator P − KE . For each Fock state, the total kinetic energy sums over all the constituent particles l in that state, P − β = l p − l . The kinetic energy of the quark is p − q = ( p 2 ⊥,q + m 2 q )/p + q and that of the gluon is p − g = p 2 ⊥,g /p + g . The number of basis states N tot for the Fock space |q + |qg is therefore This is the number that controls the overall numerical complexity of the calculation.
In the numerical simulation, we take L ⊥ = 50 GeV −1 (= 9.87 fm) and N ⊥ = 16. Exceptions are separately noted. This translates into a rather large lattice spacing a ⊥ in physical units. In order to stay safe from lattice effects, we must use rather small values of g 2μ and m g in physical units to stay close enough to the continuum, i.e., with Q s a ⊥ π. However, since the actual physical behavior of the system only depends on dimensionless combinations of the parameters, one can directly reinterpret our results as valid for larger values (in physical units) of g 2μ on a correspondingly smaller (in physical units) lattice size L ⊥ . The main purpose of this paper is the development of the numerical method, and while we quote values for the parameters in physical units for convenience, the exact values of the parameters should not be interpreted as precisely matching a specific collision system.

Gluon emission and absorption matrix elements
In the basis space, the quark state is represented as a column vector c(x + ). The interaction operator V I (x + ) is represented as a matrix, which we denote as V(x + ). Each matrix element encodes the transition amplitude between two basis states, Recall that the interaction operator contains two terms, V(x + ) = V qg + V A (x + ) (see discussions in Sec. II A). In constructing the basis representation, we have discretized the transverse and the longitudinal spaces, and the light-front Hamiltonian is quantized on the same discrete space (see further details in Appendix B 3). We write out the matrix element of V qg in the basis representation in this section, and we discuss that of V A (x + ) in the next section. The V qg operator acts between the |q and the |qg sectors. In the following expressions, p l = (p + l , p x l , p y l ) is the three momentum of particle l. The symbol β l denotes the collective quantum numbers defined in Eq. (18), and the relation between the integer (half-integer) momentum quantum numbers k l and their associated momenta p l are given by Eqs. (19) and (21). The interaction operator is The matrix element for a transition from a |q state to a |qg state reads and that for a gluon absorption process is the Hermitian con- Here, u(p, λ) is the spinor of the fermion, and µ (p, λ) is the polarization vector of the vector boson. Their expressions can be found in Appendix. B 2. We use the subscripts "Q" and "q" to distinguish between the quark in the |q state and that in the |qg state. For convenience, let us define the longitudinal momentum fraction of the gluon as z ≡ p + g /p + Q , so that p + g = zp + Q and p + q = (1 − z)p + Q . Let us also define the momentum difference between the quark (gluon) in the |qg state and the quark in the |q state as The spinor-polarization vector contraction parts of the matrix elements in Eqs. (25) and (26) are summarized in Table II in Appendix B 2. They depend on the relative center-of-mass momentum, instead of separately on the single particle transverse momenta p ⊥,l . The energy difference from the phase factor in Eq. (23) also depends on ∆ m = | ∆ m |, Thus the matrix element of V qg does not depend separately on the individual momenta of the particles but on the transferred momentum.
The periodic boundary condition implemented on the transverse momentum grid should also apply to the determination of the momentum conservation δ 2 ( p ⊥,Q − p ⊥,q − p ⊥,g ) on the lattice and the calculation of the transferred momenta ∆ q and ∆ g . Due to the periodicity, p i Q and p i q + p i g (i = x, y) are equal if either they have the same value or they are different by a period in the transverse momentum space, 2λ UV . Consequently, a transition process on the lattice could correspond to more than one different physical process, so one must decide which copy of the periodical momentum space lattice should be used to evaluate the momentum differences ∆ q and ∆ g that determine the matrix element. For example, a |qg state with the quark and the gluon each carrying a transverse momentum close to the boundary λ UV can merge to a |q state with a large total transverse momentum close to 2λ UV , which is outside of the fundamental Brillouin zone. However, on a periodic lattice, we could interpret the same gluon as having a transverse momentum just beyond the opposite boundary −λ UV , merging with a quark close to λ UV into a quark with a momentum close to zero. With the first interpretation, the momentum difference vectors ∆ q and ∆ g point in the same direction, whereas for the second one, they are opposite. Thus, the relative center-ofmass momentum ∆ m = −(1 − z) ∆ q + z ∆ g that determines the matrix element and light-front energy difference will be very different with the two interpretations.
To get rid of ambiguities due to the periodicity, we choose the following prescription. We always use the value of p ⊥,Q within the fundamental Brillouin zone as p ⊥,Q in calculating the quark momentum transfer ∆ q . We then use the value of the momentum sum p ⊥,q + p ⊥,g (which might lie outside of the fundamental Brillouin zone) as p ⊥,Q in calculating the gluon momentum transfer ∆ g . For the configuration discussed above, this corresponds to the second interpretation of a back-to-back |qg state merging into a small momentum |q state. The reason for this choice is precisely to maintain this interpretation of back-to-back splitting and merging, which is the physically most relevant process for the physical situations we are interested in. We discuss the periodic boundary condition and explain our prescription in detail in Appendix C.

Background field interaction matrix elements
The V A (x + ) term is introduced by the chosen background field, and it contains two parts, one acting on the quark and the other on the gluon: Here, the symbol β l denotes the collective quantum numbers defined in Eq. (18), and the relation between the integer (halfinteger) momentum quantum numbers k l and their associated momenta p l are given by Eqs. (19) and (21). The V A (x + ) term does not contain the quantum gauge field and therefore does not directly connect different Fock sectors, so matrix elements of the type qg| V A |q and q| V A |qg are zero. The background field does not change the particle's longitudinal momentum p + l either, so the matrix elements between two |qg states from different K-segments are also zero.
The background field is local in coordinate space, so it is convenient to evaluate the matrix element in the coordinate basis. The matrix element for a transition from a |q basis state to another |q basis state reads The collective basis numberβ l is defined in Eq. (22), and the relation between the basis numbers (k + l , n x l , and n y l ) and their associated momenta/locations (p + l , r x l , and r y l ) are given by Eqs. (19) and (20). The matrix element for a transition from a |g basis state to another |g basis state reads β g (k + g , n x g , n y g , λ g , c g ) V A (x + ) β g (p + g , n x g , n y g , λ g , c g ) = −i2g f ac g c g A a + ( r ⊥,g , x + )δ λ g ,λ g δ k + g ,k + g δ n x g ,n x g δ n y g ,n y g .
The background field A a + (= A −,a /2) is generated from the sampled color charges on the same discretized transverse lattice of the Fock state. The longitudinal dimension of the color charge in x + (note that this is the light-front time of the incident quark) is taken to consist of N η independent layers [33]. The color charge, as well as the generated background field, extend from 0 to L η along x + . Thus each layer has a thickness of τ = L η /N η , with the n τ -th (n τ = 1, 2, . . . , N η ) layer spanning x + = [(n τ − 1)τ, n τ τ]. The correlation relation of the color charge as defined in Eq. (9) now takes a discrete form ρ a (n x , n y , n τ )ρ b (n x , n y , n τ ) = g 2μ2 δ ab δ n x ,n x δ n y ,n y a 2 ⊥ δ n τ ,n τ τ .
The Kronecker delta dividing the discrete resolution would become the Dirac delta in Eq. (9) in the continuous limits of a ⊥ → 0 and τ → 0. For generality, we allow the time step δx + to be smaller than the layer thickness τ; this allows one to continuously go from scattering off a large coherent (independent of x + ) background field to scattering off independent scattering centers represented by separate layers in x + .

Time evolution in the basis
We now look at the time evolution in this basis representation. The solution of the time-evolution equation, Eq. (12), acquires a matrix form. In each time step, the evolution reads where V(x + ) is the interaction matrix in the basis representation, and we have already discussed its matrix element V ββ (x + ) in Secs. II C 2 and II C 3. We could now select a suitable numerical method that takes advantage of the interaction matrix's structure. First, we notice that since our interaction matrix is a sum of two terms, V I (x + ) = V qg,I (x + ) + V A,I (x + ), we can decompose the evolution over an infinitesimally short interval into two successive operations Then we use different numerical methods for the two different kinds of interactions.
The gluon emission/absorption operator V qg is off-diagonal in the Fock space and is thus challenging to exponentiate. Therefore we use the fourth-order Runge-Kutta (RK4) method to calculate its contribution in the time evolution, The explicit form of the RK4 operator U RK4 can be found in Appendix. D. In terms of computational complexity, the RK4 method on the basis space is, in principle, O(N 2 tot ), but it is more like O(N tot ) for V qg . That is because the gluon emission/absorption interaction is nonzero only when the momentum is conserved, so the matrix V qg is very sparse. In practice, we organize the numerical computation to iterate over only the matrix elements allowed by momentum conservation, which achieves this O(N tot ) complexity.
On the contrary, the interaction with the background field V A (x + ) is diagonal in the Fock space: it does not cause transitions between |q and |qg sectors. Moreover, the background field in our simulation is eikonal, meaning that the interaction is diagonal in coordinate space and in helicity space. One only needs to exponentiate a N c × N c color matrix to achieve a unitary evolution over a time step, which can be calculated analytically with the Cayley-Hamilton theorem [37]. Therefore, it is feasible to do the calculation in the exponential form by Fourier transforming the wave function into coordinate space and then back again as Here, F = F ( p ⊥ → r ⊥ ) and F −1 = F −1 ( r ⊥ → p ⊥ ) are the Fourier and the inverse Fourier transformation operators, respectively (see further details in Appendix. B 3). Note that the kinetic energy operator is diagonal in momentum, not coordinate space. Thus the kinetic energy phase part of the interaction picture interaction needs to be evaluated in momentum, not coordinate space. The computational complexity of the kinetic energy part is O(N tot ). The (inverse) Fourier transform is carried out through the fast Fourier transform algorithm, which has a complexity of O(N tot log N tot ), and the interaction in coordinate space is O(N tot ). Thus the overall complexity of the background field interaction is O(N tot log N tot ). The full evolution for each time step combines the two contributions as The total computational complexity of each time step is, therefore, O(N tot log N tot ), much more efficient than the O(N 2 tot ) operations that a momentum space interaction with the background field would be. Thus splitting the interaction into two successive steps by Eq. (35) and using a Fourier transform for the background field allow for a very efficient time-evolution algorithm.

III. RESULTS
By carrying out the explicit time evolution of the state, we are able to access the information about its time development as a function of x + . In this section, we study the time evolution of the quark state by looking into its longitudinal momentum, transverse momentum, helicity, and color.
We simulate three different cases. In the first case, the interaction contains just the gluon emission/absorption term V = V qg ; in the second case, the interaction contains just the background field term V = V A . Finally, we consider the full In the cases with nonzero transitions between the |q and the |qg sectors, we start with an initial condition as a single quark state with a definite color, helicity and momentum. When studying the effect just from the background field, i.e., no transitions between the |q and the |qg sectors, we choose a superposition of a single |q and a single |qg state as the initial state to study their respective evolutions under the interaction. These initial states do not correspond exactly to those in a physical high-energy scattering process, where the quark would have already developed a gluon cloud before the interaction. However, it enables us to test the physical effects of the different parts of the Hamiltonian, and our numerical method, in a cleaner and more tractable setup.

A. Gluon emission and absorption
The interaction V qg excites transitions between the |q and the |qg sectors. This effect is intertwined with the phase rotation generated by the free part of the Hamiltonian P − KE , which in the interaction picture is manifested by the time evolution of the interaction matrix V qg,I (x + ) = e i 1 This phase factor leads to a decoherence between emissions separated by a long enough light-front time. To understand the effects from the gluon emissions/absorptions and the phase factor separately, we run the simulations in two cases: with the phase factor, in which we take V I (x + ) as V qg,I (x + ); and without the phase factor, in which we take V I (x + ) as V qg .
We first study the evolution of the quark state in the longitudinal momentum p + phase space. Figure 2 shows the evolution of the probabilities of different p + states, including the |q sector and the K-segments of the |qg sector characterized by the gluon longitudinal momentum fraction z. The probability of each p + state sums over all states in the transverse momentum space, helicity space, and color space. The initial state of the quark is a single quark state with p ⊥,Q = 0 ⊥ , p + Q = P + = 8.5 GeV, light-front helicity λ Q = 1/2, and color c Q = 1. In Fig. 2(a), in the absence of the phase factor, the system oscillates between the initial |q state and all the p + states in the |qg sector. In addition, those different p + states oscillate with the same frequency but with different amplitudes. In Fig. 2(b), with the phase factor restored, the probability for each p + state behaves as a damped oscillation.
To understand the oscillational patterns observed in the simulation via Fig. 2, we study a simplified two-mode problem analytically. Let us consider the state in a two-dimensional   vector space, corresponding to the two Fock sectors. The state vector reads The interaction operator V qg, By solving the time evolution equation as Eq. (11), we obtain the probabilities of the states as sinusoidal functions of the evolution time: For convenience, we have defined w ≡ |u|, The parameter w corresponds to the magnitude of the gluon emission and absorption term, and ∆ is the energy difference arising from the phase factor. The oscillation frequency depends on both the matrix element and the energy difference, as seen in the expression of η. The oscillation amplitude depends on the ratio of the two terms w 2 /η 2 . This two-mode process is essentially the Rabi oscillation, with a Rabi frequency of 2w and a detuning of ∆ [38,39].
The solution of the two-mode problem in Eq. (41) helps understand the evolution of the extended |q + |qg state in the basis space, which is essentially an N tot -mode problem. Let us consider the transition between one |q state and n different |qg states, and the interaction operator is given by The simulation without the phase factor corresponds to only considering the part V qg of the Hamiltonian. In this case, there are only two nonzero eigenvalues, w ± = ± |u 1 | 2 + |u 2 | 2 + . . . + |u n | 2 , which are opposite to each other. Thus the situation is very similar to the two-mode problem. As a result, each basis state oscillates with the same frequency w = |w ± |, although the amplitudes of those oscillations could be different, depending on the values of the interaction matrix elements u i . The probability for each p + state, summing over different transverse momentum modes, therefore also oscillates with the same frequency.
In the full calculation, u i s are the matrix elements of V qg in Eqs. (25) and (26), and they depend on the transferred momentum. The frequency w is dominated by the most significant transition mode, so it is approximately w ∝ λ 2 UV /P + , where P + is the total longitudinal momenta of the state, and λ UV is the largest allowed transverse momentum on the lattice.
The probability of the quark staying in the |q sector at (a) various N ⊥ (P + = 8.5 GeV) and (b) various P + (N ⊥ = 16). The interaction contains just the gluon emission/absorption term V(x + ) = V qg , and phase factor is not included. The initial state of the quark is a single quark state with p ⊥,Q = 0 ⊥ , p + Q = P + , light-front helicity λ Q = 1/2, and color c Q = 1. Parameters in these simulations: L ⊥ = 50 GeV −1 , L η = 50 GeV −1 , m g = 0.1 GeV, m q = 0.02 GeV, K = 8.5. The duration of each time step in the simulation is δx + = 0.39 GeV −1 . Figure 3 shows the evolution of probability of the |q sector at different λ UV (= N ⊥ π/L ⊥ ) by taking different N ⊥ at a fixed L ⊥ , and at different P + . The dependence of the oscillation frequency on P + and λ UV indeed agrees with the expectation w ∝ λ 2 UV /P + . When the phase factor is restored, this corresponds to including both the V qg and the P − KE terms in the Schrödinger picture light-front Hamiltonian P − . Unlike in the case without the phase factor, there are now n + 1 different eigenvalues. Each basis state is, in essence, a superposition of different eigenstates. The summation over these states leads to decoherence, which appears in Fig. 2(b) as a damped oscillation. The probability of each p + state approaches an asymptotic value, which is related to the matrix elements of the Hamiltonian.
We present the probability distribution of the quark state in the p + space after the evolution in Fig. 4. The initial state is a single quark state, and those |qg states with different p + configurations emerge during the evolution. As we see from the result, the gluon emission/absorption process favors |qg states with either small or large z, the gluon longitudi- FIG. 4. The probability of the quark state at different p + configurations after the evolution, with various quark masses. The interaction is V(x + ) = V qg , and the phase factor is not included, i.e., The initial state of the quark is a single quark state with p ⊥,Q = 0 ⊥ , p + Q = P + , light-front helicity λ Q = 1/2, and color c Q = 1. Parameters in these simulations: nal momentum fraction. The dependencies on z and the quark mass can be understood by examining the spinor-polarization vector contractions in the matrix elements of V qg , as in Table. II in Appendix B 3. Quark light-front-helicity-conserving transitions to both gluon polarization states, , are enhanced at small gluon momentum fraction z. Overall emissions in this soft gluon limit, where the emission matrix element is independent of the gluon polarization, are the most likely ones. For large gluon momentum fraction z, on the other hand, the only surviving quark light-front-helicity-conserving emissions are the ones where also the gluon has the same helicity as the quark, [↓→↑↓], are proportional to the quark mass, and heavily weight large values of z, which can be seen in a comparison of the different mass results in Fig. 4.
Next, we study the evolution of the quark state in the transverse momentum space. Figures 5 and 6 demonstrate the probability distributions in the transverse momentum plane for successive times. The transverse momentum distributions are shown separately for the quark in the |q sector and the quark and the gluon in the |qg sector. Now that we do not have a background field, the total transverse momentum is conserved. Thus, the quark in the |q sector stays in its initial momentum state p ⊥,Q = 0 ⊥ . The quark and the gluon in the |qg sector are back-to-back in momentum and have distributions that are symmetric around the origin, apart from the edges of the discrete transverse momentum lattice, where rotational invariance is lost. The distributions without the phase factor are shown in Fig. 5. Here, the emitted quark and gluon both favor large transverse momentum modes, as we see in the sequential distributions in Figs. 5(b) and 5(c). The probabilities of different transverse momentum modes in the |qg sec-tor oscillate coherently, so they maintain their relative magnitudes while rising and falling as functions of x + through the evolution, as seen in Fig. 2(a). When the phase factor is included, the emitted quark and gluon show a changing concentric circular pattern in transverse momentum space, as we see in the sequential distributions in Figs. 6(b) and 6(c). As we have discussed earlier in the context of the evolution of different p + states, here different transverse momentum states in the |qg sector are also different superpositions of the eigenstates. Thus, the probabilities of different transverse momentum modes in the |qg sector do not oscillate coherently, and their relative magnitudes change through the evolution. Additionally, the oscillation frequency of each eigenstate depends on the change of the light-front energy and the value of V qg , both depending on the transferred momentum squared. This explains why the pattern reflecting the relative magnitudes among different transverse momentum states is azimuthally symmetric and centered at the initial momentum mode of the quark. The states at later times exhibit artificial effects from the periodic boundaries, and they are not presented here.
To see the effect of the phase factor more clearly, we take the ratio of the probability distribution with the phase factor over that without the phase factor. Since both the V qg interaction and the phase factors are azimuthally symmetric in the transferred p ⊥ plane, we analyze the evolution of the p ⊥ distribution at p ⊥ = | p ⊥ |, arg p ⊥ = 0, π. We set the initial state as a single quark state with p ⊥,Q = 0 ⊥ , and run the simulations with and without the phase factors at various P + . The results are shown in Fig. 7.
In the left panels, Figs. 7(a), 7(c), and 7(e), the probability distributions of the quark in the |qg sector are shown as a function of p ⊥ = | p ⊥ | (arg p ⊥ = 0, π) at a sequence of x + s, with P + = 85 GeV, 8.5 GeV, 4.25 GeV, respectively. The distributions with the phase factor, as in the solid lines, show oscillational patterns, compared to those without the phase factor, as in the dashed lines. In the plots of the ratio of the probability with the phase factor over that without the phase factor, as in Figs. 7(b), 7(d), and 7(f), there is a peak around zero momentum transfer, and it gets narrower over time. (One exception is the x + = 25 GeV curve at P + = 4.25 GeV in Fig. 7(f): there is a dip instead of a peak in the center. But this is caused by artificial reflections from the periodic boundary, so we should neglect it for the purpose of this discussion.) By comparing the three different P + cases, one can see that the peak narrows faster at a smaller P + . Moreover, the peak develops at a rate inversely proportional to P + , which can be seen by comparing the x + = 12.5 GeV (x + = 25 GeV) curve in Fig. 7(d) to the x + = 6.25 GeV (x + = 12.5 GeV) curve in Fig. 7(f). This is because a smaller P + leads to a larger kinetic energy P − KE ∝ 1/P + , making the decoherence faster. This behavior is a demonstration of the familiar effect leading to Fermi's golden rule. At late times x + → ∞, the only allowed transitions are the ones that conserve the light-front energy P − . This energy conservation is enforced by the phase factor, canceling the energy nonconserving transitions, even when they are favored by large transition matrix elements.
We then look at the evolution of the quark state in color space, as in Fig. 8. The initial state here is a single quark with  color index c Q = 2. Only six of the |qg color states are allowed in the transitions due to color conservation. Without the phase factor, the probabilities of those states oscillate over time, as in Fig. 8(a). The oscillation is suppressed when the phase factor is restored, as in Fig. 8(b). This oscillation and its suppression have the same reason as in the p + distribution shown in Fig. 2. Without the phase factor, the probability of each momentum mode oscillates coherently, so the probability of each color state, which is a summation over all the momentum modes, also oscillates coherently. However, with the phase factor, different momentum modes oscillate with different frequencies, eventually going out of phase. Thus, the probability of a color state, which sums over all the momentum modes, even acquiring an oscillation initially, could not maintain it.
Lastly, we examine the evolution of the quark state in he-licity phase space. The results are presented in Fig. 9. As the evolution time increases, states in the |qg sector appear. Since the initial state is a single quark state with λ Q = 1/2, the produced |qg states favor the {λ q = 1/2, λ g = ±1} configurations in which the quark helicity is preserved. The transition to the {λ q = −1/2, λ g = 1} state is weighted by the quark mass, which is relatively small in this case. The {λ q = −1/2, λ g = −1} state is not allowed. Very much like the evolution of probability distribution in the color space, the probabilities of those helicity states oscillate over time when the phase factor is not included, as in Fig. 9(a), and the oscillations are suppressed when the phase factor is restored, as in Fig. 9(b).
From the above results and discussions, we see that the evolution with the gluon emission/absorption interaction contains two contributions: the transition between the |q and the |qg sectors by V qg , and a phase rotation by P − KE . This interaction preserves the system's total momentum, and it changes the distribution of the state in both the p + and the p ⊥ spaces, as well as in color and helicity spaces. Without the phase factor, the transitions happen as coherent oscillations between different states, but the phase factor causes the transitions to decohere.

B. Interaction with background field
In this section, we study the effect from the background field without gluon emission or absorption. The background field interacts with the |q and the |qg sectors separately, and does not in itself cause transitions between them. The interaction with just the |q sector was previously studied with the tBLFQ approach in Ref. [24]. The background field in the simulation has P + A = 0, so it does not change the p + configuration of the system. The nonzero component of the background field is A − , which couples to the J + current of the fermion field, so the light-front helicity of the quark state is not affected either. The background field only affects the distributions in transverse momentum space and in color space.
We present the evolution of the quark state in transverse momentum space in two cases: one with a relatively weaker field with g 2μ = 0.018 GeV 3/2 in Fig. 10, and the other with a relatively stronger field with g 2μ = 0.108 GeV 3/2 in Fig. 11. In both cases, the initial state is a superposition of a |q state with p + Q = P + = 8.5 GeV, p ⊥,Q = 0 ⊥ , helicity λ Q = 1/2, color index c Q = 1 and a |qg state with p + q = 0.5 GeV, p + g = 8 GeV, p ⊥,q = p ⊥,g = 0 ⊥ , helicity λ q = 1/2, λ g = 1, and color index c q = 1, c g = 1. The basis coefficient for each of the two is 1/ √ 2. The total evolution time of the presented results is L η = 25 GeV −1 . The typical transverse momentum that the particles obtained from the background field is characterized by the saturation scale Q s , as defined in Eq. (10). In both simulations, the values of Q s are far below the UV cutoff of the grid λ UV = π/a ⊥ so that the calculated result is close to the continuum limit and away from the lattice effects. The values of the dimensionless quantity Q s a ⊥ in the two cases are 0.13 and 0.78, respectively, both sufficiently smaller than π. As we see in Figs. 10 and 11, the majority of the occupied momentum modes are still away from the boundary of the transverse momentum lattice by the end of the evolution.
Under the interaction with the background field, both the initial |q state and the initial |qg state transfer to other momentum modes within their Fock sector. This momentum transfer is more obvious with the stronger field in Fig. 11 compared to that in Fig. 10. The circular pattern resulting from the phase factor appears in the transverse momentum distribution. Because of its relatively small longitudinal momentum p + q = 0.5 GeV, compared to the total P + = 8.5 GeV of the system, the quark in the |qg has a more significant phase rotation from the phase factors. Thus, the circular pattern is most noticeable for the quark in the |qg sector when the background field is weak, as in Fig. 10(b). By comparing the transverse momentum distribution of the gluon and that of the quark, one sees, especially with the stronger field in Fig. 11, the effect of Casimir scaling; because C A > C F the gluon gets a stronger momentum kick from the background field than the quark.
Then, we look at the evolution of the quark state in color space, as in Fig. 12. The initial state is a superposition of a |q state with color c Q = 1 and a |qg state with color c q = 1, c g = 1. The basis coefficient for each of the two is 1/ √ 2. The interaction V A carries out a color rotation within the |q and within the |qg sector, separately. In the two cases with and without the phase factor, as in Figs. 12(a) and 12(b), all color states emerge during the evolution, and the state approaches a uniform color distribution in each Fock sector.

Cross sections
The interaction of a particle with the background field is usually quantified in terms of the cross section for scattering off the field. The study in Ref. [24] calculated the cross section of a pure |q state under the CGC background field. In this work we study the cross section of a pure |qg state. These studies would get us prepared for calculating the cross section of a QCD eigenstate in the |q + |qg Fock space in the future.
The cross section of a process is defined as the sum of the , and phase factor is included. The initial state is a superposition of a |q state with p + Q = P + = 8.5 GeV, p ⊥,Q = 0 ⊥ , helicity λ Q = 1/2, color index c Q = 1, and a |qg state with p + q = 0.5 GeV, p + g = 8 GeV, p ⊥,q = p ⊥,g = 0 ⊥ , helicity λ q = 1/2, λ g = 1, color index c q = 1, c g = 1. The basis coefficient for each of the two is 1/ √ 2. From left to right, the transverse momentum distributions of the particle are shown at increasing light-front time instances. The number at the bottom of each panel is the total probability of the plotted states. Parameters in the simulation: m g = 0.1 GeV, N ⊥ = 16, L ⊥ = 50 GeV −1 , g 2μ = 0.018 GeV 3/2 , m q = 0.02 GeV. The duration of each background field layer is τ = 12.5 GeV −1 and that of the each time step in the simulation is δx + = 0.39 GeV −1 . For the rightmost panels, which are at the last of the evolution, the total evolution time is L η = 25 GeV −1 ; the value of the dimensionless quantity Q s a ⊥ [Q s is defined in Eq. (10)] is 0.13. squares of the transition amplitudes, Here, ψ i stands for the initial state and φ f the final state; φ f sums over the phase space of the final state. The S in the equation is the evolution operator from the initial state to the final state.
In the usual case corresponding to a physical scattering experiment, the time evolution happens over an infinite interval from x + = −∞ to x + = ∞. For a finite-size target, this allows for an incoming quark to develop a cloud of gluons before the target and for the Fock states of the scattered particle to reorganize after the target through the V qg interaction. In our explicit numerical time-evolution procedure, such an infinite time evolution would not be feasible. Instead, we initialize our system in a specific Fock state at the time x + = 0 and study the evolution within the target color field for a finite time L η . Thus the calculation we are doing here does not actually correspond to scattering, and the quantity defined by Eq. (43) should not be interpreted as a usual cross section. Studying a physical scattering process is possible with the same timeevolution algorithm. However, it requires using initial condi-    11. The evolution of the transverse momentum distributions of (a) the quark in the |q sector, (b) the quark in |qg sector, and (c) the gluon in |qg sector. The interaction contains just the background interaction term V(x + ) = V A (x + ), and phase factor is included. The initial state is a superposition of a |q state with p + Q = P + = 8.5 GeV, p ⊥,Q = 0 ⊥ , helicity λ Q = 1/2, color index c Q = 1, and a |qg state with p + q = 0.5 GeV, p + g = 8 GeV, p ⊥,q = p ⊥,g = 0 ⊥ , helicity λ q = 1/2, λ g = 1, color index c q = 1, c g = 1. The basis coefficient for each of the two is 1/ √ 2. From left to right, the transverse momentum distributions of the particle are shown at increasing light-front time instances. The number at the bottom of each panel is the total probability of the plotted states. Parameters in those panels: m g = 0.1 GeV, N ⊥ = 16, L ⊥ = 50 GeV −1 , g 2μ = 0.108 GeV 3/2 , m q = 0.02 GeV. The duration of each background field layer is τ = 12.5 GeV −1 and that of the each time step in the simulation is δx + = 0.39 GeV −1 . For the rightmost panels, which are at the last of the evolution, the total evolution time is L η = 25 GeV −1 ; the value of the dimensionless quantity Q s a ⊥ [Q s is defined in Eq. (10)] is 0.78.
tions at x + = 0 with a Fock state with a fully developed gluon cloud that corresponds to an incoming quark at x + = −∞, projecting out to similar scattering states at the end of the target. In this paper, we focus on understanding the interaction within the target and leave the description of the correct asymptotic states to future work. Note that this issue did not concern the earlier tBLFQ calculation with the bare quark in Ref. [24], since in the absence of gluon radiation the time development between x + = ±∞ and the target is trivial; thus, the results of that work could indeed be understood as quark-nucleus scattering cross sections.
In evaluating the total cross section, one should average over the color charge density ρ of the target as in Eq. (9), Here, the . . . stands for a configuration average of the background field. The total cross section includes a projection to the final state at the amplitude level, and a summation over all possible states at the cross section level [31,40,41]. Using the unitarity of the S -matrix, i.e.,the optical theorem, the total cross section can also be expressed in terms of the expectation value of the diagonal elements of the scattering amplitude, i.e.,in terms of the imaginary part of the forward elastic amplitude.
Since the background field interacts with both the quark and the gluon in the |qg state, we first study their respective effects. In the eikonal limit of p + = ∞, both the single quark cross section dσ q / d 2 b and the single gluon cross section dσ g / d 2 b reduce to traces of Wilson lines and can be written in terms of the charge density g 2μ , the interaction duration L η , and the IR cutoff m g [31]. The total cross section of a single quark interacting with the background field is (see Appendix E for detailed derivations of Wilson line expectation values) and that of a single gluon is where C F = (N 2 c − 1)/(2N c ) = 4/3 and C A = N c = 3. Here, one uses the representation in terms of the forward elastic amplitude, which in this case is the expectation value of a single Wilson line. Note that in the CGC picture, this means that the total cross section is the expectation value of a nonsinglet Wilson line operator, and as a consequence, it very strongly depends on the IR cutoff provided by m g . Thinking differentially in terms of the momentum transfer from the target, it includes, besides the finitek ⊥ cross section, which results from the Fourier transform of the color singlet dipole operator, the part at k ⊥ = 0 ⊥ that is not singlet (see similar considerations in, e.g., Refs. [42,43]). Now let us look at the cross section of a |qg state. The same background field would interact with the quark and the gluon, and the total cross section in the eikonal limit reads dσ qg,tot The calculation of this product of quark and gluon Wilson lines is discussed in Appendix E. Here, f qg (ξ) is a correlation function between the quark and the gluon f qg (ξ) = [7 cos(ξ/2) + cos(3ξ/2)]/8 (see derivation in Appendix E). In the argument of f qg as in Eq. (47), K 1 is the modified Bessel function of the second kind, and x ⊥ and y ⊥ are the transverse coordinates of the quark and the gluon. Unlike the cross sections of the single particle, which are independent of the transverse coordinates, the quark-gluon cross section has a nontrivial dependence on their difference | x ⊥ − y ⊥ |. In the limit | x ⊥ − y ⊥ | → ∞, the Wilson lines seen by the quark and the gluon become uncorrelated. In this limit, Bessel function K 1 approaches zero, and with f (0) = 1, the quark-gluon cross section reduces to dσ qg,tot This is just the product of the single quark and the single gluon cross sections, i.e., the case where the quark and the gluon interact with uncorrelated background fields separately. We ran the simulations with a single |qg initial state under the interaction with the background field at various g 2μ , and we calculated the total cross sections according to Eq. (44). Following the above discussions, we studied four different cases of the interaction: (1) the background field interacts with just the quark, i.e., V(x + ) = V A,q (x + ), (2) the background field interacts with just the gluon, i.e., V(x + ) = V A,g (x + ); (3) the same background field interacts with both the quark and the gluon, i.e., V( (4) different background fields interact with the quark and the gluon, i.e., V(x + ) = V A,q (x + ) + V A ,g (x + ), where A and A are independently generated background fields in the simulation. We present the results in Figs. 13, 14, and 15. In these simulations, the initial state is a single quark-gluon state with p ⊥,q = p ⊥,g = 0 ⊥ , light-front helicity λ q = 1/2, λ g = 1, and color c q = 1, c g = 1. The phase factor is not included, which is equivalent to taking P + = ∞. As studied in Ref. [24], for the evolution of a single quark state with the chosen background field, the total cross sections at finite P + do not show noticeable differences from the P + = ∞ case. We find it also true for the quark-gluon state by running simulations with various P + .   Figure 13 shows the calculated cross sections for the first two cases. The result of the background field interacting with just the quark (gluon) in the |qg state agrees with the eikonal expectation of a single quark (gluon) separately scattering on the background field in Eq. (45) [Eq. (46)], as one would expect. These calculations help check the correctness of our numerical calculations, and they might also be helpful to study processes involving a single quark or gluon. Figure 14 shows the calculated cross section for the latter two cases. In Fig. 14(a), the total quark-gluon cross section of the quark and the gluon interacting with different background fields agrees with the eikonal prediction of the uncorrelated scattering in Eq. (48), as one would expect. However, even in the case where the quark and the gluon inter- FIG. 14. The total cross sections of the |qg state, as functions of g 2μ , evaluated at various N ⊥ . (a) Quark and gluon interact with different background fields, i.e., V(x + ) = V A,q (x + ) + V A ,g (x + ), in which A and A are independently generated background fields for each simulation; (b) Quark and gluon interact with the same background field, i.e., V(x + ) = V A,q (x + ) + V A,g (x + ). In both panels, the solid line is the uncorrelated eikonal prediction as calculated from Eq. (48). The initial state is a quark-gluon state with p ⊥,q = p ⊥,g = 0 ⊥ , lightfront helicity λ q = 1/2, λ g = 1, and color c q = 1, c g = 1. The phase factor is not included in these simulations, i.e., V I (x + ) = V(x + ). Each data point of the total cross section is calculated according to the definition in Eq. act with the same background field, which is more likely to happen in a dressed quark scattering process, the total cross section agrees with this uncorrelated prediction in Eq. (48) as well. In other words, the correlation between the quark and the gluon through interacting to the same background field is too small to be noticeable in the total cross section. To see this quantitatively, at the strength g 2μ where the correlation is strong, the cross section is already close to its black disc limit of σ tot → 2, so the correlation is of little account (see more discussions in Appendix E) .
To get an impression of the relative magnitude of the four cases discussed above, we put them together in Fig. 15 for comparison. The cross section as a function of g 2μ saturates . The data points in this figure are taken from the N ⊥ = 16 results of the four panels in Figs. 13 and 14, leaving out the uncertainties. The red solid, blue dashed, and yellow dotted lines are the eikonal predictions as calculated from Eqs. (45), (46) and (48), respectively. The initial state is a single quark-gluon state with p ⊥,q = p ⊥,g = 0 ⊥ , light-front helicity λ q = 1/2, λ g = 1, and color c q = 1, c g = 1. The phase factor is not included in these simulations, i.e., V I (x + ) = V(x + ). Parameters in these simulations: most rapidly for a |qg state, second for a gluon state, and last for a quark state, also seen from their corresponding eikonal expectation in Eqs. (45), (46), and (48). From the above results and discussions, the physical picture is that the interaction with the background field changes the distribution in transverse momentum space and color space. We also see that the cross section of a |qg state agrees with the eikonal expectation, and the correlation between the quark and the gluon is significantly suppressed in the total cross section defined by Eq. (43).

C. Emission, absorption and background field
Having studied the gluon emission/absorption and the background field separately in the Secs. III A and III B, we now put the two together to have the full interaction V( We consider the initial state of the quark as a single quark state with p ⊥,Q = 0 ⊥ , light-front helicity λ Q = 1/2, and color c Q = 1. The transition probabilities of the quark to other states are shown in Fig. 16. The probability of the quark in its initial state is in the yellow solid line, that of the other |q states is in the blue dashed line, and that of the |qg states is in the red dotted line. When the background field is absent, the probabilities of the states in the |q sector that are different from the initial state are always 0, as shown in Fig. 16(a). With the full interaction, the result shows the combined effects from the gluon emission/absorption and the interaction with the background field. When the background field is relatively  16. The probability of the quark staying in its initial state and the transition probabilities to other states during the evolution. The initial state of the quark is a single quark state with p ⊥,Q = 0 ⊥ , p + Q = P + = 8.5 GeV, light-front helicity λ Q = 1/2, and color c Q = 1. The probability of the quark in its initial state is in the yellow solid line, that of the other states also in the |q is in the blue dashed line, and that of the |qg states is in the red dotted line. From top to bottom: (a) without background field, i.e., g 2μ = 0, (b) with a relatively weak background field, g 2μ = 0.018 GeV 3/2 , and (c) with a relatively strong background field, g 2μ = 0.144 GeV 3/2 . The simulations in the left panels do not include the phase factors, and those in the right panels do. Parameters in these panels: N ⊥ = 16, L ⊥ = 50 GeV −1 , L η = 50 GeV −1 , N η = 4, m g = 0.1 GeV, m q = 0.02 GeV. The duration of each background field layer is τ = 12.5 GeV −1 and that of the each time step in the simulation is δx + = 0.39 GeV −1 .
weak, the result resembles the case with emission and absorption only, but different |q states also emerge, see Fig. 16(b). With a stronger background field, the probability of different |q states is larger, see Fig. 16(c).
The evolution of the quark state in the p + phase space is shown in Fig. 17. The result is very similar to that without the background field in Fig. 2, since the change in p + results from V qg and not from the background field interaction.
The evolution of the quark state in the transverse momentum space is shown in Fig. 18 at g 2μ = 0.018 GeV 3/2 and in Fig. 19 at g 2μ = 0.144 GeV 3/2 . Circular patterns appear as a result of the phase rotation, similar to those in the cases with the gluon emission/absorption in Fig. 6 and those with the background field in Figs. 10 and 11. In addition, transi-  tioning to other momentum modes in both the |q and the |qg sectors appear, resulting from the interaction with the background field. This effect is more obvious with the stronger field in Fig. 19 compared to that in Fig. 18.
The evolution of the quark state in the color phase space is shown in Fig. 20. The initial state is a bare quark with color index c Q = 2. The V qg interaction allows the transition to six of the |qg color states, as we have seen in Fig. 8. The V A interaction allows the color transitions within the |q sector and within the |qg sector, as we have seen in Fig. 12. As a result, all color states emerge during the evolution in Fig. 20. Similar to the evolution with just the V qg interaction, the probabilities of those states oscillate in the simulation without the phase factor, as in Fig. 20(a), and the oscillation is suppressed when the phase factor is restored as in Fig. 20(b).
The evolution of the quark state in helicity space is shown in Fig. 21. The result is very similar to that in the V qg evolution in Fig. 9, since the change in helicity results from V qg and not from the background field interaction.
To sum up this section, we have studied the evolution with the full interaction V qg + V A , where the former is in charge of gluon emission/absorption, and the latter term controls the transitions within each of the |q and the |qg sector. By adjusting the relative magnitude of the two, one could access different physics regimes. From the nonperturbative time evolution, we investigate the combined effects from the full interaction in the quark phase space, including the longitudinal momentum, the transverse momentum, helicity, and color spaces. By adjusting the strength of the background field, one is able to change the relative importance of the gluon emission and absorption, and the color decoherence and momentum broadening due to the background field. Our results overall are consistent with the expectations from having the two different kinds of interactions separately. The interaction contains both the gluon emission/absorption and the background interaction term V(x + ) = V qg + V A (x + ), and the phase factor is included. The initial state of the quark is a single quark state with p ⊥,Q = 0 ⊥ , light-front helicity λ Q = 1/2, color index c Q = 1. From left to right: the transverse momentum distributions of the particle are shown at increasing light-front time instances. The number at the bottom of each panel is the total probability of the plotted states. Parameters in the simulation: m g = 0.1 GeV, N ⊥ = 16, L ⊥ = 50 GeV −1 , g 2μ = 0.018 GeV 3/2 , m q = 0.02 GeV. The duration of each background field layer is τ = 12.5 GeV −1 and that of the each time step in the simulation is δx + = 0.39 GeV −1 . For the rightmost panels, which are at the last of the evolution, the total evolution time is L η = 25 GeV −1 ; the value of the dimensionless quantity Q s a ⊥ [Q s is defined in Eq. (10)] is 0.13.   Probability

IV. CONCLUSIONS AND OUTLOOK
In this work, we developed a numerical implementation of the time-evolution Hamiltonian formalism, tBLFQ, for the interactions of a |q + |qg system with a target color field. Our formulation enables us to access the wave function of the quark at any intermediate time during the evolution, and to continuously tune the relative importance of the interaction with the target field, and of the gluon emissions and absorptions, without taking any parametric limits.
We carried out explicit time evolutions of the quark as a quantum state inside the background color field. Our calculation enables us to access explicitly the time evolution of the transverse and longitudinal momentum, color, and helic-ity of the scattering partons. The light-front Hamiltonian of our system consists of three parts: the kinetic energy term, which leads to a phase rotation of the state, the interaction with the background field, and the gluon emission/absorption. We studied these effects both individually and in combination. The simulations were done for three different cases: the gluon emission/absorption alone, the interaction with the background field alone, and the full interaction that combines the previous two terms. We also compared the processes with and without the phase rotation from the kinetic energy term. Overall, in the limiting cases, the results correspond qualitatively and quantitatively to what one could expect based on general physical arguments, or explicit calculations. We therefore believe that our numerical method is now well tested and robust to be applied to several different physical situations.
In this paper we have focused on developing and testing the numerical method. In the future, as discussed in the Introduction, our goal is to apply this numerical method to different physical situations, such as jet quenching in a hot plasma and a high-energy scattering with subeikonal effects. In this work, we use a single |q state or a single |qg state with definite momentum to study the dynamical process in a simplified yet clean picture. Specific physical applications require initial conditions that are matched to the studied physical system, and calculations of the physical observables that are of interest. For the case of high-energy scattering, one needs as an initial condition a dressed quark state formulated in a way that is consistent with our truncation of the Fock space. In addition to the perturbative calculation of this state, another possibility would be to solve the eigenvalue equation with the QCD Hamiltonian in our truncated Fock space. In this work, we take the background field of the nucleus as the MV model, and keep the dominant field component at high energy (A − ) in our calculation. For the purposes of understanding subeikonal effects and the role of spin in high-energy scattering, it would be interesting to generalize this to a background field with transverse components [44]. In a separate physical situation from that of high-energy scattering, our calculation provides a systematic way to study the interactions of an energetic parton in a colored medium, which is the situation in jet quenching, when a highly energetic parton interacts with a colored medium. Many calculations of jet quenching are done in the approximation of independent static scattering centers. We hope that our formulation would provide for a way to generalize this and enable an understanding of jet quenching in a more general nonperturbatively strong gluonic field configuration, such as the one provided by the pre-equilibrium glasma fields in the initial stage of a heavy ion collision. The light-front coordinates are defined as (x + , x − , x 1 , x 2 ), where x + = x 0 + x 3 is the light-front time, x − = x 0 − x 3 the longitudinal coordinate, and x ⊥ = (x 1 , x 2 ) the transverse coordinates. In this paper, we also use "x" and "y" as the transverse indices, and they should be understood the same as the indices "1" and "2" introduced here. The nonvanishing elements of the metric tensor are The Dirac matrices are four unitary traceless 4 × 4 matrices, where,σ Appendix B: The light-front Hamiltonian

Derivation of the light-front QCD Hamiltonian with a background field
In this section, we derive the light-front QCD Hamiltonian according to Ref. [28] but with an additional background field. The QCD Lagrangian with a background field is given in Eq. (1), The equation of motion for the gauge field gives the color-Maxwell equation, with the current density J κ s ≡ f sac F κµ a C c µ + Ψγ κ T s Ψ. In the light-cone gauge of A + a = A + a = 0, the κ = + component of Eq. (B2) does not contain time derivatives and can be written as By disregarding the zero modes [45], one inverts the equation to We define the free solutionÃ µ a such that lim g→0 A µ a =Ã µ a . According to Eq. (B4), the free field reads The equation of motion for the fermion field gives the color-Dirac equation, We now separate the dynamical component of the fermion field by introducing projectors Λ ± = γ 0 γ ± /2. The projected spinors are thereby Ψ ± = Λ ± Ψ, and we obtain a coupled set of spinor equations from Eq. (B6), Equation (B8) does not contain time derivatives and can be written as a constraint relation, By substituting Eq. (B9) into Eq. (B7), we get In analogy to the free solutionÃ, we define the free spinor Ψ =Ψ + +Ψ − with It is also easy to see thatΨ ± = Λ ±Ψ . The conjugate momenta are We now turn to the construction of the canonical Hamiltonian density through a Legendre transformation, It is convenient to add a total derivative −∂ κ (F κ+ s A s + ) to the Hamiltonian P − = 2P + , We eliminate the light-front time derivatives of the fields by applying the equations of motions in Eqs. (B2) and (B6), and rewrite the full light-front Hamiltonian in terms of only the "tilde" variables defined in Eqs. (B5) and (B11). We introduce the current density of free fields solutionJ µ a in analogy to J µ a , asJ µ s ≡ f sac F µκ aC c κ +Ψγ µ T sΨ , and notice that their "+" components are the same, Finally, we get the light-front Hamiltonian with the background field as The two terms in the first line are the kinetic energy for the gauge field, the background field, and the fermion field. The four terms in the second line can be written collectively as gJ µ a C a µ , which include the three-gluon interaction and the vertex interaction; the latter is responsible for the gluon emission and quark-antiquark-pair-production processes. The two terms in the third line are the instantaneous-gluon interaction and the four-gluon interaction respectively. The last line contains the instantaneous-fermion interaction. For each interaction involving the gluon field, it also involves the background field. Since we are interested in the interactions introduced by the background field to the quark but not the dynamics of the background field itself, we thereby neglect the kinetic energy of the background field and its self interaction in this work.
In the text, we drop the tilde on all variables to simplify the notations, but their meanings are not changed.

Spin and polarization
We use the following spinor representation. The u, v spinors are defined as, and v(p, λ = The polarization vectors for gluon are defined as where ⊥ ± = (1, ±i)/ √ 2. The spinor-polarization vector contraction partū(p Q , λ Q ) γ µ u(p q , λ q ) µ (p g , λ g ) and its complex conjugates are summarized in Table II for different helicity configurations.

Quantization in a discrete space
We consider that the system is contained in a box of finite volume Ω = 2L(2L ⊥ ) 2 . We have introduced two artificial length parameters, L in the longitudinal direction and L ⊥ in transverse directions. In the longitudinal direction, −L ≤ x − ≤ L, we impose periodic boundary conditions for bosons and antiperiodic boundary conditions for fermions such that the longitudinal momentum space is discretized as, , . . . , ∞ for fermions , 2π L k + , with k + = 1, 2, . . . , ∞ for bosons . (B20) In the transverse dimension, −L ⊥ ≤ x 1 , x 2 ≤ L ⊥ , we impose the periodic boundary conditions and discretize the space into 2N ⊥ × 2N ⊥ grids. The corresponding momentum space is also II. Spinor-polarization vector contraction for different helicity configurations. For any transverse two-dimensional vector, p ⊥ = (p x , p y ), define p R ≡ p x + ip y , and p L ≡ p x − ip y . As defined in Sec. II C 2, z ≡ p + g /p + Q is the longitudinal momentum fraction of the gluon, and ∆ m is the relative center-of-mass momentum defined in Eq. (28).

(B25)
The mode expansion for field operators on such a discrete momentum basis is where p · x = p + x − /2 − p ⊥ · x ⊥ is the 3-product for the spatial components of p µ and x µ . Each single particle state is specified by five quantum numbers,ᾱ = {k + , k 1 , k 2 , λ} and c (a), where λ is the light-front helicity, and c (a) is the color index.
Note that this is the same with the basis number β = {ᾱ, c} defined in our basis representation. The creation operators b † α,c , d † α,c and a † α,a create quarks, antiquarks, and gluons with their corresponding quantum numbers, respectively. They obey the following commutation and anticommutation relations: The fields obey the standard equal-light-front-time commutation relations, and here we write it out for the dynamical fields: in which Λ + = γ 0 γ + /2 is the same light-front projector introduced in Appendix. B 1, and with i, j = 1, 2, and (x) is the sign function. A single quark basis state is defined as The basis in the transverse coordinate space that is related to it through Fourier transformation is defined as We define single gluon basis states as and β g (k + g , n 1 g , n 2 g , λ g , c g ) = k 1 g ,k 2 g e i(n 1  22. An example of a quark (denoted as "q") and a gluon (denoted as "g") transferring into/from a quark (denoted as "Q") with their momenta satisfying k q + k g = k Q . This is a nonproblematic case when we do not need to worry about the periodicity. The grids inside the fundamental Brillouin zone, i.e., those with momentum numbers in the range of [−N ⊥ , N ⊥ − 1], are in solid lines. The grids outside this range are in dashed lines. In (a), particles are marked at their momentum quantum numbers assigned on the lattice. In (b), particles are marked at their momentum quantum numbers used to calculate the transferred momenta ∆ q and ∆ g .
ways in the fundamental Brillouin zone. We resolve this ambiguity by making consistent choices in matching the physical process and the process calculated on the lattice. Let us first look into each of the three momentum-conserved cases of the qg ↔ Q transition, i.e., k q + k g = k Q , k q + k g = k Q + 2N ⊥ , and k q + k g = k Q − 2N ⊥ , separately.
1. k q + k g = k Q Since the sum k q + k g is already inside the lattice range [−N ⊥ , N ⊥ − 1], we take k tot = k q + k g = k Q directly and calculate the transferred momenta as ∆ q ≡ k q − k tot and ∆ g ≡ k g − k tot . An example is shown in Fig. 22.
2. k q + k g = k Q + 2N ⊥ In this situation, the sum k q + k g exceeds the positive boundary of the lattice. This could happen when both k q and k g are large and positive, as the example illustrated in Fig. 23(a). There could be more than one choice in applying the periodic boundary conditions to the momentum quantum numbers. We choose to bring the gluon to the opposite direction as k g → k g − 2N ⊥ . Therefore we calculate the transferred momenta as ∆ q ≡ k q − k Q and ∆ g ≡ k g − (k q + k g ). This prescription is shown in Fig. 23(b). The corresponding physical process is a quark and a gluon carrying large but opposite momenta transforming into/from a quark carrying a small momentum. There could be alternative ways in applying the periodic boundary conditions, as shown in Figs. 23(c) and 23(d). The process shown in Fig. 23(c) is obtained by bringing the quark q one period below, k q → k q − 2N ⊥ . In this interpretation, a quark and a gluon, carrying opposite momentum, transfer into/from a quark carrying a small momentum. Differently, the process shown in Fig. 23(d) is obtained by bringing the quark Q one period above, k Q → k Q +2N ⊥ . In this interpretation, a quark and a gluon, each carrying a positive momentum, transfer into/from a quark carrying a larger positive momentum.
3. k q + k g = k Q − 2N ⊥ This is very similar to the previous situation where k q +k g = k Q +2N ⊥ . This could happen when both k q and k g are large and negative, as the example illustrated in Fig. 24(a). We choose to bring the gluon to the opposite direction as k g → k g + 2N ⊥ [see Fig. 24(b)]. Therefore, we calculate the transferred momenta as ∆ q ≡ k q − k Q and ∆ g ≡ k g − (k q + k g ). Two alternative ways in applying the periodic boundary conditions are shown in Figs. 24(c) and 24(d).
Our choices for all three cases discussed above can be summarized into one as It is generalized to the two-dimensional transverse space in Eq. (27). With this prescription, we could, on the lattice of one fundamental Brillouin zone, maintain the interpretation 24. An example of a quark (denoted as "q") and a gluon (denoted as "g") transferring into/from a quark (denoted as "Q") with their momenta satisfying k q + k g = k Q − 2N ⊥ . The grids inside the fundamental Brillouin zone, i.e., those with momentum numbers in the range of [−N ⊥ , N ⊥ − 1], are in solid lines. The grids outside this range are in dashed lines. In (a), particles are marked at their momentum quantum numbers assigned on the lattice. In (b), particles are marked at their momentum quantum numbers used to calculate the transferred momenta ∆ q and ∆ g . In (c) and (d), two other choices in applying the periodic boundary conditions are shown.
of back-to-back splitting and merging, which is physically the most significant process in the qg ↔ Q transition.
ψ; x + + δx + I = ψ; x + I + Here, the states with the superscript 1 or 2 are the "trial" states that are evaluated at the midpoint and the end point. In the high-energy limit of P + → ∞, V qg,I (x + ) loses its dependence on the light-front time and reduces to V qg . In this case, we could write U RK4 in a collective form. By defining λ ≡ −i/2V qg , the Runge-Kutta algorithm reduces to To see the stability of this method, we can plot |U RK4 (λδx + )| in the complex plane of λδx + . The stability boundary defined by the contour |U RK4 (λδx + )| = 1 is shown in Fig. 25. Note that in this case λ is effectively purely imaginary, and one sees from the plot that the method is very close to unitary for a large range of Im[λδx + ]. In this Appendix, we derive the Wilson line of a quarkgluon state in the eikonal limit and discuss its behavior with regard to the total scattering cross section.
To begin with, consider a quark or a gluon propagating through the background of a classical color field. In the eikonal limit, the momentum of the particle is approximated as P µ = (P + ≈ √ s, P − = 0, P ⊥ = 0) and likewise for the background field P µ A = (P + A = 0, P − A ≈ √ s, P A,⊥ = 0). In such circumstances, the interaction Hamiltonian in the interaction picture is equivalent to that in the Schrödinger picture, V I (x + ) = V(x + ), since the phase factor e ±i/2P − x + reduces to 1. The evolution of the quark interacting with the background field for a finite distance in light-front time, x + = [0, L η ], is written in terms of a fundamental Wilson line, where A µ = a t a A a µ and t a are the SU(3) generators in the fundamental representation. Similarly, the evolution of the gluon in the eikonal limit is described by the adjoint Wilson line, where A µ = a T a A a µ and (T a ) bc = −i f abc are the SU(3) generators in the adjoint representation.
Next, we consider the scattering of a quark-gluon state, in which the quark and the gluon interact with the same background field simultaneously. The scattering amplitude is simply the tensor product of the quark and the gluon Wilson lines, U qg (0, L η ; x ⊥ , y ⊥ ) = U F (0, L η ; x ⊥ ) ⊗ U A (0, L η ; y ⊥ ) . (E3) Physical observables such as the cross section could be determined from the Wilson line averaged over the background field configurations, which is essentially the scattering amplitude. Note that the dimension of the Wilson line is the same as that of the particle's color space. In calculating the total scattering cross section of the particle state l, one should sum over the final color states and average over the initial color states, (E4) where U l is the Wilson line, i, f the color indices, and N l the dimension of the color space of particle l. The trace Tr is over the color indices. The . . . stands for a configuration average of the background field.
(E9) The solution is U R 1 (0, L η ; x ⊥ ) β 1 α 1 U R 2 (0, L η ; y ⊥ ) β 2 α 2 =Ū R 1 (0, L η ; x ⊥ )Ū R 2 (0, L η ; y ⊥ ) For a quark-gluon state scattering on a background field with constantμ, Eq. (E10) becomes, In calculating the total cross section, it is the real part of the trace of the averaged Wilson line that matters as in Eq. (E4), The product of the quark and the gluon Wilson linesŪ qg plotted as a function of the dimensionless quantity g 2μ L η /m g at various m g r, where r = | x ⊥ − y ⊥ |, according to Eq. (E12). The part without the quark-gluon correlation,Ū FŪA = exp −g 4μ2 L η /(8πm 2 g )(C F + C A ) is plotted in the solid line.
so we are interested in the following expression: and it is plotted in Fig. 26. It is a periodic function with a period of 4π and oscillates between 1 and −1. In Eq. (E12), this term depends on the dimensionless quantity g 2μ L η /m g , just as theŪ FŪA term, but it also depends on the separation between the quark and the gluon, r ≡ | x ⊥ − y ⊥ |. The smaller the value of r is, the faster f qg deviates from 1 as a function of g 2μ L η /m g , suggesting that the correlation is stronger when the quark and the gluon are closer. One could also see this in the limit of infinite separation where r = ∞, the correlation becomes f qg (0) = 1.
The contribution from f qg (ξ) as a correction to theŪ FŪA term insideŪ qg is actually very small. Even in the strongest correlation case where m g r = 0, the first node of f qg = 0 occurs at g 2μ L η /m g = 2π, where the value ofŪ FŪA already reduces to 0.0011. We present the plots ofŪ FŪA and the correlated Wilson lineŪ qg as functions of the dimensionless quantity g 2μ L η /m g at various m g r in Fig. 27. ThoseŪ qg curves barely deviate fromŪ FŪA , even in the zero separation case. Indeed, the influence from the correlation function f qg is not very noticeable inŪ qg as we have discussed. From here, one could expect that f qg has little influence in the total cross section of a |qg state interacting with the background field.