Covariant Chiral Kinetic Equation in Non-Abelian Gauge field from “covariant gradient expansion”

In recent years, there has been a considerable amount of theoretical work on the chiral kinetic theory (CKT) in relativistic heavy ion collisions. The CKT aims to incorporate the chiral anomaly into kinetic theory and provide a consistent formalism to describe various novel chiral effects, e.g., chiral magnetic effect Vilenkin:1980fu ; Kharzeev:2007jp ; Fukushima:2008xe , chiral vortical effect Vilenkin:1978hb ; Kharzeev:2007tn ; Erdmenger:2008rm ; Banerjee:2008th , chiral separation effect Son:2004tq ; Metlitski:2005pr and so on, which are all associated with the chiral anomaly. Recent progress on chiral effects and chiral kinetic theory in relativistic heavy ion collisions can be found in the reviews such as Kharzeev:2013ffa ; Kharzeev:2015znc ; Liu:2020ymh ; Gao:2020vbh ; Gao:2020pfu . The chiral kinetic equation has been derived from various methods, such as semiclassical approach Duval:2005vn ; Wong:2011nt ; Son:2012wh ; Stephanov:2012ki ; Stone:2013sga ; Dwivedi:2013dea ; Akamatsu:2014yza ; Chen:2014cla ; Manuel:2014dza ; Hayata:2017ihy , Wigner function formalism Gao:2012ix ; Chen:2012ca ; Hidaka:2016yjf ; Huang:2018wdl ; Gao:2018wmr ; Liu:2018xip , effective field theory Son:2012zy ; Carignano:2018gqt ; Lin:2019ytz ; Carignano:2019zsh and world-line approach Mueller:2017lzw ; Mueller:2017arw ; Mueller:2019gjj . The numerical simulation based on chiral kinetic equation can be found in Refs. Sun:2016nig ; Sun:2016mvh ; Sun:2017xhx ; Sun:2018idn ; Sun:2018bjl ; Zhou:2018rkh ; Zhou:2019jag ; Liu:2019krs .

Despite all these development, so far most of the literature focuses on the CKT in Abelian gauge field. Only very restricted work Stone:2013sga ; Akamatsu:2014yza ; Hayata:2017ihy ; Mueller:2019gjj had discussed the CKT in non-Abelian gauge field. However, as we all know, the dynamics of the produced quark-gluon plasma in relativistic heavy ion collisions are mainly determined by quantum chromodynamics — non-Abelian $SU(3)$ gauge field. Especially, in the small $x$ physics, the initial state in relativistic nucleus-nucleus collisions can be described as a classical coherent non-Abelian gauge field configuration called the color glass condensateGribov:1984tu ; Mueller:1985wy ; McLerran:1993ka ; McLerran:1994vd ; Iancu:2003xm . It still remains an open question how the decoherence from the classical color field to the quark gluon plasma takes place. In order to address these problems, we need generalize the CKT in Abelian gauge field to the one in non-Abelian gauge field.

In this paper, we will be dedicated to deriving the chiral kinetic equation in $SU(N)$ gauge field from the quantum transport theory Heinz:1983nx ; Elze:1986hq ; Elze:1986qd ; Elze:1989un ; Ochs:1998qj based on the Wigner functions from quantum gauge field theory. In Sec.2, we review the Wigner function formalism given in Refs. Heinz:1983nx ; Elze:1986hq ; Elze:1986qd ; Elze:1989un and present some results in Ref.Ochs:1998qj that would be useful for our present work. In Sec.3, we apply the “covariant gradient expansion” given in Elze:1986hq ; Elze:1986qd ; Elze:1989un ; Ochs:1998qj to expanding the Wigner equations for massless fermions up to the first order and disentangle the Wigner equations by the method developed in the Abelian case in Ref.Gao:2018wmr . We find that only the timelike component of the Wigner functions is independent and all other spacelike components can be derivative from timelike component directly. Such result is very similar to the Abelian case and reduces the Wigner equations greatly. We present the covariant chiral kinetic equation for this independent Wigner function in 8-dimensional form, i.e., 4-dimensional momentum space and 4-dimensional coordinate space. In comparison with the Abelian case, the extra constraint equation appears in non-Abelian case. In Sec.4, we decompose the results further in the color space and find that the color singlet phase-space distribution function and multiplet ones are totally coupled with each other. In Sec.5, We discuss the modified Lorentz transformation of the distribution function in phase space when we define it in different reference frames. With the results in previous sections, we calculate the vector and axial currents induced by color field and vorticity in Sec.6. It turns out that the non-Abelian chiral anomaly can be derived directly from the 4-dimentional Berry curvature in the vacuum contribution of the color singlet Wigner function. With specific distribution near global equilibrium, we can obtain the non-Abelian counterparts of chiral magnetic effect and chiral vortical effect. Finally, we summarize the paper in Sec.7.

In this work, we use the convention for the metric $g^{\mu\nu}=\mathrm{diag}(1,-1,-1,-1)$, Levi-Civita tensor $\epsilon^{0123}=1$. We choose natural units such that $\hbar=c=1$ except for the cases when we want to display $\hbar$ dependence to clarify the perturbative expansion.

In quantum transport theory, the gauge invariant density matrix for spin-1/2 quarks is defined as Heinz:1983nx ; Elze:1986hq ; Elze:1986qd

$\displaystyle\rho\left(x+\frac{y}{2},x-\frac{y}{2}\right)=\bar{\psi}\left(x+\frac{y}{2}\right)U\left(x+\frac{y}{2},x\right)\otimes U\left(x,x-\frac{y}{2}\right)\psi\left(x-\frac{y}{2}\right).$

(1)

where the direct product is over both spinor and color indices. The element of density matrix with specific color and spinor indices is given by

$\displaystyle\rho_{\alpha\beta}^{ij}\left(x+\frac{y}{2},x-\frac{y}{2}\right)=\bar{\psi}_{\beta}^{j^{\prime}}\left(x+\frac{y}{2}\right)U^{j^{\prime}j}\left(x+\frac{y}{2},x\right)U^{ii^{\prime}}\left(x,x-\frac{y}{2}\right)\psi^{i^{\prime}}_{\alpha}\left(x-\frac{y}{2}\right).$

(2)

where $\alpha,\beta$ denote spinor indices, $i,i^{\prime},j,j^{\prime}$ mean color indices in fundamental representation and $U^{j^{\prime}j}$ or $U^{i^{\prime}i}$ is the Wilson line or gauge link

$\displaystyle U^{ij}(x,y)=\left[{P}\exp\left(\frac{ig}{\hbar}\int_{y}^{x}dz^{\mu}A_{\mu}(z)\right)\right]^{ij}$

(3)

which is necessary to keep the operator gauge invariant. In the definition of Wilson line $P$ denotes path ordering of the operator and the integral in the exponent is taken along the straight path from $x$ to $y$. The gauge field potential is defined by $A_{\mu}=A^{a}_{\mu}t^{a}$, with the $N^{2}-1$ hermitian generators of $SU(N)$ in the fundamental representation satisfying

$\displaystyle\textrm{Tr}\,t^{a}=0,\ \ \ \left[t^{a},t^{b}\right]=if^{abc}t^{c},\ \ \ \left\{t^{a},t^{b}\right\}=\frac{1}{N}\delta^{ab}{\bf 1}+d^{abc}t^{c}.$

(4)

For non-Abelian gauge field, the covariant derivative in the fundamental representation is defined as,

$\displaystyle D_{\mu}(x)$

$\displaystyle=$

$\displaystyle\partial_{\mu}-\frac{ig}{\hbar}A_{\mu}(x),$

(5)

and the field strength tensor follows as

$\displaystyle F_{\mu\nu}(x)$

$\displaystyle\equiv$

$\displaystyle F_{\mu\nu}^{a}t^{a}=-\frac{\hbar}{ig}\left[D_{\mu},D_{\nu}\right]=\partial_{\mu}A_{\nu}(x)-\partial_{\nu}A_{\mu}(x)-\frac{ig}{\hbar}\left[A_{\mu}(x),A_{\nu}(x)\right].$

(6)

The Wigner operator $\hat{W}(x,p)$ is related to the gauge invariant density matrix by Fourier transformation

$\displaystyle\hat{W}(x,p)=\int\frac{d^{4}y}{(2\pi)^{4}}e^{-ip\cdot y}\rho\left(x+\frac{y}{2},x-\frac{y}{2}\right),$

(7)

and the Wigner function is defined as ensemble averaging of the Wigner operator

$\displaystyle W(x,p)=\langle\hat{W}(x,p)\rangle.$

(8)

In our present work, we will concentrate on the quark matter under a purely classical external non-Abelian gauge field, in which ordinary matrix multiplication rules in spinor space or color space suffice and the Wigner equations will not generate the so-called BBGKY-hierarchy Groot:1980 and can be closed by itself

			$\displaystyle\left[m-\gamma^{\mu}\left(p_{\mu}+\frac{1}{2}i\mathscr{D}_{\mu}(x)\right)\right]W(x,p)$		(9)
		$\displaystyle=$	$\displaystyle\frac{ig}{2}\gamma^{\mu}\partial^{\nu}_{p}\left\{\int_{0}^{1}ds{\,}\frac{1+s}{2}\left[e^{-\frac{1}{2}is\Delta}F_{\mu\nu}(x)\right]W(x,p)\right.$
			$\displaystyle\left.\hskip 42.67912pt+W(x,p)\int_{0}^{1}ds{\,}\frac{1-s}{2}\left[e^{\frac{1}{2}is\Delta}F_{\mu\nu}(x)\right]\right\},$

together with the hermitian adjoint equation

			$\displaystyle W(x,p)\left[m-\gamma^{\mu}\left(p_{\mu}-\frac{1}{2}i\mathscr{D}^{\dagger}_{\mu}(x)\right)\right]$		(10)
		$\displaystyle=$	$\displaystyle-\frac{ig}{2}\partial^{\nu}_{p}\left\{\int_{0}^{1}ds{\,}\frac{1-s}{2}\left[e^{-\frac{1}{2}is\Delta}F_{\mu\nu}(x)\right]W(x,p)\right.$
			$\displaystyle\left.\hskip 42.67912pt+W(x,p)\int_{0}^{1}ds{\,}\frac{1+s}{2}\left[e^{\frac{1}{2}is\Delta}F_{\mu\nu}(x)\right]\right\}\gamma^{\mu},$

where we have introduced the definition of covariant derivative in the adjoint representation for a second-rank tensor $\mathcal{T}(x)$ in color space by

$\displaystyle\mathscr{D}_{\mu}(x)\mathcal{T}(x)\equiv\left[{D}_{\mu}(x),\mathcal{T}(x)\right]=\partial_{\mu}^{x}\mathcal{T}(x)-\frac{ig}{\hbar}\left[A_{\mu}(x),\mathcal{T}(x)\right],$

(11)

and $\Delta\equiv\partial_{p}\cdot\mathscr{D}(x)$ with $\mathscr{D}(x)$ only acting on $F_{\mu\nu}$ and $\partial_{p}$ always on $W$ after or in front of it. It should be noted that in the definition of the Wigner function given by Eq. (8) and the Wigner equations (9) and (10) there is no normal ordering in the Wigner matrix because we did not make any manipulation on the order of the quark field. It has been demonstrated in Gao:2019zhk ; Fang:2020com that this plays a central role to give rise to the chiral anomaly in the quantum kinetic theory.

If we take the convention in Ochs:1998qj , momentum derivatives standing to the right of the Wigner function are defined in the sense of partial integration as

$\displaystyle W(x,p)\partial^{\nu_{1}}_{p}\cdots\partial^{\nu_{k}}_{p}\equiv(-1)^{k}\partial^{\nu_{k}}_{p}\cdots\partial^{\nu_{1}}_{p}W(x,p),$

(12)

and define generalized non-local momentum and derivative operators $\Pi_{\mu}$ and $G_{\mu}$ as

	$\displaystyle\Pi_{\mu}$	$\displaystyle=$	$\displaystyle p_{\mu}+\frac{g}{2}\int_{0}^{1}ds{\,}\left(e^{-\frac{1}{2}is\Delta}F_{\mu\nu}(x)\right)is\partial^{\nu}_{p},$
	$\displaystyle G_{\mu}$	$\displaystyle=$	$\displaystyle D_{\mu}+\frac{g}{2}\int_{0}^{1}ds{\,}\left(e^{-\frac{1}{2}is\Delta}F_{\mu\nu}(x)\right)\partial^{\nu}_{p},$		(13)

the Wigner equations can be cast into a more compact form Ochs:1998qj ,

	$\displaystyle 2mW(x,p)$	$\displaystyle=$	$\displaystyle\gamma^{\mu}\left(\left\{\Pi_{\mu},W(x,p)\right\}+i\left[G_{\mu},W(x,p)\right]\right),$		(14)
	$\displaystyle 2mW(x,p)$	$\displaystyle=$	$\displaystyle\left(\left\{\Pi_{\mu},W(x,p)\right\}-i\left[G_{\mu},W(x,p)\right]\right)\gamma^{\mu}.$		(15)

Adding or subtracting the two equations above gives

	$\displaystyle 4mW(x,p)$	$\displaystyle=$	$\displaystyle\left\{\gamma^{\mu},\left\{\Pi_{\mu},W(x,p)\right\}\right\}+i\left[\gamma^{\mu},\left[G_{\mu},W(x,p)\right]\right],$		(16)
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\left[\gamma^{\mu},\left\{\Pi_{\mu},W(x,p)\right\}\right]+i\left\{\gamma^{\mu},\left[G_{\mu},W(x,p)\right]\right\}.$		(17)

In spinor space, we can decompose the Wigner function into

$\displaystyle W=\frac{1}{4}\left[{{\mathscr{F}}}+i\gamma^{5}{{\mathscr{P}}}+\gamma^{\mu}{{\mathscr{V}}}_{\mu}+\gamma^{\mu}\gamma^{5}{{\mathscr{A}}}_{\mu}+\frac{1}{2}\sigma^{\mu\nu}{{\mathscr{S}}}_{\mu\nu}\right].$

(18)

In this work, we will restrict ourselves to the massless or chiral fermions. In consequence, if we introduce a chirality basis via

$\displaystyle{\mathscr{J}}^{\mu}_{s}$

$\displaystyle=$

$\displaystyle\frac{1}{2}\left({\mathscr{V}}^{\mu}+s{\mathscr{A}}^{\mu}\right),$

(19)

where $s=+1/-1$ denotes right-handed/left-handed component, the equations for the chiral Wigner function ${\mathscr{J}}^{\mu}_{s}$ will decouple from all the other components of the Wigner function and each other as well, which leads to

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\left\{\Pi_{\mu},{\mathscr{J}}^{\mu}_{s}\right\},$		(20)
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\left[G_{\mu},{\mathscr{J}}^{\mu}_{s}\right],$		(21)
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\left\{\Pi^{\mu},{\mathscr{J}}^{\nu}_{s}\right\}-\left\{\Pi^{\nu},{\mathscr{J}}^{\mu}_{s}\right\}{+}s\hbar\epsilon^{\mu\nu\alpha\beta}\left[G_{\alpha},{\mathscr{J}}_{s\beta}\right],$		(22)

where we have recovered the $\hbar$ dependence before the generalized derivative operators in the last equation in order to make perturbative expansion in the following section. These Wigner equations will be the starting point of our present work in the following. For brevity, we will suppress the subscript $s$ of the left-hand or right-hand Wigner function ${\mathscr{J}}^{\mu}_{s}$ in the subsequent sections and reinstate it when it is necessary.

In the Abelian plasma, the disentanglement theorem of Wigner functions has been demonstrated in Ref. Gao:2018wmr , which tell us that up to any order of $\hbar$ among four components of Wigner functions $\mathscr{J}^{\mu}$ only the timelike component is independent and satisfies only one independent Wigner equation, the other spatial components can be totally fixed from this independent Wigner function and the Wigner equations for them are all satisfied automatically. Now let us try to generalize this disentanglement formalism from Abelian gauge field to non-Abelian gauge field. In order to achieve this goal, we will resort to the “covariant gradient expansion” proposed in Refs. Elze:1986qd ; Elze:1989un ; Ochs:1998qj . In this expansion scheme, when we have one extra covariant derivative $D_{\mu}$ or $\mathscr{D}_{\mu}$, we will have one extra higher order contribution. The “covariant gradient expansion” preserves gauge invariance order by order automatically. Actually we can trace such expansion in powers of $\hbar$, e.g., in the Wigner equations (22) and the generalized non-local momentum and derivative operators

	$\displaystyle\Pi_{\mu}$	$\displaystyle=$	$\displaystyle\sum_{k=0}^{\infty}\hbar^{k}\Pi_{\mu}^{(k)}=p_{\mu}-\hbar\frac{ig}{2}\sum_{k=0}^{\infty}\left(-\frac{i\hbar}{2}\right)^{k}\frac{k+1}{(k+2)!}\left[\left(\partial_{p}\cdot{\mathscr{D}}\right)^{k}F_{\nu\mu}\right]\partial_{p}^{\nu}$		(23)
	$\displaystyle G_{\mu}$	$\displaystyle=$	$\displaystyle\sum_{k=0}^{\infty}\hbar^{k}G_{\mu}^{(k)}=D_{\mu}-\frac{g}{2}\sum_{k=0}^{\infty}\left(-\frac{i\hbar}{2}\right)^{k}\frac{1}{(k+1)!}\left[\left(\partial_{p}\cdot{\mathscr{D}}\right)^{k}F_{\nu\mu}\right]\partial_{p}^{\nu}.$		(24)

Up to the second order of $\hbar$, the non-local operators $\Pi_{\mu}$ and $G_{\mu}$ are given by

	$\displaystyle\Pi_{\mu}^{(0)}$	$\displaystyle=$	$\displaystyle p_{\mu},\ \ \ \Pi_{\mu}^{(1)}=\frac{ig}{4}F_{\mu\nu}\partial^{\nu}_{p},\ \ \Pi_{\mu}^{(2)}=\frac{g}{12}\left[\left(\partial_{p}\cdot{\mathscr{D}}\right)F_{\mu\nu}\right]\partial_{p}^{\nu},$		(25)
	$\displaystyle G_{\mu}^{(0)}$	$\displaystyle=$	$\displaystyle D_{\mu}+\frac{g}{2}F_{\mu\nu}\partial^{\nu}_{p},\ \ G_{\mu}^{(1)}=-\frac{ig}{8}\left[\left(\partial_{p}\cdot{\mathscr{D}}\right)F_{\mu\nu}\right]\partial_{p}^{\nu}.$		(26)

We can also expand the Wigner operator as

$\displaystyle W(x,p)=\sum_{k=0}^{\infty}\hbar^{k}W^{(k)}(x,p).$

(27)

However it should be noted that the “covariant gradient expansion” is not completely identical to an expansion in powers of $\hbar$ for non-Abelian gauge field which had been pointed out in Elze:1986qd ; Elze:1989un ; Ochs:1998qj though it is identical for Abelian gauge field. In non-Abelian case, there is an extra gauge potential $A_{\mu}$ with $ig/\hbar$ in the covariant derivative $D_{\mu}$ or $\mathscr{D}_{\mu}$ in Eqs. (23) and (24) while there only exist ordinary derivative $\partial_{\mu}^{x}$ in the Abelian case.

In order to disentangle the Wigner equations further, it is convenient to introduce time-like 4-vector $n^{\mu}$ with normalization $n^{2}=1$. For simplicity we assume $n^{\mu}$ is a constant vector. With the auxiliary vector $n^{\mu}$, we can decompose any vector $X^{\mu}$ into the component parallel to $n^{\mu}$ and the other components perpendicular to $n^{\mu}$,

$$X^{\mu}=X_{n}n^{\mu}+\bar{X}^{\mu},$$

(28)

where $X_{n}=X\cdot n$ and $\bar{X}^{\mu}=\Delta^{\mu\nu}X_{\nu}$ with $\Delta^{\mu\nu}=g^{\mu\nu}-n^{\mu}n^{\nu}$. The gauge field tensor $F^{\mu\nu}$ can be also decomposed into

$$F^{\mu\nu}=E^{\mu}n^{\nu}-E^{\nu}n^{\mu}-\bar{\epsilon}^{\mu\nu\sigma}B_{\sigma}$$

(29)

with

$$E^{\mu}=F^{\mu\nu}n_{\nu},\ \ B^{\mu}=\frac{1}{2}\bar{\epsilon}^{\mu\rho\sigma}F_{\rho\sigma},$$

(30)

where for notational convenience we have defined $\bar{\epsilon}^{\mu\alpha\beta}=\epsilon^{\mu\nu\alpha\beta}n_{\nu}$.

Now we can decompose the Wigner functions and Wigner equations along the direction $n^{\mu}$ order by order. The leading order or the zeroth order result is very simple

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle p_{n}{\mathscr{J}}_{n}^{(0)}+\bar{p}_{\mu}\bar{{\mathscr{J}}}^{(0)\mu},$		(31)
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\left[G_{n}^{(0)},{\mathscr{J}}_{n}^{(0)}\right]+\left[\bar{G}_{\mu}^{(0)},\bar{{\mathscr{J}}}^{(0)\mu}\right],$		(32)
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\bar{p}^{\mu}{\mathscr{J}}_{n}^{(0)}-p_{n}\bar{{\mathscr{J}}}^{(0)\mu},$		(33)
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\bar{p}^{\mu}\bar{{\mathscr{J}}}^{(0)\nu}-\bar{p}^{\nu}\bar{{\mathscr{J}}}^{(0)\mu}.$		(34)

From Eq.(33), we can express the space-like component $\bar{{\mathscr{J}}}^{(0)\mu}$ in terms of ${\mathscr{J}}_{n}^{(0)}$

$\displaystyle\bar{{\mathscr{J}}}^{(0)\mu}$

$\displaystyle=$

$\displaystyle\bar{p}^{\mu}\frac{{\mathscr{J}}_{n}^{(0)}}{p_{n}}.$

(35)

Substituting this relation into Eqs.(31) gives rise to the on-shell condition

$\displaystyle p^{2}\frac{{\mathscr{J}}_{n}^{(0)}}{p_{n}}$

$\displaystyle=$

$\displaystyle 0,$

(36)

which means ${{\mathscr{J}}_{n}^{(0)}}/{p_{n}}$ must be proportional to the Dirac delta function $\delta(p^{2})$

$\displaystyle\frac{{\mathscr{J}}_{n}^{(0)}}{p_{n}}$

$\displaystyle=$

$\displaystyle f^{(0)}\delta(p^{2}),$

(37)

where $f^{(0)}$ can be regarded as the usual particle distribution function in four-dimensional momentum space and four-dimensional coordinate space. It must be non-singular function at $p^{2}=0$. Putting Eqs.(37) and (35) together, we get the full Wigner function of the zeroth order

$\displaystyle{\mathscr{J}}^{(0)\mu}$

$\displaystyle=$

$\displaystyle p^{\mu}f^{(0)}\delta(p^{2}).$

(38)

The transport equation satisfied by $f^{(0)}$ can be obtained from Eq.(32)

$\displaystyle 0$

$\displaystyle=$

$\displaystyle\left[G_{\mu}^{(0)},p^{\mu}f^{(0)}\delta(p^{2})\right].$

(39)

It is obvious that Eq.(34) is automatically satisfied with the expression (35).

The next-to-leading order or the first order equations are given by

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle 2p_{n}{\mathscr{J}}^{(1)}_{n}+2\bar{p}_{\mu}\bar{{\mathscr{J}}}^{(1)\mu}+\left\{\Pi_{n}^{(1)},{\mathscr{J}}^{(0)}_{n}\right\}+\left\{\bar{\Pi}_{\mu}^{(1)},\bar{{\mathscr{J}}}^{(0)\mu}\right\},$		(40)
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\left[G_{n}^{(0)},{\mathscr{J}}^{{(1)}}_{n}\right]+\left[\bar{G}_{\mu}^{(0)},\bar{{\mathscr{J}}}^{{(1)}\mu}\right]+\left[G_{n}^{(1)},{\mathscr{J}}^{{(0)}}_{n}\right]+\left[\bar{G}_{\mu}^{(1)},\bar{{\mathscr{J}}}^{{(0)}\mu}\right],\hskip 8.0pt$		(41)
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle 2\bar{p}^{\mu}{\mathscr{J}}^{(1)}_{n}-2p_{n}\bar{{\mathscr{J}}}^{(1)\mu}+\left\{\bar{\Pi}^{(1)\mu},{\mathscr{J}}^{(0)}_{n}\right\}-\left\{\Pi^{(1)}_{n},\bar{{\mathscr{J}}}^{(0)\mu}\right\}$		(42)
			$\displaystyle{+}s\bar{\epsilon}^{\mu\alpha\beta}\left[G_{\alpha}^{(0)},{\mathscr{J}}_{\beta}^{(0)}\right],$
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle 2\bar{p}^{\mu}\bar{{\mathscr{J}}}^{(1)\nu}-2\bar{p}^{\nu}\bar{{\mathscr{J}}}^{(1)\mu}+\left\{\bar{\Pi}^{(1)\mu},\bar{{\mathscr{J}}}^{(0)\nu}\right\}-\left\{\bar{\Pi}^{(1)\nu},\bar{{\mathscr{J}}}^{(0)\mu}\right\}$		(43)
			$\displaystyle{+}s\bar{\epsilon}^{\mu\nu\alpha}\left(\left[G_{\alpha}^{(0)},{\mathscr{J}}_{n}^{(0)}\right]-\left[G_{n}^{(0)},{\mathscr{J}}_{\alpha}^{(0)}\right]\right).$

From Eq.(42), we can express $\bar{{\mathscr{J}}}^{(1)\mu}$ in terms of ${\mathscr{J}}^{(1)}_{n}$ and ${\mathscr{J}}^{(0)}_{n}$

	$\displaystyle\bar{{\mathscr{J}}}^{(1)\mu}$	$\displaystyle=$	$\displaystyle\bar{p}^{\mu}\frac{{\mathscr{J}}^{(1)}_{n}}{p_{n}}{+}\frac{s}{2p_{n}}\bar{\epsilon}^{\mu\alpha\beta}\left[G_{\alpha}^{(0)},\bar{p}_{\beta}\frac{{\mathscr{J}}_{n}^{(0)}}{p_{n}}\right]$		(44)
			$\displaystyle+\frac{1}{2p_{n}}\left(\left\{\bar{\Pi}^{(1)\mu},p_{n}\frac{{\mathscr{J}}^{(0)}_{n}}{p_{n}}\right\}-\left\{\Pi^{(1)}_{n},\bar{p}^{\mu}\frac{{\mathscr{J}}_{n}^{(0)}}{p_{n}}\right\}\right).$

Substituting it into Eqs.(40) and (41) gives rise to the modified on-shell condition and transport equation for ${\mathscr{J}}_{n}^{(1)}$, respectively,

	$\displaystyle p^{2}\frac{{\mathscr{J}}_{n}^{(1)}}{p_{n}}$	$\displaystyle=$	$\displaystyle{-}\frac{s}{2p_{n}}\bar{\epsilon}^{\mu\alpha\beta}\bar{p}_{\mu}\left[G_{\alpha}^{(0)},\bar{p}_{\beta}\frac{{\mathscr{J}}_{n}^{(0)}}{p_{n}}\right]-{\frac{1}{2}}\left\{\Pi_{\mu}^{(1)},p^{\mu}\frac{{\mathscr{J}}^{(0)}_{n}}{p_{n}}\right\}$		(45)
			$\displaystyle-\frac{\bar{p}_{\mu}}{2p_{n}}\left(\left\{\bar{\Pi}^{(1)\mu},p_{n}\frac{{\mathscr{J}}^{(0)}_{n}}{p_{n}}\right\}-\left\{\Pi^{(1)}_{n},\bar{p}^{\mu}\frac{{\mathscr{J}}_{n}^{(0)}}{p_{n}}\right\}\right),$
	$\displaystyle\left[G_{\mu}^{(0)},p^{\mu}\frac{{\mathscr{J}}^{(1)}_{n}}{p_{n}}\right]$	$\displaystyle=$	$\displaystyle{-}\frac{s}{2}\bar{\epsilon}^{\mu\alpha\beta}\left[\bar{G}_{\mu}^{(0)},\frac{1}{p_{n}}\left[G_{\alpha}^{(0)},\bar{p}_{\beta}\frac{{\mathscr{J}}_{n}^{(0)}}{p_{n}}\right]\right]-\left[G_{\mu}^{(1)},p^{\mu}\frac{{\mathscr{J}}^{{(0)}}_{n}}{p_{n}}\right]$		(46)
			$\displaystyle-\frac{1}{2}\left[\bar{G}_{\mu}^{(0)},\frac{1}{p_{n}}\left(\left\{\bar{\Pi}^{(1)\mu},p_{n}\frac{{\mathscr{J}}^{(0)}_{n}}{p_{n}}\right\}-\left\{\Pi^{(1)}_{n},\bar{p}^{\mu}\frac{{\mathscr{J}}_{n}^{(0)}}{p_{n}}\right\}\right)\right].$

It is easy to verify that the general expression of the constraint equation (45) is given by

$\displaystyle\frac{{\mathscr{J}}_{n}^{(1)}}{p_{n}}=f^{(1)}\delta(p^{2}){+}\frac{s}{2p_{n}}\bar{\epsilon}^{\mu\alpha\beta}p_{\mu}\left\{\frac{g}{2}F_{\alpha\beta},f^{(0)}\right\}\delta^{\prime}(p^{2})+\left\{\Pi_{\mu}^{(1)},p^{\mu}f^{(0)}\right\}\delta^{\prime}(p^{2}).$

(47)

Just like $f^{(0)}$, the function $f^{(1)}$ is also a non-singular distribution function at $p^{2}=0$ in four-dimensional momentum space and four-dimensional coordinate space and can be regarded as the first order correction to $f^{(0)}$. The transport equation for $f^{(1)}$ can be directly obtained by inserting Eq.(47) into Eq.(46) and will not be presented explicitly here to avoid too lengthy equations. Putting Eqs.(44) and (47) together, we get the full Wigner function of the first order

	$\displaystyle{\mathscr{J}}^{(1)\mu}$	$\displaystyle=$	$\displaystyle p^{\mu}\left[f^{(1)}\delta(p^{2}){+}\frac{s}{2p_{n}}\bar{\epsilon}^{\nu\alpha\beta}p_{\nu}\left\{\frac{g}{2}F_{\alpha\beta},f^{(0)}\right\}\delta^{\prime}(p^{2})+\left\{\Pi_{\nu}^{(1)},p^{\nu}f^{(0)}\right\}\delta^{\prime}(p^{2})\right]$		(48)
			$\displaystyle+\frac{1}{2p_{n}}\left(\left\{\bar{\Pi}^{(1)\mu},p_{n}f^{(0)}\delta(p^{2})\right\}-\left\{\Pi^{(1)}_{n},\bar{p}^{\mu}f^{(0)}\delta(p^{2})\right\}\right)$
			$\displaystyle{+}\frac{s}{2p_{n}}\bar{\epsilon}^{\mu\alpha\beta}\left[G_{\alpha}^{(0)},\bar{p}_{\beta}f^{(0)}\delta(p^{2})\right].$

As we note in the zeroth order case, the equation (34) is automatically satisfied once we have the expression (35). Now we can check if the first order equation (43) also holds automatically by using the first order expression (44) together with Eqs.(32), (36) and (38). In consequence, after direct calculation we find that the first order equation (43) is not satisfied automatically but lead to the constraint equation for ${\mathscr{J}}_{\mu}^{(0)}$ or $f^{(0)}$

$\displaystyle 0$

$\displaystyle=$

$\displaystyle n_{\alpha}\left(\left[F^{\nu\alpha},{\mathscr{J}}^{(0)\mu}\right]+\left[F^{\alpha\mu},{\mathscr{J}}^{(0)\nu}\right]+\left[F^{\mu\nu},{\mathscr{J}}^{(0)\alpha}\right]\right).$

(49)

Because $n^{\alpha}$ is an arbitrary auxiliary vector with normalization $n^{2}=1$, the constraint equation should not depend on $n^{\alpha}$ or this equation should hold for any $n^{\alpha}$. This leads to the Lorentz covariant constraint equation

$\displaystyle\left[F^{\nu\alpha},{\mathscr{J}}^{(0)\mu}\right]+\left[F^{\alpha\mu},{\mathscr{J}}^{(0)\nu}\right]+\left[F^{\mu\nu},{\mathscr{J}}^{(0)\alpha}\right]=0,$

(50)

which is equivalent to

$\displaystyle\left[\tilde{F}^{\alpha\beta},{\mathscr{J}}^{(0)}_{\alpha}\right]=0\ \ \ {\textrm{with}\ \ \ \tilde{F}^{\alpha\beta}=\frac{1}{2}\epsilon^{\alpha\beta\mu\nu}}F_{\mu\nu}.$

(51)

In Ref. Ochs:1998qj , similar constraints for $\mathscr{F}$ and $\mathscr{S}_{\mu\nu}$ in Eq.(18) had already been obtained. Such constraints only arise in the quantum transport theory with non-Abelian gauge field. The disentanglement theorem of Wigner functions in Abelian gauge field given in Ref. Gao:2018wmr show that all these constraint equations in Abelian cases are satisfied automatically and holds up to any order of $\hbar$. We also notice that the first order equation (43) gives the constraint for the zeroth order Wigner function ${\mathscr{J}}^{(0)\mu}$ because the first order Wigner functions are totally canceled due to the antisymmetry of the equation. Hence in order to get the constraint for the first order Wigner function ${\mathscr{J}}^{(1)\mu}$, we need the second order Wigner functions and equations. The second order expression of Eq.(22) is given by

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle 2\bar{p}^{\mu}{\mathscr{J}}^{(2)}_{n}-2p_{n}\bar{{\mathscr{J}}}^{(2)\mu}$		(52)
			$\displaystyle+\left\{\bar{\Pi}^{(1)\mu},{\mathscr{J}}^{(1)}_{n}\right\}-\left\{\Pi^{(1)}_{n},\bar{{\mathscr{J}}}^{(1)\mu}\right\}+\left\{\bar{\Pi}^{(2)\mu},{\mathscr{J}}^{(0)}_{n}\right\}-\left\{\Pi^{(2)}_{n},\bar{{\mathscr{J}}}^{(0)\mu}\right\}$
			$\displaystyle{+}s\bar{\epsilon}^{\mu\alpha\beta}\left(\left[G_{\alpha}^{(0)},{\mathscr{J}}_{\beta}^{(1)}\right]+\left[G_{\alpha}^{(1)},{\mathscr{J}}_{\beta}^{(0)}\right]\right),$
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle 2\bar{p}^{\mu}\bar{{\mathscr{J}}}^{(2)\nu}-2\bar{p}^{\nu}\bar{{\mathscr{J}}}^{(2)\mu}$		(53)
			$\displaystyle+\left\{\bar{\Pi}^{(1)\mu},\bar{{\mathscr{J}}}^{(1)\nu}\right\}-\left\{\bar{\Pi}^{(1)\nu},\bar{{\mathscr{J}}}^{(1)\mu}\right\}+\left\{\bar{\Pi}^{(2)\mu},\bar{{\mathscr{J}}}^{(0)\nu}\right\}-\left\{\bar{\Pi}^{(2)\nu},\bar{{\mathscr{J}}}^{(0)\mu}\right\}$
			$\displaystyle{+}s\bar{\epsilon}^{\mu\nu\alpha}\left(\left[G_{\alpha}^{(0)},{\mathscr{J}}_{n}^{(1)}\right]-\left[G_{n}^{(0)},{\mathscr{J}}_{\alpha}^{(1)}\right]+\left[G_{\alpha}^{(1)},{\mathscr{J}}_{n}^{(0)}\right]-\left[G_{n}^{(1)},{\mathscr{J}}_{\alpha}^{(0)}\right]\right).$

From the first equation above, we can express $\bar{{\mathscr{J}}}^{(2)\mu}$ in terms of ${\mathscr{J}}^{(2)}_{n}$, ${\mathscr{J}}^{(1)}_{n}$ and ${\mathscr{J}}^{(0)}_{n}$ as

	$\displaystyle\bar{{\mathscr{J}}}^{(2)\mu}$	$\displaystyle=$	$\displaystyle\bar{p}^{\mu}\frac{{\mathscr{J}}^{(2)}_{n}}{p_{n}}{+}\frac{s}{2p_{n}}\bar{\epsilon}^{\mu\alpha\beta}\left(\left[G_{\alpha}^{(0)},{\mathscr{J}}_{\beta}^{(1)}\right]+\left[G_{\alpha}^{(1)},{\mathscr{J}}_{\beta}^{(0)}\right]\right)$		(54)
			$\displaystyle+\frac{1}{2p_{n}}\left(\left\{\bar{\Pi}^{(1)\mu},{\mathscr{J}}^{(1)}_{n}\right\}-\left\{\Pi^{(1)}_{n},\bar{{\mathscr{J}}}^{(1)\mu}\right\}\right)$
			$\displaystyle+\frac{1}{2p_{n}}\left(\left\{\bar{\Pi}^{(2)\mu},{\mathscr{J}}^{(0)}_{n}\right\}-\left\{\Pi^{(2)}_{n},\bar{{\mathscr{J}}}^{(0)\mu}\right\}\right).$

Similar to the first order, substituting it into Eq.(53) and using Eqs. (44), (45) and (46) leads to the constraint for ${\mathscr{J}}^{(1)\mu}$

$\displaystyle\left[\tilde{F}_{\alpha\beta},{\mathscr{J}}^{(1)\alpha}\right]$

$\displaystyle=$

$\displaystyle-\frac{3}{32}\left[\left[\tilde{F}_{\nu\alpha}\partial^{p}_{\beta}-\tilde{F}_{\nu\beta}\partial^{p}_{\alpha},F^{\nu\kappa}\partial^{p}_{\kappa}\right],{\mathscr{J}}^{(0)\alpha}\right].$

(55)

As we just mentioned above, these constraints are unique for non-Abelian gauge field and absent for Abelian field. Such constraints actually originate from the fact that the “covariant gradient expansion” is not completely identical to an expansion in powers of $\hbar$ for non-Abelian gauge field. One difference between non-Abelian and Abelian is the operator $G_{\mu}^{(0)}$. In the non-Abelian case, the derivative in $G_{\mu}^{(0)}$ is covariant derivative $D_{\mu}$, while in Abelian case, it is ordinary space-time derivative $\partial^{x}_{\mu}$. When we calculate high order contribution through iterative process, we will meet the commutator $[D_{\mu},D_{\nu}]=igF_{\mu\nu}/\hbar$ in non-Abelian gauge field and this term will contribute to the lower power order, but for the ordinary derivative such issue will never happen in Abelian gauge field. Actually, during our calculation of (55), we find that if we do not use the constraints for ${\mathscr{J}}^{(0)\alpha}$ in Eq. (50) or (51) beforehand, we will have the same term as the right side of Eq. (49) but with minus sign. This term from the second order equation will eventually cancel the one from the first order. Although we can not give the general proof, we expect that the third order equation of (22) will cancel the second order result (55) and so on. Adding all the contributions up to any high order, the constraint equation (22) should also be satisfied automaticaly.

Up to now, the Wigner function ${{\mathscr{J}}}^{\mu}$ is still an $N\times N$ matrix in color space. Hence it is necessary to decompose the Wigner function into color singlet and multiplet components:

$\displaystyle{{\mathscr{J}}}_{\mu}(x,p)$

$\displaystyle=$

$\displaystyle{\mathscr{J}}_{\mu}^{I}(x,p){\bf{1}}+{{\mathscr{J}}}^{a}_{\mu}(x,p)t^{a},$

(56)

with

$\displaystyle{\mathscr{J}}_{\mu}^{I}(x,p)$

$\displaystyle=$

$\displaystyle\frac{1}{N}\textrm{tr}{{\mathscr{J}}}_{\mu}(x,p),\,\ \ \ {{\mathscr{J}}}^{a}_{\mu}(x,p)={2}\textrm{tr}\,t^{a}{{\mathscr{J}}}_{\mu}(x,p).$

(57)

It should be noted that we use upper index “$I$” to denote singlet component. Similarly, we can decompose the operators into the color singlet and multiplet contributions:

$\displaystyle G_{\mu}^{(0)}$

$\displaystyle=$

$\displaystyle D_{\mu}+G_{\mu}^{(0)a}t^{a},\ \ \ \Pi_{\mu}^{(1)}=\Pi_{\mu}^{(1)a}t^{a},\ \ \ G_{\mu}^{(1)}=G_{\mu}^{(1)a}t^{a},$

(58)

where

$\displaystyle G_{\mu}^{(0)a}=\frac{g}{2}F_{\mu\nu}^{a}\partial_{p}^{\nu},\ \ \ \Pi_{\mu}^{(1)a}=\frac{ig}{4}F_{\mu\nu}^{a}\partial_{p}^{\nu},\ \ \ G_{\mu}^{(1)a}=-\frac{ig}{8}\left(\mathscr{D}^{ac}_{\lambda}F^{c}_{\mu\nu}\right)\partial_{p}^{\lambda}\partial_{p}^{\nu},$

(59)

with ${\mathscr{D}}^{ac}_{\lambda}=\delta^{ca}\partial^{x}_{\lambda}+gf^{bca}A^{b}_{\lambda}/\hbar$. With such decomposition, the singlet and multiplet components of Wigner functions at the zeroth order can be derived from Eqs.(38)

$\displaystyle{\mathscr{J}}^{(0){I}\mu}$

$\displaystyle=$

$\displaystyle p^{\mu}f^{(0){I}}\delta(p^{2}),\ \ \ {\mathscr{J}}^{(0)a\mu}=p^{\mu}f^{(0)a}\delta(p^{2}),$

(60)

which satisfy the coupled transport equations

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\partial^{x}_{\mu}{\mathscr{J}}^{(0){I}\mu}+\frac{1}{N}G_{\mu}^{(0)a}{\mathscr{J}}^{(0)a\mu},$		(61)
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle{\mathscr{D}}^{ac}_{\mu}{\mathscr{J}}^{(0)c\mu}+2G_{\mu}^{(0)a}{\mathscr{J}}^{(0){I}\mu}+d^{bca}G_{\mu}^{(0)b}{\mathscr{J}}^{(0)c\mu}.$		(62)

Similarly but more complicatedly, the color decomposition of first order Wigner functions can be derived from Eq (48)

	$\displaystyle{\mathscr{J}}^{(1){I}\mu}$	$\displaystyle=$	$\displaystyle p^{\mu}f^{(1){I}}\delta(p^{2}){-}\frac{s}{2}\epsilon^{\mu\nu\alpha\beta}p_{\nu}\frac{g}{2N}F^{a}_{\alpha\beta}f^{(0)a}\delta^{\prime}(p^{2})$		(63)
			$\displaystyle{+}\frac{s}{2p_{n}}\bar{\epsilon}^{\mu\alpha\beta}p_{\beta}\left(\partial^{x}_{\alpha}f^{(0){I}}+\frac{1}{N}G_{\alpha}^{(0)a}f^{(0)a}\right)\delta(p^{2}),$
	$\displaystyle{\mathscr{J}}^{(1)a\mu}$	$\displaystyle=$	$\displaystyle p^{\mu}f^{(1)a}\delta(p^{2}){-}s\epsilon^{\mu\nu\alpha\beta}p_{\nu}\left(\frac{g}{2}F^{a}_{\alpha\beta}f^{(0){I}}+\frac{1}{2}d^{bca}\frac{g}{2}F^{b}_{\alpha\beta}f^{(0)c}\right)\delta^{\prime}(p^{2})$		(64)
			$\displaystyle{+}\frac{s}{2p_{n}}\bar{\epsilon}^{\mu\alpha\beta}p_{\beta}\left(\frac{}{}{\mathscr{D}}^{ac}_{\alpha}f^{(0)c}+2G_{\alpha}^{(0)a}f^{(0){I}}+d^{bca}G_{\alpha}^{(0)b}f^{(0)c}\frac{}{}\right)\delta(p^{2})$
			$\displaystyle+\frac{1}{2p_{n}}if^{bca}\left(\bar{\Pi}^{(1)b\mu}\left[p_{n}f^{(0)c}\delta(p^{2})\right]-\Pi^{(1)b}_{n}\left[\bar{p}^{\mu}f^{(0)c}\delta(p^{2})\right]\frac{}{}\right)$
			$\displaystyle+if^{bca}p^{\mu}\left[\Pi_{\nu}^{(1)b}\left(p^{\nu}f^{(0)c}\right)\right]\delta^{\prime}(p^{2}),$

which satisfy the corresponding transport equations

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle\partial^{x}_{\mu}{\mathscr{J}}^{(1){I}\mu}+\frac{1}{N}G_{\mu}^{(0)a}{\mathscr{J}}^{(1)a\mu},$		(65)
	$\displaystyle 0$	$\displaystyle=$	$\displaystyle{\mathscr{D}}^{ac}_{\mu}{\mathscr{J}}^{(1)c\mu}+2G_{\mu}^{(0)a}{\mathscr{J}}^{(1){I}\mu}+d^{bca}G_{\mu}^{(0)b}{\mathscr{J}}^{(1)c\mu}+if^{bca}G_{\mu}^{(1)b}{\mathscr{J}}^{(0)c\mu}.$		(66)

In order to attain all the results above, we have used the Eq.(4) repeatedly. We note that the singlet distribution $f^{(0){I}}$ and multiplet distribution $f^{(0){a}}$ are totally coupled with each other even in the zeroth order transport equation, which displays the much complexity for non-Abelian chiral kinetic equation, in comparison with chiral kinetic equation in Abelian gauge field.