# Supplementary Material for
“The role of electron polarization on electron driven spin diffusion”

Alessandro Chessari1, Samuel F. Cousin2, Sami Jannin1 and Quentin Stern1
1Univ. Lyon, CNRS, ENS Lyon, UCBL, Université de Lyon, CRMN UMR 5082, 69100
Villeurbanne, France
2Aix Marseille Univ, CNRS, ICR, 13397 Marseille, France

##  I Polarization quantification

Fig. 1(a) shows the proton polarization building up under DNP for positive and
negative DNP, monitored using 0.1° pulses. At the same time, the signal
broadens and sharpens for positive and negative DNP, respectively, as visible
in Fig. 1(b). This is the consequence of intense radiation damping which
causes a change in apparent nuclear $T_{2}$ [1]. Because the signal loss
during the spectrometer dead time is not negligible, the change in apparent
nuclear $T_{2}$ induces an underestimation and an overestimation of the
positive and negative polarization, respectively. The negative polarization is
so strongly overestimated that is reaches an nonphysical absolute value above
$100\%$. The gap between the positive and negative DNP polarization is mostly
an artefact, although a precise quantification is not available in the
literature. However, we know from previous works that the steady-state DNP
polarization is on the other of $70\%$[2], which corresponds roughly to the
average between the negative and positive polarization recorded here.

(a)

(b)

Figure 1: (1(a)) Proton polarization along time under DNP in absolute value
both for positive and negative DNP. (3(b)) Full width at half maximum (FWMH)
of the signal along time for both for positive and negative DNP. The
(constant) FWHM of the thermal equilibrium signal is shown for comparison.

##  II Determination of the saturation width

All HypRes experiments in this study were recorded using a saturation scheme
consisting of 19 sinc excitation pulses of 100 $\mu s$ with a nutation
frequency $\omega_{1}/2\pi=51$ kHz, descritized in 1000 points and separated
by random delays on the order of $\approx 10$ ms (10.100, 100.105, 100.023,
100.157, 100.101, 100.011, 100.073, 100.054, 100.097, 100.047, 100.026,
100.084, 100.065, 100.105, 100.023, 100.157, 100.101 and 100.011 ms). The sinc
pulse shape was generated using Topspin’s built-in functions and chosen to
excite over a bandwidth of 200 kHz, with $n=10$. The profile of the saturation
scheme was measured experimentally and confirmed by simulation.

Experimentally, we first recorded a reference spectrum without saturation
(faint gray lines on Fig. 2) after 10 s of DNP at 1.2 K. Then we repeated the
measurement using the saturation scheme described above varying the carrier
frequency of the saturation pulses from -300 to +300 MHz and obtained the
traces represented in black on Fig. 2(a). The blue line on Fig. 2(a) is the
integral which shows that the saturation width has indeed a bandwidth of 200
kHz.

The excitation profile of the pulse was simulated using a single-spin Hilbert
space with the time dependent Hamiltonian

$\hat{H}(t)=\Delta\Omega\hat{I}_{z}+\omega_{1}A(t)\\{\hat{I}_{x}\cos{\phi(t)}+\hat{I}_{y}\sin{\phi(t)}\\}$
(1)

where $\Delta\Omega$, $\omega_{1}$, $A(t)$ and $\phi(t)$ are the offset
between the Larmor frequency of the spin and the carrier frequency of the
pulse, the nutation frequency of the pulse and the time dependent pulse
amplitude (between 0 and 1) and the phase of the pulse generated by Topspin.
The initial state of the density matrix was assumed to be
$\hat{\rho}_{0}=\hat{I}_{z}$. It was propagated under the time dependent
Hamiltonian during time steps of 0.1 $\mu$s (corresponding to the
discretization of the pulse) using the sandwich formula
$\hat{\rho}(t+dt)=\exp{(+i\hat{H}(t)dt)}\hat{\rho}(t)\exp{(-i\hat{H}(t)dt)}$.
The remaining polarization along the $z$-axis after the pulse was computed
with the trace $P_{z}=\textrm{Tr}(\hat{I}_{z}\hat{\rho})/|I|$. The simulation
was repeated for 300 offset frequencies $\Delta\Omega$ between -300 and +300
MHz. Plotting the remaining polarization $P_{z}(\Delta\Omega)$ against the
offset frequency $\Delta\Omega$ gives the profile of the remaining
polarization along the $z$-axis for a single saturation pulse. To take into
account the fact that the saturation scheme consists of 19 pulses, we assumed
that the magnetization in the transverse plan decays during the delay between
the saturation pulses and so the remaining polarization along the $z$-axis
after the $k$th pulse is
$P_{z}^{final}(\Delta\Omega)=P_{z}(\Delta\Omega)^{k}$. The resulting simulated
saturation profile is shown on Fig. 2(b). It confirms that the excitation
scheme saturates the spins from -100 kHz to +100 kHz.

(a)

(b)

Figure 2: Measured (2(a)) and simulated (2(b)) saturation profile of the train
of sinc pulses. The faint grey line on both plots represents the NMR signal
before saturation. The black signals in (2(a)) are the measured signals as a
function of the offset of the saturation pulses. The bold blue line in (2(a))
is the integral of the measured signals in black. The bold blue line in (2(b))
is the simulated remaining polarization along the $z$-axis after saturation.

##  III Transition rate

The Hamiltonian of the system, in the $\ket{\alpha,\beta}$,
$\ket{\beta,\alpha}$ subspace, consists of two terms:
$\hat{H}(t)=\hat{H}_{D}+\hat{H}_{HFI}(t)$, where

$\displaystyle\hat{H}_{D}=-\frac{1}{2}D_{12}\hat{\sigma}_{x}$
$\displaystyle\hat{H}_{HF}(t)=\frac{1}{2}\Delta(\bar{P}+P^{\prime}(t))\hat{\sigma}_{z}$
(2)

The first term describes the nucleus-nucleus dipolar coupling, with the
coupling given by

$\displaystyle D_{12}$ $\displaystyle=$ $\displaystyle
D_{0}\frac{1}{a_{12}^{3}}\frac{1-3\cos^{2}\theta_{12}}{2}$ $\displaystyle
D_{0}$ $\displaystyle=$ $\displaystyle\frac{\mu_{0}}{4\pi}\hbar\gamma_{n}^{2}$
(3)

where $a_{12}$ the distance between the nuclei and $\cos\theta_{12}$ the angle
between the $z$ direction and the vector linking the two nuclei. The
$\hat{H}_{HF}(t)$ term describes the hyperfine coupling between the two nuclei
and the electron, which is treated semi-classically [3], the hyperfine
coupling is given by

$\displaystyle A_{z,i}$ $\displaystyle=$ $\displaystyle
A_{0}\frac{1}{r_{i}^{3}}\frac{1-3\cos^{2}\theta_{i}}{2}$ $\displaystyle A_{0}$
$\displaystyle=$ $\displaystyle\frac{\mu_{0}}{4\pi}\hbar\gamma_{n}\gamma_{e}$
(4)

where $r_{i}$ and $\theta_{i}$ are the distance and the angle between the
$i-th$ nucleus and the electron. We also introduced an additional parameter
$\Delta=A_{z,1}-A_{z,2}$ which represents the energy splitting between the
nuclear levels in unit of the average electron polarization. We can describe
our system in the interacting picture, given by
$\hat{H}_{0}(t)=\hat{H}_{HF}(t)$ and $\hat{V}=\hat{H}_{D}$. In this
representation we have
$\hat{H}^{(I)}(t)=\hat{U}_{0}^{\dagger}\hat{H}\hat{U}_{0}+i(\partial_{t}\hat{U}_{0}^{\dagger})\hat{U}_{0}$,
with

$\displaystyle\hat{U}_{0}(t)$ $\displaystyle=$
$\displaystyle\text{exp}\left\\{-i\int_{0}^{t}d\tau\hat{H}_{0}(\tau)\right\\}$

The full Hamiltonian in the interacting picture reduces to

$\displaystyle\hat{H}^{(I)}(t)=\hat{V}^{(I)}(t)=\hat{U}_{0}^{\dagger}\hat{V}\hat{U}_{0}=-\frac{1}{2}D_{12}\hat{\sigma}_{x}^{(I)}(t)=-\frac{1}{2}D_{12}[h(\tau)e^{i\Delta\bar{P}\tau}\hat{\sigma}_{+}+h^{*}(\tau)e^{-i\Delta\bar{P}\tau}\hat{\sigma}_{-}]$
(5)

with $h(\tau)=e^{i\Delta\int_{0}^{\tau}P^{\prime}(s)ds}$, where we have used

$\displaystyle\hat{\sigma}_{\pm}^{(I)}(t)$
$\displaystyle=\hat{U}_{0}^{\dagger}(t)\hat{\sigma}_{\pm}\hat{U}_{0}=e^{\pm
i\Delta[\bar{P}t+\int_{0}^{t}P^{\prime}(\tau)d\tau]}\hat{\sigma}_{\pm}$ (6)
$\displaystyle\hat{\sigma}_{x}^{(I)}(\tau)$
$\displaystyle=\hat{\sigma}_{+}^{(I)}(t)+\hat{\sigma}_{-}^{(I)}(t)$ (7)

From the Schrödinger equation, we get for the time evolution operator, for the
first-order in the perturbation parameter $D_{12}$

$\displaystyle\hat{U}^{(I)}(t)$ $\displaystyle=$ $\displaystyle
1-i\int_{0}^{t}d\tau\hat{V}^{(I)}(\tau)+\mathcal{O}(D_{12}^{2})$ (8)

To obtain the transition rate probability between the state
$\ket{\alpha,\beta}$ and $\ket{\beta,\alpha}$ we study the evolution of the
projection on $\ket{\beta,\alpha}$ of the evolution of a state prepared in
$\ket{\alpha,\beta}$, i. e.

$\displaystyle
c_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}(t)=\braket{\beta,\alpha}{\hat{U}(t)}{\alpha,\beta}=\tensor{[}_{I}]{\braket{\beta,\alpha}{\hat{U}^{(I)}(t)}{\alpha,\beta}}{{}_{I}}$
(9)

Since the perturbation is off diagonal, we can sum over all the final state
and, making use of the completeness relation
$\sum_{\ket{m}=\ket{\alpha,\beta},\ket{\beta,\alpha}}\ket{m}\bra{m}=1$, write
for the transition probability

$\displaystyle P_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}$
$\displaystyle=|c_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}(t)|^{2}=\sum\limits_{\ket{m}=\ket{\alpha,\beta},\ket{\beta,\alpha}}|c_{\ket{\alpha,\beta}\rightarrow\ket{m}}|^{2}=\sum\limits_{\ket{m}=\ket{\alpha,\beta},\ket{\beta,\alpha}}\tensor{[}_{I}]{\braket{\alpha,\beta}{\hat{U}^{(I)}(t)}{m}}{{}_{I}}\tensor{[}_{I}]{\braket{m}{\hat{U}^{(I)}(t)}{\alpha,\beta}}{{}_{I}}$
(10)
$\displaystyle=\text{Tr}\\{\rho_{0}^{(I)}\hat{U}^{(I)}(t)\hat{U}^{(I)}(t)\\}=\frac{D_{12}^{2}}{4}\int_{0}^{t}\int_{0}^{t}dt^{\prime}dt^{\prime\prime}\langle\hat{\sigma}_{x}(t^{\prime})\hat{\sigma}_{x}(t^{\prime\prime})\rangle+\mathcal{O}(D_{12}^{3})$
(11)

where the last average is calculated with respect the density matrix
$\hat{\rho}_{0}=\ket{\alpha,\beta}\bra{\alpha,\beta}$. We can now perform the
variable change $s=\frac{1}{2}(t^{\prime}+t^{\prime\prime})$,
$\tau=t^{\prime\prime}-t^{\prime}$

$\displaystyle
P_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}=\frac{D_{12}^{2}}{4}\int_{0}^{t}ds\int_{-B(s)}^{B(s)}d\tau\langle\hat{\sigma}_{x}(s+\tau/2)\hat{\sigma}_{x}(s-\tau/2)\rangle$
(12)

with

$B(s)=\begin{cases*}s-t&if $s\leqslant\frac{t}{2}$,\\\ s&if
$s>\frac{t}{2}$\end{cases*}$ (13)

Now, since the time evolution operator contains a random process, the average
in Eq.(12) has to be intended also as an average over all the realization of
the process $P^{\prime}(t)$. With this in mind, we can exploit the time
translation symmetry of the Hamiltonian and, in the limit of large $t$, the
integration gives

$\displaystyle
W_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}=\frac{d}{dt}P_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}=\frac{D_{12}^{2}}{4}\int_{-\infty}^{\infty}d\tau\langle\hat{\sigma}_{x}(\tau)\hat{\sigma}_{x}(0)\rangle_{P^{\prime}}$
(14)

The result is further simplified by the algebra of the Pauli matrices,

$\displaystyle
W_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}=\frac{D_{12}^{2}}{4}\int_{-\infty}^{\infty}d\tau\langle\hat{\sigma}_{-}(\tau)\hat{\sigma}_{+}(0)\rangle_{P^{\prime}}$
(15)

To calculate this correlator, we need to break down the calculation into two
steps that will be addressed in the two following sections:

  * •

First, we need to average over the density matrix which did not take into
account until now the effects of the environment. In particular, the density
matrix is $\rho_{0}$ only at $t\rightarrow-\infty$, and its actual time
evolution can be obtained tracing out the environmental degrees of freedom. To
do this, we will assume a particular realization of the process
$P^{\prime}(t)$ in the following section. The key part of this section will be
the a Lindblad Master equation and the quantum regression theorem

  * •

Then, we need to do also the averaging over all the processes taking into
account the stochastic behaviour of the electronic polarization, i.e. the
electronic flip-flop

##  IV Nuclear Spin Relaxation

Following [4] we will first describe generally how, for a system coupled to a
bath, the environment degrees of freedom can be traced out. Our goal here is
to describe nuclear relaxation for a pair of nuclear spins in a bath of
surrounding nuclei coupled to the system via the dipolar coupling. We will not
consider any relaxation mechanism leading our subsystem to a paired nuclear
spin state since these states are split away in energy due to the Zeeman
interacton. To do so, the nuclear bath will be described with a Jaynes-
Cummings model, a well-known model in quantum optics [5]. We will show that
within this simple model we can derive the Lindblad operator and the Lindblad
Master equation. Since we are considering a reduced system with a basis
composed of only unpaired nuclear spins, our result will be slightly different
from the form in which the Lindblad Master Equation as found in the NMR
literature [6, 7, 8] since the Lindblad operator will be now applied to
$\hat{\sigma}_{+}=\hat{I}^{+}_{1}\hat{I}^{-}_{2}$ and
$\hat{\sigma}_{-}=\hat{I}^{+}_{2}\hat{I}^{-}_{1}$ rather than on the single
spin raising and lowering operators $\hat{I}^{+}_{k}$ and $\hat{I}^{-}_{k}$.
Let us assume that our total system lives in the Hilbert space
$\mathcal{H}_{T}$ divided into the environment space living in
$\mathcal{H}_{E}$ and our system living in $\mathcal{H}$. The Liouville-Von
Neumann equation gives us the evolution of the density matrix

$\displaystyle\dot{\hat{\rho}}_{T}(t)=-i[\hat{H}_{T},\hat{\rho}_{T}(t)]$ (16)

where the total system Hamiltonian $\hat{H}_{T}$ can be written as
$\hat{H}_{T}=\hat{H}\otimes\mathds{1}_{E}+\mathds{1}\otimes\hat{H}_{E}+\alpha\hat{H}_{SE}$,
being $\hat{H}_{E}\in\mathcal{L}[\mathcal{H}_{E}]$ the Hamiltonian of the
environment and $\hat{H}_{SE}\in\mathcal{L}[\mathcal{H}_{T}]$ the Hamiltonian
of the system-environment interaction, which is characterised by the constant
$\alpha$. We can define the interacting picture by $H+H_{E}$, such that the
time evolution of the total density matrix is given by

$\displaystyle\frac{d}{dt}\hat{\rho}^{(I)}_{T}(t)$ $\displaystyle=$
$\displaystyle-i\alpha[\hat{H}_{SE}^{(I)}(t),\hat{\rho}^{(I)}_{T}(t)]$ (17)

This equation can be formally solved to get

$\displaystyle\hat{\rho}_{T}^{(I)}(t)$ $\displaystyle=$
$\displaystyle\hat{\rho}_{T}^{(I)}(0)-i\alpha\int_{0}^{t}d\tau[\hat{H}_{SE}^{(I)}(\tau),\hat{\rho}^{(I)}_{T}(\tau)]$
(18)

Then, inserting Eq. (17) into Eq. (18)

$\displaystyle\frac{d\hat{\rho}^{(I)}_{T}(t)}{dt}$ $\displaystyle=$
$\displaystyle-i\alpha[\hat{H}^{(I)}_{SE}(t),\hat{\rho}^{(I)}_{T}(0)]-\alpha^{2}\int_{0}^{t}d\tau[\hat{H}_{SE}^{(I)}(t),[\hat{H}_{SE}^{(I)}(\tau),\hat{\rho}_{T}^{(I)}(\tau)]]$
(19)

Inserting in the latter again the zero-th order expansion for
$\rho^{(I)}(\tau)$ around $t$ we get

$\displaystyle\frac{d\hat{\rho}^{(I)}_{T}(t)}{dt}$ $\displaystyle=$
$\displaystyle-i\alpha[\hat{H}^{(I)}_{SE}(t),\hat{\rho}^{(I)}_{T}(0)]-\alpha^{2}\int_{0}^{t}d\tau[\hat{H}^{(I)}_{SE}(t),[\hat{H}_{SE}^{(I)}(\tau),\hat{\rho}^{(I)}_{T}(t)]]+\mathcal{O}(\alpha^{3})$
(20)

Finally, ignoring high-orders terms, the evolution of the density matrix of
the system may be obtained tracing out the environmental degrees of freedom

$\displaystyle\frac{d\hat{\rho}^{(I)}(t)}{dt}$ $\displaystyle=$
$\displaystyle-i\alpha\text{Tr}_{E}[\hat{H}_{SE}^{(I)}(t),\hat{\rho}^{(I)}_{T}(0)]-\alpha^{2}\int_{0}^{t}d\tau\text{Tr}_{E}[\hat{H}_{SE}^{(I)}(t),[\hat{H}_{SE}^{(I)}(\tau),\hat{\rho}^{(I)}_{T}(t)]]$
(21)

Now, the interacting system-environment Hamiltonian can be written, without
loss of generalities as

$\displaystyle\hat{H}_{SE}=\sum\limits_{i}S_{i}\otimes E_{i}$ (22)

where $S_{i}$ and $E_{i}$ are operators acting on the system and on the
environment, respectively. Given this decomposition, Eq. (17) can be always be
rewritten as [4]

$\displaystyle\frac{d\hat{\rho}^{(I)}(t)}{dt}$ $\displaystyle=$
$\displaystyle-\alpha^{2}\int_{0}^{t}d\tau\text{Tr}_{E}[\hat{H}_{SE}^{(I)}(t),[\hat{H}_{SE}^{(I)}(\tau),\hat{\rho}^{(I)}_{T}(t)]]$
(23)

We may also suppose that the system and the environment are noncorrelated
during all the time evolution, i.e. that the correlation and relaxation time
scale are much smaller than that of the system. In this regime we can assume
the environment state to be stationary and decoupled from the system state,
then $\hat{\rho}_{E}^{(I)}=\hat{\rho}_{E}(0)$ and
$\hat{\rho}_{T}^{(I)}=\hat{\rho}^{(I)}(t)\otimes\hat{\rho}_{E}(0)$. Extending
the upper integration limit to infinity and by the change of variable
$\tau\rightarrow t-\tau$ leads to

$\displaystyle\frac{d\hat{\rho}^{(I)}(t)}{dt}$ $\displaystyle=$
$\displaystyle-\alpha^{2}\int_{0}^{\infty}d\tau\text{Tr}_{E}[\hat{H}_{SE}^{(I)}(t),[\hat{H}_{SE}^{(I)}(t-\tau),\hat{\rho}^{(I)}(t)\otimes\hat{\rho}_{E}(0)]]$
(24)

At this point we need to choose the environment Hamiltonian and the coupling
term between the system and the bath for our specific case. The simplest idea
is to consider all the surrounding nuclei as a set of quantum harmonic
oscillators, following the so-called Jaynes-Cummings model. The bath
Hamiltonian is then
$\hat{H}_{E}=\sum\limits_{j}\nu_{j}\hat{a}^{\dagger}_{j}\hat{a}_{j}$ [9],
where $\nu_{j}$ are the frequency and $\hat{a}^{\dagger}_{j}$ the annihilation
operator for the $j$th-oscillator. We assume that the environment induces
transitions in the system by coupling the raising and lowering operators
$\hat{\sigma}_{+}$ and $\hat{\sigma}_{-}$. On the environment the
corresponding effect is the annihilation and creation of an energy mode of the
oscillators. The system-environment interaction may then be written as follows

$\displaystyle\hat{H}_{SE}=\sum\limits_{j}g_{j}(\hat{a}^{\dagger}_{j}\hat{\sigma}_{-}+\hat{a}_{j}\hat{\sigma}_{+})=\hat{F}\hat{Q}^{\dagger}+h.c.$
(25)

where $g_{j}$ represents the coupling strength between the oscillators and the
system. The operator $F$ and $Q$ are defined as follows

$\hat{F}=\hat{\sigma}_{-}\quad\hat{Q}=\sum\limits_{k}g_{k}\hat{a}$ (26)

That, in the interacting picture of $H+H_{E}$, writes

$\hat{F}^{(I)}(t)=e^{-i\omega_{*}t}\hat{F}=e^{-i\Delta(\bar{P}t+\int_{0}^{t}P^{\prime}(\tau)d\tau)}\hat{\sigma}_{-}\qquad\qquad\hat{Q}^{(I)}(t)=\sum\limits_{k}g_{k}\hat{a}_{k}(t)=\sum\limits_{k}g_{k}\hat{a}_{k}e^{-i\nu_{k}t}$
(27)

where, to make the notation easier, we have defined
$\omega_{*}=\Delta\bar{P}+\frac{1}{t}\int_{0}^{t}P^{\prime}(\tau)d\tau$.
Taking into account the following relations, where $n_{B}$ is the Bose-
Einstein distribution

$\displaystyle\text{Tr}_{E}\left\\{\hat{\rho}_{E}\hat{a}_{k}(t)\hat{a}_{k^{\prime}}(t-\tau)\right\\}=0$
$\displaystyle\text{Tr}_{E}\left\\{\hat{\rho}_{E}\hat{a}_{k}^{\dagger}(t)\hat{a}_{k^{\prime}}(t-\tau)\right\\}=\delta_{kk^{\prime}}\langle\hat{a}^{\dagger}_{k}\hat{a}_{k}\rangle
e^{+i\nu_{k}t}e^{-i\nu_{k^{\prime}}(t-\tau)}=\delta_{kk^{\prime}}\langle\hat{a}^{\dagger}_{k}\hat{a}_{k}\rangle
e^{+i\nu_{k}\tau}=\delta_{kk^{\prime}}n_{B}(\nu_{k})e^{i\nu_{k}\tau}$ (28)
$\displaystyle\text{Tr}_{E}\left\\{\hat{\rho}_{E}\hat{a}_{k}(t)a_{k^{\prime}}^{\dagger}(t-\tau)\right\\}=\delta_{kk^{\prime}}\langle\hat{a}_{k}\hat{a}^{\dagger}_{k}\rangle
e^{-i\nu_{k}t}e^{+i\nu_{k^{\prime}}(t-\tau)}=\delta_{kk^{\prime}}\langle\hat{a}_{k}\hat{a}^{\dagger}_{k}\rangle
e^{-i\nu_{k}\tau}=\delta_{kk^{\prime}}(1+n_{B}(\nu_{k}))e^{-i\nu_{k}\tau}$
(29)
$\displaystyle\text{Tr}_{E}\left\\{\hat{\rho}_{E}\hat{a}^{\dagger}_{k}(t)a_{k^{\prime}}^{\dagger}(t-\tau)\right\\}=0$
(30)

Eq. (24) becomes (here we avoid the superscript indicating the interacting
picture for simplicity)

$\frac{d\hat{\rho}(t)}{dt}=\int_{0}^{\infty}d\tau\\{\hat{F}^{{\dagger}}(t-\tau)\rho(t)\hat{F}(t)\text{Tr}_{E}[\hat{Q}(t-\tau)\rho_{E}\hat{Q}^{{\dagger}}(t)]-\hat{F}(t)\hat{F}^{\dagger}(t-\tau)\rho(t)\text{Tr}_{E}[\hat{Q}^{\dagger}(t)\hat{Q}(t-\tau)\rho_{E}]+\hat{F}(t-\tau)\rho(t)\hat{F}^{\dagger}(t)\text{Tr}_{E}[\hat{Q}^{\dagger}(t-\tau)\rho_{E}\hat{Q}(t)]-\hat{F}^{\dagger}(t)\hat{F}(t-\tau)\rho(t)\text{Tr}_{E}[\hat{Q}(t)\hat{Q}^{\dagger}(t-\tau)\rho_{E}]+h.c\\}$
(31)

that further simplify to

$\frac{d\hat{\rho}(t)}{dt}=\int_{0}^{\infty}d\tau\\{[\hat{F}^{{\dagger}}(t-\tau)\rho(t)\hat{F}(t)-\hat{F}(t)\hat{F}^{\dagger}(t-\tau)\rho(t)]\langle\hat{Q}^{\dagger}(t)\hat{Q}(t-\tau)\rangle_{E}+[\hat{F}(t-\tau)\rho(t)\hat{F}^{\dagger}(t)-\hat{F}^{\dagger}(t)\hat{F}(t-\tau)\rho(t)]\langle\hat{Q}(t)\hat{Q}^{\dagger}(t-\tau)\rangle_{E}+h.c\\}$
(32)

and finally to

$\displaystyle\frac{d\hat{\rho}(t)}{dt}=$
$\displaystyle[2\hat{F}^{{\dagger}}(t)\rho(t)\hat{F}(t)-{\hat{F}(t)\hat{F}^{\dagger}(t),\rho(t)}]\int_{0}^{\infty}d\tau
e^{-i\omega_{*}\tau}\langle\hat{Q}^{\dagger}(t)\hat{Q}(t-\tau)\rangle_{E}$
(33)
$\displaystyle+[2\hat{F}(t)\rho(t)\hat{F}^{\dagger}(t)-{\hat{F}^{\dagger}(t)\hat{F}(t),\rho(t)}]\int_{0}^{\infty}d\tau
e^{i\omega_{*}\tau}\langle\hat{Q}(t)\hat{Q}^{\dagger}(t-\tau)\rangle_{E}$ (34)

Now, introducing the rate constants

$\displaystyle\Gamma_{+}(\omega_{*})$ $\displaystyle=4\int_{0}^{\infty}d\tau
e^{-i\omega_{*}\tau}\langle\hat{Q}^{\dagger}(t)\hat{Q}(t-\tau)\rangle_{E}=4\pi
g(\omega_{*})^{2}_{*}\rho_{D}(\omega_{*})n_{B}(\omega_{*})$ (35)
$\displaystyle\Gamma_{-}(\omega_{*})$ $\displaystyle=4\int_{0}^{\infty}d\tau
e^{i\omega_{*}\tau}\langle\hat{Q}(t)\hat{Q}^{\dagger}(t-\tau)\rangle_{E}=4\pi
g(\omega_{*})^{2}_{*}\rho_{D}(\omega_{*})(1+n_{B}(\omega_{*}))$ (36)

where we have considered the levels of the bath as a continuum of states with
a density of states give by $\rho_{D}$. The result, back in the Schrödinger
picture is

$\displaystyle\frac{d\hat{\rho}}{dt}$ $\displaystyle=$
$\displaystyle-i[\hat{H}(t),\hat{\rho}(t)]+\frac{\Gamma^{*}_{+}}{2}\mathcal{D}(\hat{\sigma}_{+})\hat{\rho}(t)+\frac{\Gamma^{*}_{-}}{2}\mathcal{D}(\hat{\sigma}_{-})\hat{\rho}(t)$
(37)

where we have defined $\Gamma^{*}_{+}=\Gamma_{+}(\omega_{*})$,
$\Gamma^{*}_{-}=\Gamma_{-}(\omega_{*})$ and the dissipator

$\displaystyle\mathcal{D}(\hat{X})\hat{\rho}$ $\displaystyle=$
$\displaystyle\hat{X}\hat{\rho}\hat{X}^{\dagger}-\frac{1}{2}\\{\hat{X}^{\dagger}\hat{X},\hat{\rho}\\}$
(38)

With a similar procedure, considering a coupling with $\hat{\sigma}_{z}$ in
Eq. (25) (we will not enter the details of the dephasing term, we will suppose
that mechanisms that can lead to such a coupling exist) we can finally get the
Master equation in the Lindblad form $\partial_{t}\rho(t)=\mathcal{L}\rho(t)$
where

$\displaystyle\mathcal{L}\rho(t)$ $\displaystyle=$
$\displaystyle-i[\hat{H}(t),\hat{\rho}(t)]+\frac{\Gamma^{*}_{+}}{2}\mathcal{D}(\hat{\sigma}_{+})\hat{\rho}(t)+\frac{\Gamma^{*}_{-}}{2}\mathcal{D}(\hat{\sigma}_{-})\hat{\rho}(t)+\Gamma_{2}\mathcal{D}(\hat{\sigma}_{z})\hat{\rho}(t)$
(39)

This Master Equation is the same as the Lindblad Master equation as found in
recent NMR literature [6, 7, 8], reduced in the
$\ket{\alpha,\beta},\ket{\beta,\alpha}$ subspace. Our simple model used to
derive it has all the essential ingredients to describe relaxation: the
nuclear relaxation constants $\Gamma^{*}_{+}$, $\Gamma^{*}_{-}$ and
$\Gamma_{2}$. Any other relaxation mechanisms can be taken into account in the
same way we did so far and will lead to corrections to the relaxation rates.
Now, we can exploit the quantum regression theorem [9][5] that states that, if
the equations of motions for the expectation values of a set of system
operators are expressed as

$\displaystyle\frac{d}{dt}\langle\hat{A}_{i}(t)\rangle=\sum\limits_{j}G_{ij}\langle\hat{A}_{j}(t)\rangle$
(40)

then for the correlations hold

$\displaystyle\frac{d}{d\tau}\langle\hat{A}_{i}(t+\tau)\hat{A}_{k}(t)\rangle=\sum\limits_{j}G_{ij}\langle\hat{A}_{j}(t+\tau)\hat{A}_{k}(t)\rangle$
(41)

In the weak coupling regime $D_{12}\ll\Gamma_{2}$ [10], the evolution of
$\langle\hat{\sigma}_{-}(t)\rangle$ is given, through the Lindblad Master
equation presented in Eq. (39)
$\partial_{t}\rho^{(I)}(t)=\mathcal{L}^{(I)}\rho^{(I)}(t)$ (we put ourself in
the interaction picture given by $H_{0}$ as in the last section for
simplicity)

$\displaystyle\frac{d}{dt}\langle\hat{\sigma}_{-}(t)\rangle$
$\displaystyle=\text{Tr}\left\\{\frac{d\hat{\rho}^{(I)}(t)}{dt}\hat{\sigma}_{-}^{(I)}(t)\right\\}$
(42)
$\displaystyle=\text{Tr}\left\\{\left[-i[\hat{V}^{(I)}(t),\hat{\rho}^{(I)}(t)]+\frac{\Gamma^{*}_{+}}{2}\mathcal{D}(\hat{\sigma}_{+}^{(I)})\hat{\rho}^{(I)}(t)+\frac{\Gamma^{*}_{-}}{2}\mathcal{D}(\hat{\sigma}_{-}^{(I)})\hat{\rho}^{(I)}(t)+\Gamma_{2}\mathcal{D}(\hat{\sigma}_{z}^{(I)})\hat{\rho}^{(I)}(t)\right]\hat{\sigma}_{-}^{(I)}(t)\right\\}$
(43)
$\displaystyle=-\left(\frac{\Gamma^{*}_{+}}{2}+\frac{\Gamma^{*}_{-}}{2}+\Gamma_{2}\right)\langle\hat{\sigma}_{-}(t)\rangle$
(44)

By means of (7) and of (41), neglecting $\Gamma_{+}^{*}$ and $\Gamma_{-}^{*}$
in front of $\Gamma_{2}$ (i.e., assuming that nuclear spin-spin relaxation is
fast in front of nuclear spin-lattice relaxation)

$\displaystyle\langle\hat{\sigma}_{-}(\tau)\hat{\sigma}_{+}(0)\rangle$
$\displaystyle=$ $\displaystyle
e^{-\Gamma_{2}|\tau|}e^{-i\Delta\bar{P}\tau-i\Delta\int_{0}^{\tau}P^{\prime}(s)ds}$
(46)

where we have set the initial value
$\langle\hat{\sigma}_{-}(0)\hat{\sigma}_{+}(0)\rangle=1$, that corresponds to
our choice for the density matrix
$\hat{\rho}_{0}=\ket{\alpha,\beta}\bra{\alpha,\beta}$. Similarly, we obtain
that all the others correlations are null.

##  V Electronic flip-flop

Finally, we will consider the electron motion not as a deterministic signal,
but rather as undergoing stochastic flip-flop, i.e. instantaneously changing
in its polarisation, along the $z$ axis with average rate $\Gamma_{c}$. If
$T_{1e}$ and $T_{2e}$ are the spin-lattice and spin-spin relaxation time of
the electron, respectively, the $\Gamma_{c}$ in an insulating solid is [11]

$\Gamma_{c}=\frac{1}{T_{1e}}+\frac{1}{T_{2e}}$ (47)

In our conditions, the $1/T_{2e}$ is large in front of $1/T_{1e}$ due to
strong electron-electron interactions and so we have $\Gamma_{c}=1/T_{2e}$.
However, the derivation presented here holds regardless of the origin of the
fluctuations of the electron spin state. Being $\bar{P}$ the average electron
polarization, we can write

$P(t)=P^{\prime}(t)+\bar{P}$ (48)

being $P^{\prime}(t)$ the unbiased process. From [3] we have that the power
spectral density of the unbiased process is

$S_{P^{\prime}}(\omega)=(1-\bar{P}^{2})\frac{2\Gamma_{c}}{\omega^{2}+\Gamma_{c}^{2}}$
(49)

The autocorrelation function for the process $\Delta P^{\prime}(t)$ is then
[3]

$\displaystyle R(\tau)=\Delta^{2}\langle P^{\prime}(0)P^{\prime}(\tau)\rangle$
$\displaystyle=$ $\displaystyle\Delta^{2}(1-\bar{P}^{2})e^{-\Gamma_{c}|\tau|}$
(50)

Since the electron motion is stochastic, we are interested in the expectation
value in Eq. (15) has to be taken as an ensemble average over all the possible
realization of the signal $P^{\prime}(t)$ from (46). Hence, we are interested
in the quantity
$\bigl{<}f(\tau)\bigr{>}_{P^{\prime}}=\bigl{<}e^{-i\Delta\int_{0}^{\tau}P^{\prime}(t)dt}\bigr{>}_{P^{\prime}}$.
We can start from the differential equation $\dot{f}(\tau)=-i\Delta
P^{\prime}(\tau)f(\tau)$, which, after a single iteration, assumes the form
[12]

$\displaystyle\partial_{\tau}f(\tau)$ $\displaystyle=$ $\displaystyle-i\Delta
P^{\prime}(\tau)-\Delta^{2}\int_{0}^{\tau}P^{\prime}(\tau)P^{\prime}(t)f(t)dt$

exploiting the time reversal symmetry of $P^{\prime}(t)$ we get

$\displaystyle\frac{\partial\langle f(\tau)\rangle}{\partial\tau}$
$\displaystyle=$
$\displaystyle-\int_{0}^{\infty}\theta(\tau-t)R(\tau-t)\langle
f(\tau)\rangle_{P^{\prime}}dt$ (51)

in the Laplace domain, being $\tilde{F}(s)$ the Fourier transform of $\langle
f(\tau)\rangle_{P^{\prime}}$, the latter equation reads
$i\tilde{F}(s)=-\tilde{R}(s)\tilde{F}(s)$, such that

$\displaystyle\tilde{F}(s)=\frac{1}{\tilde{R}(s)+s}$ (52)

Being $\tilde{R}(s)=\Delta^{2}\frac{1-\bar{P}^{2}}{s+\Gamma_{c}}$, we can
easily obtain the Fourier transform of $\langle f(\tau)\rangle_{P^{\prime}}$
if we use the identity for the Fourier transform of time-symmetric functions
$F(\omega)=2\text{Re}\left\\{\tilde{F}(s=i\omega)\right\\}$

$\displaystyle F(\omega)$ $\displaystyle=$
$\displaystyle\frac{2\Gamma_{c}(1-\bar{P}^{2})\Delta^{2}}{(\omega^{2}-(1-\bar{P}^{2})\Delta^{2})^{2}+\omega^{2}\Gamma_{c}^{2}}$
(53)

Exploiting the convolution theorem, Eq. (15) becomes

$\displaystyle W_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}$
$\displaystyle=\frac{D_{12}^{2}}{4}\int_{-\infty}^{\infty}e^{-i\Delta\bar{P}\tau-\Gamma_{2}|\tau|}\langle
f(\tau)\rangle_{P^{\prime}}d\tau$ (54)
$\displaystyle=\frac{D_{12}^{2}}{4}\frac{1}{2\pi}\int_{-\infty}^{\infty}S_{\Gamma_{2}}(\omega^{\prime}-\Delta\bar{P})F(\omega^{\prime})d\omega^{\prime}$
(55)

Where we have introduced the Fourier transform of $e^{-\Gamma_{2}|t|}$, namely

$\displaystyle S_{\Gamma_{2}}(\omega)$ $\displaystyle=$
$\displaystyle\frac{2\Gamma_{2}}{\omega^{2}+\Gamma_{2}^{2}}$ (56)

that is, a Lorentzian shape with FWHM$=2\Gamma_{2}$.
Now, it is clear [13] that
$\frac{1}{2\pi}\int_{-\infty}^{\infty}S_{\Gamma_{2}}(\omega^{\prime}-\Delta\bar{P})F(\omega^{\prime})d\omega^{\prime}=2\text{Re}\left\\{\tilde{F}(s=i\Delta\bar{P}+\Gamma_{2})\right\\}$,
such that [12]

$\displaystyle
W_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}=\frac{D_{12}^{2}}{2}\frac{\Gamma_{c}(1-\bar{P}^{2})\Delta^{2}+\Gamma_{2}(\bar{\Gamma}^{2}+\Delta^{2})}{((1-2\bar{P}^{2})\Delta^{2}+\Gamma_{2}\bar{\Gamma})^{2}+\Delta^{2}\bar{P}^{2}(\Gamma_{2}+\bar{\Gamma})^{2}}$
(57)

where we have introduced $\bar{\Gamma}=\Gamma_{2}+\Gamma_{c}$. The diffusion
coefficient is then given by
$D=W_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}a_{12}^{2}$. The final
result is shown in Fig. 3(a).

For $\bar{P}=0.6$, Fig. 3(a) shows that spin diffusion is particularly
efficient for $r<0.5$ nm. To explain this local maximum inside the spin
diffusion barrier, let us examine a specific case. In Fig. 3(b), we show the
diffusion coefficient for $\bar{P}=0.6$ highlighting the position of the local
maximum, at $r=$0.34\text{\,}\mathrm{n}\mathrm{m}$$. This maximum is the
result of the convolution integral of Eq. (55), which is shown in Fig. 4(b).
Indeed the maximum is reached whenever $F(\omega)$, that represents the
spectrum of the system under the evolution of the dynamic part of
$\hat{H}_{HF}(t)$, is tuned such as to have a maximum in $\Delta\bar{P}$, that
is the spectrum under the static part of the hyperfine coupling.

(a)

(b)

Figure 3: In (3(a)) the diffusion coefficient in function of the distance from
the electron at different electron polarization. In (3(b)) the diffusion
coefficient and its limit cases at $\bar{P}=0.6$, where we have marked the
diffusion coefficient at $r=$0.34\text{\,}\mathrm{n}\mathrm{m}$$. Both plots
are obtained with the following parameters:
$\Gamma_{2}=$9\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$,
$\Gamma_{c}=$2\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$,
$a_{12}=$0.66\text{\,}\mathrm{n}\mathrm{m}$$.

(a)

(b)

Figure 4: In (4(a)) the plot of $F(\omega)$ in function of $\omega$ for
different values of the electron distance. In (4(b)) the plot of $F(\omega)$
and $S_{\Gamma_{2}}(\omega-\Delta\bar{P})$ in function of $\omega$ at a
distance from the electron of $r=$0.34\text{\,}\mathrm{n}\mathrm{m}$$. Both
plots are obtained with the following parameters: $\bar{P}=0.6$,
$\Gamma_{2}=$9\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$,
$\Gamma_{c}=$2\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$,
$a_{12}=$0.66\text{\,}\mathrm{n}\mathrm{m}$$.

As limit cases, we can consider

  * •

The static field limit, $\Gamma_{c}\rightarrow\infty$ and $\bar{P}=1$, for
which

$\displaystyle W_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}$
$\displaystyle=$
$\displaystyle\frac{D_{12}^{2}}{2}\frac{\Gamma_{2}}{\Gamma_{2}^{2}+\Delta^{2}},\qquad\Gamma_{c}\rightarrow\infty$
(58)

in agreement with the literature [14][15]

  * •

The high polarization limit, $\bar{P}\gg 0$. In the limit
$\Gamma_{2}<\Gamma_{c}$ the function $S_{\Gamma_{2}}(\omega-\Delta\bar{P})$
selects only frequencies $\omega\approx\Delta\bar{P}$, then for $\Delta>0$ (i.
e., not too far from the electron, such that
$(1-P^{2})\Delta^{2}\ll\omega^{2}$) we can write

$\displaystyle F(\omega)$ $\displaystyle=$
$\displaystyle\frac{(1-\bar{P}^{2})\Delta^{2}}{\omega^{2}}\frac{2\Gamma_{c}}{\omega^{2}+\Gamma_{c}^{2}},\qquad\bar{P}\gg
0,\,\Delta>0$ (59)

in agreement with the expression found by E. P. Horvitz [3]. In this regime we
can take $\omega^{-2}\approx(\Delta\bar{P})^{-2}$ in the integral (55) and the
final result is given by

$\displaystyle
W_{\ket{\alpha,\beta}\rightarrow\ket{\beta,\alpha}}=\frac{D_{12}^{2}}{4}\frac{1}{2\pi}\int_{-\infty}^{\infty}S_{\Gamma_{2}}(\omega^{\prime}-\Delta\bar{P})F(\omega^{\prime})d\omega^{\prime}=\frac{D_{12}^{2}}{2}\frac{1-\bar{P}^{2}}{\bar{P}^{2}}\frac{\bar{\Gamma}}{\Delta^{2}\bar{P}^{2}+\bar{\Gamma}^{2}},\quad\bar{P}\gg
0,\,\Delta>0$ (60)

The results for both cases are shown in Fig. 3(b), where $D_{stat}$
corresponds to the static limit case $\Gamma_{c}\rightarrow 0$, $\bar{P}=1$
and $D_{dyn}$ is the result for Horvitz’s approximation.

## References

  * Elliott  _et al._ [2021] S. J. Elliott, Q. Stern, M. Ceillier, T. El Daraï, S. F. Cousin, O. Cala, and S. Jannin, Progress in Nuclear Magnetic Resonance Spectroscopy 126, 59 (2021).
  * Stern  _et al._ [2021] Q. Stern, S. F. Cousin, F. Mentink-Vigier, A. C. Pinon, S. J. Elliott, O. Cala, and S. Jannin, Science Advances 7, eabf5735 (2021).
  * Horvitz [1971] E. P. Horvitz, Physical Review B 3, 2868 (1971).
  * Manzano [2020] D. Manzano, AIP Advances 10, 025106 (2020).
  * D.F.Walls [1994] D.F.Walls,  _Quantom Optics_ (Springer, 1994).
  * Bengs and Levitt [2020] C. Bengs and M. H. Levitt, Journal of Magnetic Resonance 310, 106645 (2020).
  * Karabanov  _et al._ [2014] A. Karabanov, G. Kwiatkowski, and W. Köckenberger, Molecular Physics 112, 1838 (2014).
  * Karabanov  _et al._ [2018] A. Karabanov, D. Wiśniewski, F. Raimondi, I. Lesanovsky, and W. Köckenberger, Physical Review A 97, 031404 (2018).
  * Silveri  _et al._ [2017] M. P. Silveri, J. A. Tuorila, E. V. Thuneberg, and G. S. Paraoanu, Reports on Progress in Physics 80 (2017).
  * Milburn [1994] D. F. W. G. J. Milburn,  _Quantum optics_ (Springer, 1994).
  * Pell  _et al._ [2019] A. J. Pell, G. Pintacuda, and C. P. Grey, Progress in nuclear magnetic resonance spectroscopy 111, 1 (2019).
  * Li  _et al._ [2013] J. Li, M. P. Silveri, K. S. Kumar, J. M. Pirkkalainen, A. Vepsäläinen, W. C. Chien, J. Tuorila, M. A. Sillanpää, P. J. Hakonen, E. V. Thuneberg, and G. S. Paraoanu, Nature Communications 4, 1 (2013).
  * Wódkiewicz  _et al._ [1984] K. Wódkiewicz, B. W. Shore, and J. H. Eberly, Phys. Rev. A 30, 2390 (1984).
  * Suter and Ernst [1985] D. Suter and R. Ernst, Physical Review B 32, 5608 (1985).
  * Karabanov  _et al._ [2015] A. Karabanov, D. Wiśniewski, I. Lesanovsky, and W. Köckenberger, Phys. Rev. Lett. 115, 020404 (2015).


# Simultaneous Determination of Fragmentation Functions and Test on Momentum
Sum Rule

Jun Gao1,2, ChongYang Liu1,2, XiaoMin Shen1,2,3, Hongxi Xing4,5,6, Yuxiang
Zhao6,7,8,9 1 INPAC, Shanghai Key Laboratory for Particle Physics and
Cosmology
& School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China
2Key Laboratory for Particle Astrophysics and Cosmology (MOE), Shanghai
200240, China
3Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany
4Key Laboratory of Atomic and Subatomic Structure and Quantum Control (MOE),
Guangdong Basic Research Center of Excellence for Structure and Fundamental
Interactions of Matter, Institute of Quantum Matter, South China Normal
University, Guangzhou 510006, China
5Guangdong-Hong Kong Joint Laboratory of Quantum Matter, Guangdong Provincial
Key Laboratory of Nuclear Science, Southern Nuclear Science Computing Center,
South China Normal University, Guangzhou 510006, China
6Southern Center for Nuclear-Science Theory (SCNT), Institute of Modern
Physics, Chinese Academy of Sciences, Huizhou 516000, China
7Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou, Gansu
730000, China
8University of Chinese Academy of Sciences, Beijing 100049, China
9Key Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle
Physics, Central China Normal University, Wuhan 430079, China

###### Abstract

We perform a simultaneous global analysis of hadron fragmentation functions
(FFs) to various charged hadrons at next-to-leading order in QCD. The world
data set includes results from electron-positron single-inclusive
annihilation, semi-inclusive deep inelastic scattering, as well as proton-
proton collisions including jet fragmentation measurements which lead to
strong constraints on the gluon fragmentations. By carefully selecting hadron
kinematics to ensure the validity of QCD factorization and the convergence of
perturbative calculations, we achieve a satisfying best fit with
$\chi^{2}/$d.o.f.$=0.90$, in the simultaneous extraction of FFs for light
charged hadrons ($\pi^{\pm}$, $K^{\pm}$ and $p/\bar{p}$). The total momentum
of $u$, $d$ quarks and gluon carried by light charged hadrons have been
determined precisely. That urges future precision measurements on
fragmentation to neutral hadrons, which are crucial for the test of
fundamental sum rules in QCD fragmentation.

Introduction.– Single inclusive hadron production has been extensively
measured in electron-positron single-inclusive annihilation (SIA), semi-
inclusive deep inelastic scattering (SIDIS), as well as proton-proton
collisions ($pp$) with the collision energy ranging from $\sqrt{s}=2.23$ GeV
at BEPCII Ablikim et al. (2023) to $\sqrt{s}=13$ TeV at the Large Hadron
Collider Aaij et al. (2023). Such processes can be factorized in terms of the
hadron fragmentation functions (FFs) within the framework of factorization
theorem of quantum chromodynamics (QCD) Collins et al. (1989). The FFs are
crucial for comprehending color confinement in QCD, as they not only represent
the final state hadronization of quarks and gluons into color singlet hadrons,
but also play a vital role in accurately probing the internal nucleon
structure at the level of specific quark flavors in the initial state Metz and
Vossen (2016), which becomes especially critical in the era of electron-ion
colliders Accardi et al. (2016); Anderle et al. (2021), because of the
unprecedented demanding on the precise determination of various parton
distribution functions (PDFs). Moreover, the determination of FFs provides
essential baseline information in the identification of fundamental transport
properties of quark-gluon plasma created in heavy ion collisions Burke et al.
(2014).

Analogous to PDFs, the fundamental properties of FFs are rooted in their
physical interpretation as number densities and the associated sum rules,
which stem from the conservation laws in QCD Collins and Soper (1982). Recent
investigations have renewed the interest regarding the validity of these
properties for both single-hadron Collins and Rogers (2023) and dihadron FFs
Pitonyak et al. (2023) based on their operator definitions. The information
encoded in these properties of FFs is essential for understanding the
mechanism of hadronization. However, due to their nonperturbative nature and
the complexity involved in defining all out-states Collins and Rogers (2023),
FFs cannot be computed and their fundamental properties cannot be justified
from first principles. Therefore, it is of paramount importance to rigorously
determine the FFs through global analyses of worldwide data using the
collinear factorization formalism. Significant efforts have been dedicated to
the phenomenological extraction of FFs, utilizing various data samples, as
demonstrated by the works of DSS de Florian et al. (2007), HKN Hirai et al.
(2007), AKK Albino et al. (2008), NNFF Bertone et al. (2018), MAPFF Khalek et
al. (2021), and JAM Moffat et al. (2021). Nevertheless, one of the least
explored directions is the simultaneous fit by including different identified
final state hadrons.

In this Letter, we present novel results on FFs for light charged hadrons,
specifically $\pi^{\pm}$, $K^{\pm}$ and $p/\bar{p}$. Our analysis represents a
simultaneous determination at next-to-leading order (NLO) in QCD using global
datasets of SIA, SIDIS, and $pp$ collisions. Several innovations are
highlighted in our analysis. A stringent selection criterion has been
implemented on the kinematics of the fragmentation processes to ensure the
validity of leading twist collinear factorization and the associated
perturbative calculations of QCD. Additionally, we have incorporated residual
theory uncertainties into the analysis. In contrast to using Mellin
transforms, we have utilized a fast interpolation technique for the
calculations of the cross sections, significantly increasing the efficiency of
the global fit. The $pp$ datasets consist of measurements on production cross
section ratios of various charged hadrons, as well as measurements on charged
hadron production in jets at the LHC. The latter datasets provide strong
constraints on the gluon FFs. To the best of our knowledge, our work
represents the first joint determination of FFs of various charged hadrons
from a global analysis, including data from SIA, SIDIS and $pp$ collisions. It
also marks the first inclusion of jet fragmentation measurements in a global
analysis for light charged hadrons. The comprehensive analysis provides a
state-of-the-art determination of FFs, allowing for a test on the fundamental
law of momentum sum rule.

Theoretical setup and data characteristics.– The global analysis of FFs
requires a parametrization form at the initial scale $Q_{0}$. We take the
following form for charged hadrons

$zD^{h}_{i}\left(z,Q_{0}\right)=z^{\alpha^{h}_{i}}(1-z)^{\beta^{h}_{i}}\exp\left(\sum_{n=0}^{m}a^{h}_{i,n}(\sqrt{z})^{n}\right),$
(1)

where $\\{\alpha,\beta,a_{n}\\}$ are free parameters in the fit, and $i$ and
$h$ represent the flavor of the parton and the hadron, respectively. The
$Q_{0}$ is set to be 5 GeV and a zero-mass scheme is utilized for heavy quarks
with active quark-flavors $n_{f}=5$. One advantage of the above
parametrization form is that the FFs are positively defined, eliminating the
need for additional positivity constraints. The light quark fragmentation is
categorized into favored and unfavored components according to the hadron
constituents, while isospin or flavor asymmetry is allowed. For instance, we
assume $D^{\pi^{+}}_{u}$ and $D^{\pi^{+}}_{\bar{d}}$ have the same shape but
independent normalization, similar for $D^{K^{+}}_{u}$ and
$D^{K^{+}}_{\bar{s}}$. Furthermore, the FFs of negative-charged hadrons are
related to the positive ones via charge conjugation. The degree of polynomials
$m$ has been increased till no significant improvements of fit are observed,
with the final values varying from 0 to 2 depending on the flavors of parton
and hadron. The total number of free parameters is 63 for $\pi^{+}$, $K^{+}$,
and $p$ together.

The FFs are evolved to higher scales using two-loop time-like splitting
kernels to maintain consistency with the NLO analysis. The splitting functions
were calculated in Refs. Stratmann and Vogelsang (1997) and have been
implemented in HOPPET Salam and Rojo (2009), which is employed in the
analysis. Theoretical calculations of the differential cross sections are
carried out at NLO in QCD with the FMNLO program Liu et al. (2023), and are
accelerated with the interpolation grid and fast convolution algorithms. This
allows for an efficient scan over the parameter space of large dimension for
hundreds of thousands of times. For calculations involving initial hadrons,
CT14 NLO parton distribution functions Dulat et al. (2016) are used with
$\alpha_{S}(M_{Z})=0.118$ Zyla et al. (2020). The central values of the
renormalization and fragmentation scales ($\mu_{R,0}$ and $\mu_{D,0}$) are set
to the momentum transfer $Q$ for both SIA and SIDIS. In the case of $pp$
collisions, the central values of the factorization scale ($\mu_{F,0}$) and
renormalization scales are set to the sum of the transverse mass of all final
state particles divided by 2. For the fragmentation scale, its central value
is set to the maximum of the transverse momentum of all final state particles
for inclusive hadron production, and to the transverse momentum of the jet
multiplied by the jet cone size for fragmentation inside the jet Kaufmann et
al. (2015). Theoretical uncertainties are included in the covariance matrix of
$\chi^{2}$ calculations, and are assumed to be fully correlated among points
in each subset of the data. These uncertainties are estimated by the half
width of the envelope of theoretical predictions of the 9 scale combinations
of $\mu_{F}/\mu_{F,0}=\mu_{R}/\mu_{R,0}=\\{1/2,1,2\\}$ and
$\mu_{D}/\mu_{D,0}=\\{1/2,1,2\\}$. The impact of different choices of the
various nominal scales is minimal as long as we include theoretical
uncertainties.

A detailed explanation of the experimental data sets used in the analysis is
described in the following. For $pp$ collisions, we initially incorporate
measurements on unidentified charged hadron production from fragmentation
inside jets by CMS and ATLAS Aaboud et al. (2019); Sirunyan et al. (2018); Aad
et al. (2021); Sirunyan et al. (2022); Aaboud et al. (2018); Aad et al. (2011,
2019). These measurements can be categorized into three groups, involving the
use of an isolated photon or a $Z$ boson recoiling against the fragmented
parton, or using the clustered jet as a reference of the fragmented parton. It
is assumed that the measured cross sections of unidentified charged hadrons
are a combination of charged pion, kaon and proton, while the residual
contribution from other charged hadrons is negligible and can be safely
disregarded. Similarly, LHCb has conducted separate measurements of
$\pi^{\pm}$, $K^{\pm}$ and $p/\bar{p}$ production from fragmentation inside
jets for $Z+$jet process at LHC 13 TeV Aaij et al. (2023). Additionally, we
include inclusive hadron production measurements from ALICE Abelev et al.
(2014); Acharya et al. (2020, 2021) and STAR Agakishiev et al. (2012). Only
ratios of production cross sections of different charged hadrons or of
different collision energies are considered to minimize the impact of
normalization uncertainties. In the analysis of SIA, we incorporate a
comprehensive set of data from TASSO, TPC below the $Z$-pole Braunschweig et
al. (1989); Aihara et al. (1988), as well as from OPAL, ALEPH, DELPHI, and SLD
at the $Z$-pole Akers et al. (1994); Buskulic et al. (1995); Abreu et al.
(1998); Abe et al. (2004), and from OPAL and DELPHI above the $Z$-pole
Abbiendi et al. (2003); Abreu et al. (2000). These measurements encompass the
production of $\pi^{\pm}$, $K^{\pm}$ and $p/\bar{p}$ separately, except for
OPAL at a collision energy of 202 GeV. For SIDIS, we utilize data on the total
rate and charge asymmetry of production of unidentified charged hadrons from
H1 and ZEUS at high $Q^{2}$ Aaron et al. (2007, 2009); Abramowicz et al.
(2010). There are also measurements on the production of identified charged
hadrons from COMPASS at relatively low $Q^{2}$ with isoscalar (06$I$) or
proton (16$p$) targets Adolph et al. (2017a, b); Pierre (2019). Only the data
with the highest inelasticity and the two subsets with largest Bjorken-$x$ are
included.

Strict selection criteria are applied to the kinematics of data points.
Specifically, we exclusively select data points corresponding to momentum
fractions $z>0.01$ at LO except for single inclusive hadron production in $pp$
collisions. Additionally, it is required that $p_{T,h}(E_{h})>$ 4 GeV for data
from $pp$ collisions (SIA and SIDIS), with the hadron energy being measured in
the Breit frame for SIDIS. Also note that for COMPASS data we did not apply
the $E_{h}$ cut since the $Q$ value can be as low as 3.7 GeV.

All the aforementioned data can be categorized into subsets based on the range
of jet $p_{T}$ for jet fragmentation, $Q^{2}$ for SIDIS, and collision energy
for inclusive hadron production at $pp$ collisions and SIA. In total, there
are 138 subsets. As an illustrative example, we present the $z$ coverage for
all subsets of the ATLAS jet fragmentation measurements after the kinematic
selections in Fig. 1. They include different $p_{T}$ bins from the inclusive
jet production at 7 and 5.02 TeV, as well as from both the central and forward
jet in dijet production at 13 TeV. For each $p_{T}$ bin, the left end of the
bar indicates the lower limit of $z$ of the subset, and the width of each
colored region represents the relative proportion of jet flavors, namely $u$,
$d$-quark, gluon and other partons, shown in linear scale. It is clearly shown
that the gluon fragmentation is dominant for jet production at low-$p_{T}$,
e.g., 200 GeV, and at small rapidity. This highlights the strong constraints
from jet productions on the gluon FFs.

Figure 1:  Coverage on the momentum fraction $z$ for all subsets of the ATLAS
jet fragmentation measurements Aaboud et al. (2018); Aad et al. (2011, 2019)
after the kinematic selections. The left end of the bar indicates the lower
limit of $z$ of the subset and the width of each colored region represents the
relative portion of jet flavors, including jets from $u$, $d$-quarks, gluon,
and other partons, shown in linear scale.

The fit.– The log-likelihood functions $\chi^{2}$ are calculated for each
subset using predictions from the prescribed theory and the covariance
matrices constructed from both experimental and theoretical uncertainties. A
best fit for the parameters of the FFs is determined by minimizing the total
$\chi^{2}$ using the MINUIT program James and Roos (1975). In Table. 1, we
provide a summary of various results demonstrating the quality of the best
fit. The global $\chi^{2}$ is 1231.5 units for a total of 1370 data points,
indicating overall good agreement between theory and data. The
$\chi^{2}/N_{pt}$ values are all below 1 for groups of data of $pp$ collision,
SIA and SIDIS. For individual experiments, only the ATLAS $Z/\gamma+$ jet and
the ZEUS unidentified charged hadron production measurements show slightly
worse agreement with $\chi^{2}/N_{pt}>2$.

Experiments | $N_{pt}$ | $\chi^{2}$ | $\chi^{2}/N_{pt}$
---|---|---|---
ATLAS jets † | 446 | 350.8 | 0.79
ATLAS $Z/\gamma$+jet † | 15 | 31.8 | 2.12
CMS $Z/\gamma$+jet † | 15 | 17.3 | 1.15
LHCb $Z$+jet | 20 | 30.6 | 1.53
ALICE inc. hadron | 147 | 150.6 | 1.02
STAR inc. hadron | 60 | 42.2 | 0.70
$pp$ sum | 703 | 623.3 | 0.89
TASSO | 8 | 7.0 | 0.88
TPC | 12 | 11.6 | 0.97
OPAL | 20 | 16.3 | 0.81
OPAL (202 GeV) † | 17 | 24.2 | 1.42
ALEPH | 42 | 31.4 | 0.75
DELPHI | 78 | 36.4 | 0.47
DELPHI (189 GeV) | 9 | 15.3 | 1.70
SLD | 198 | 211.6 | 1.07
SIA sum | 384 | 353.8 | 0.92
H1 † | 16 | 12.5 | 0.78
H1 (asy.) † | 14 | 12.2 | 0.87
ZEUS † | 32 | 65.5 | 2.05
COMPASS (06$I$) | 124 | 107.3 | 0.87
COMPASS (16$p$) | 97 | 56.8 | 0.59
SIDIS sum | 283 | 254.4 | 0.90
Global total | 1370 | 1231.5 | 0.90
Table 1:  The number of data points, $\chi^{2}$, and $\chi^{2}/N_{pt}$ for the
global data sets, groups of data from $pp$ collision, from SIA, and from
SIDIS. The values are also shown for individual experiments. Data sets for
production of unidentified charged hadrons are marked with a dagger.

It is important to assess the agreement for each of the 138 subsets in the
best fit. To account for variations in the number of data points, we utilize
an effective Gaussian variable $S_{E}(\chi^{2},N_{pt})$, given by the
equation:

$S_{E}=\frac{(18N_{pt})^{3/2}}{18N_{pt}+1}\left\\{{6\over
6-\ln(\chi^{2}/N_{pt})}-{9N_{pt}\over 9N_{pt}-1}\right\\},$ (2)

which follows a normal distribution if $N_{pt}$ is not too small Hou et al.
(2021). The histogram of $S_{E}$ from our best fit closely resembles a
Gaussian distribution, with a mean of -0.33 and a standard deviation of 1.43.
The majority of the subsets (132 out of 138) have $S_{E}$ values smaller than
2, indicating good agreement. The deviation of the histogram from the standard
normal distribution suggests a potential underestimation of experimental or
theoretical uncertainties by an average factor of 1.43. This motivates a
choice of tolerance of $\Delta\chi^{2}=1.43^{2}\approx 2$ in our determination
of uncertainties of the FFs using the Hessian method Hou et al. (2021). In
conclusion, we arrive at a best fit of the $\pi^{+}$, $K^{+}$ and $p$ FFs,
along with 126 Hessian error FFs, two for each of the eigenvector directions.
The uncertainties of any observables can be calculated using predictions from
the 127 FFs and the asymmetric Hessian formula Lai et al. (2010).

Resulting FFs.– Our newly obtained NLO FFs are compared to those from NNFF and
DSS for $u$, $d$, $s$-quark and gluon at $Q=5$ GeV in Fig. 2. For simplicity,
we only show the FFs with positively and negatively charged particles summed
over for a specific final state hadron. The DSS fits have a limit of 0.05 in
the momentum fraction, with $\pi^{\pm}$, $K^{\pm}$ and $p/\bar{p}$ results
from DSS21 Borsa et al. (2022), DSS17 de Florian et al. (2017) and DSS07 de
Florian et al. (2007), respectively.

The NNFF sets used are NNFF1.0 Bertone et al. (2017). The estimated
uncertainties of FFs are also shown for NNFFs and for our fit. Reasonable
agreement can be observed between our results and DSS for FFs of $u$ and $d$
quarks to $\pi^{\pm}$, and of $u$ quark to $K^{\pm}$. However, large
discrepancies are found for FFs to protons and for FFs of gluon to all three
charged hadrons. NNFFs generally exhibit larger uncertainties and can even
become negative in some kinematic regions. Note NNFFs assume exact flavor
symmetry in parametrizations of FFs of constituent quarks to $\pi^{+}$ and
$K^{+}$. In contrast, our results show an uncertainty of 3%, 4% and 8% for FFs
of gluon to $\pi^{\pm}$ at $z=0.05,\,0.1$ and $0.3$, respectively, which are
significantly improved compared with NNFFs. The uncertainties are about 4%, 4%
and 7% for FFs of $u$-quark to $\pi^{\pm}$, $K^{\pm}$ and $p/\bar{p}$ at
$z=0.3$. The high precision of gluon FFs is mostly due to the data of jet
fragmentation at the LHC, as well as the data of inclusive hadron production
in $pp$ collisions. Note that all aforementioned analyses, including ours,
utilize a zero-mass scheme for heavy quarks, namely parametrizing FFs for
heavy quarks at their mass thresholds.

Figure 2:  Comparison of our NLO fragmentation functions to those from NNFFs
and DSS at $Q=5$ GeV. The DSS fits have a limit of 0.05 in the momentum
fraction, with $\pi^{\pm}$, $K^{\pm}$ and $p/\bar{p}$ results from DSS21,
DSS17 and DSS07, respectively. The NNFF sets used are from NNFF1.0. The
estimated uncertainties of FFs are also shown for NNFFs and for our fit. For
$\pi^{\pm}$ and $p/\bar{p}$ the left (right) panel shows results for $u$ and
$d$ ($s$ and $g$). For $K^{\pm}$ the left (right) panel shows results for $u$
and $s$ ($d$ and $g$).

With the significant reduction of uncertainty of FFs from the simultaneous
global analysis, one should be able to check the fundamental sum rules of FFs,
including momentum sum rule, charge sum rule, and particle number sum rule.
These sum rules arise from the number density interpretation of FFs and have
been pointed out to be problematic by Rogers and Collins in Ref.Collins and
Rogers (2023). Therefore, it is of critical importance to check from a data-
driven analysis whether these fundamental properties of FFs are valid. That is
especially the case for the momentum sum rule,

$\sum_{h}\int_{0}^{1}dzzD_{i}^{h}(z,Q)=1,$ (3)

due to the suppression of small-$z$ contributions. The above momentum sum
rule, summed over all hadrons, was approved for each specific flavor of quark
and gluon and remains to be valid for standard renormalization scheme like
$\overline{\rm MS}$ as we used in our analysis. We first calculate the total
momentum fraction $\langle
z\rangle_{i}^{h}=\int_{z_{min}}^{1}dzzD_{i}^{h}(z,Q)$ of parton $i$ carried by
hadron $h$ for $z_{min}<z<1$ at initial scale $Q=5$ GeV, where the choice of
the lower limit $z_{min}$ is $0.01$ for $g$, $u$ and $d$ quarks, and $0.088$
for $s$, $c$ and $b$ quarks, based on the kinematic coverage of relevant data.
Fragmentation of the latter three quarks are only directly constrained by SIA
data at $Z$-pole. The final results of $\langle z\rangle_{i}^{h}$ for light
quarks and gluon are shown in Table. 2, where the central values and
uncertainties are calculated from our best-fit and Hessian error FFs. It shows
that the three charged hadrons carry approximately 53% to 50% of the momentum
of $u$, $d$ quarks and gluon. Our analysis reveals a preference for larger FFs
of $s$ quark to $\pi^{\pm}$, with each carrying about 16% of the total
momentum of $s$ quark, mainly due to the influence of SIA data. One possible
reason is because the SIA measurements on spectrum of $\pi^{\pm}$ also include
feed-down contributions from short-lived strange hadrons Abe et al. (2004).

The total momentum carried by different sets of hadrons as functions of
$z_{min}$ for $Q$=5 (left) and 100 GeV (right) are shown in Fig. 3. The
vertical dashed lines show the corresponding lower limit of $z$ as constrained
by experimental data for $g,~{}u,~{}d$ and $s$ quark fragmentation. One can
see from the top-left figure that $\langle z\rangle_{i=g,u,d}^{h}$ for light
charged hadrons reach to a saturation region within the current experimental
coverage in contrast to the strange quark. That can provide us a reliable test
of momentum sum rule for $g,~{}u$ and $d$ quarks if the FFs to neutral hadrons
are also determined at similar precision with future measurements. For now, as
an exploratory study, one can calculate ratios of energies carried by all
hadrons and by light charged hadrons as functions of $z_{min}$ using PYTHIA8
Bierlich et al. (2022) for $q\bar{q}$ and $gg$ production in $e^{+}e^{-}$
collisions with a center of mass energy of 200 GeV. Therefore, we can estimate
total momentum carried by all hadrons, as shown in the lower panel of Fig. 3,
by applying the scaling factors derived from PYTHIA8. The central values are
slightly lower than 1 for $u$, $d$ quarks and gluon when extrapolated into
small-$z_{min}$ region. That is consistent with the momentum sum rule
considering the shown uncertainty band and additional uncertainties from the
scaling factors not counted. For strange quark the values can be well above 1
due to both the ambiguity in experimental measurements mentioned earlier and
the limited kinematic coverage of data. We leave detailed investigations of
this anomaly for a future publication.

Figure 3:  Total momentum of the partons, including $g$, $u$, $d$, and $s$ quarks, carried by charged hadrons ($\pi^{\pm}$, $K^{\pm}$, $p$ and $\bar{p}$) in the fragmentations, as functions of $z_{min}$. The central values and uncertainties are calculated from our best fit and Hessian error FFs at $Q=5$ GeV and 100 GeV. The green (red) vertical lines indicate the kinematic coverage of relevant data of constraints for $g$, $d$, $u$ ($s$). In the lower panel, the results are rescaled by ratios of energies carried by all hadrons and by charged hadrons calculated from PYTHIA8 using $q\bar{q}$ and $gg$ production in $e^{+}e^{-}$ collisions with a center of mass energy of 200 GeV.  mom. | $g(z>0.01)$ | $u(z>0.01)$ | $d(z>0.01)$ | $s(z>0.088)$
---|---|---|---|---
$\pi^{+}$ | $0.200^{+0.008}_{-0.008}$ | $0.262^{+0.017}_{-0.016}$ | $0.128^{+0.020}_{-0.019}$ | $0.161^{+0.013}_{-0.013}$
$K^{+}$ | $0.018^{+0.004}_{-0.003}$ | $0.058^{+0.005}_{-0.004}$ | $0.019^{+0.004}_{-0.004}$ | $0.015^{+0.002}_{-0.002}$
$p$ | $0.035^{+0.006}_{-0.005}$ | $0.044^{+0.004}_{-0.004}$ | $0.022^{+0.002}_{-0.002}$ | $0.015^{+0.002}_{-0.002}$
$\pi^{-}$ | $0.200^{+0.008}_{-0.008}$ | $0.128^{+0.020}_{-0.019}$ | $0.299^{+0.054}_{-0.049}$ | $0.161^{+0.013}_{-0.013}$
$K^{-}$ | $0.018^{+0.004}_{-0.003}$ | $0.019^{+0.004}_{-0.004}$ | $0.019^{+0.004}_{-0.004}$ | $0.205^{+0.014}_{-0.013}$
$\bar{p}$ | $0.035^{+0.006}_{-0.005}$ | $0.019^{+0.003}_{-0.003}$ | $0.019^{+0.003}_{-0.003}$ | $0.015^{+0.002}_{-0.002}$
Sum | $0.507^{+0.014}_{-0.013}$ | $0.531^{+0.015}_{-0.013}$ | $0.506^{+0.042}_{-0.037}$ | $0.572^{+0.029}_{-0.028}$
Table 2:  Total momentum of the partons, including $g$, $u$, $d$ and $s$
quarks, carried by various charged hadrons ($\pi^{\pm}$, $K^{\pm}$, $p$ and
$\bar{p}$) in the fragmentations. The central values and uncertainties are
calculated from our best fit and Hessian error FFs at $Q=5$ GeV. The last row
is the sum of all charged hadrons.

Conclusions.– In summary, we present a joint determination of FFs for charged
hadrons from a global analysis at NLO in QCD, including estimations of
uncertainties. Our analysis demonstrates good agreement between our best-fit
predictions and various measurements in SIA, SIDIS and $pp$ collisions. Our
work introduces several advances including both a consistent simultaneous
fitting framework and new theoretical inputs. Notably, we have included
measurements on jet fragmentation at the LHC into the global analysis,
resulting in strong constraints on the gluon FFs. Comparing our results with
previous determinations, we find significant differences, especially in the
fragmentation to protons. Discrepancies are also observed for FFs of non-
constituent quarks and gluon to charged pions. Additionally, we provide
results on the total momentum of partons carried by various charged hadrons.
Our results pave the way for future precision studies of fragmentation and QCD
at the LHC and the upcoming electron-ion colliders.

Acknowledgments. The work of JG is supported by the National Natural Science
Foundation of China (NSFC) under Grant No. 11875189. HX is supported by the
NSFC under Grants Nos. 12022512, 12035007, by the Guangdong Major Project of
Basic and Applied Basic Research Nos. 2020B0301030008, 2022A1515010683. YZ is
supported by the NSFC under Grant No. U2032105. XS is supported by the
Helmholtz-OCPC Postdoctoral Exchange Program under Grant No. ZD2022004.

## References

  * Ablikim et al. (2023) M. Ablikim et al. (BESIII), Phys. Rev. Lett. 130, 231901 (2023), eprint 2211.11253.
  * Aaij et al. (2023) R. Aaij et al. (LHCb), Phys. Rev. D 108, L031103 (2023), eprint 2208.11691.
  * Collins et al. (1989) J. C. Collins, D. E. Soper, and G. F. Sterman, Adv. Ser. Direct. High Energy Phys. 5, 1 (1989), eprint hep-ph/0409313.
  * Metz and Vossen (2016) A. Metz and A. Vossen, Prog. Part. Nucl. Phys. 91, 136 (2016), eprint 1607.02521.
  * Accardi et al. (2016) A. Accardi et al., Eur. Phys. J. A 52, 268 (2016), eprint 1212.1701.
  * Anderle et al. (2021) D. P. Anderle et al., Front. Phys. (Beijing) 16, 64701 (2021), eprint 2102.09222.
  * Burke et al. (2014) K. M. Burke et al. (JET), Phys. Rev. C 90, 014909 (2014), eprint 1312.5003.
  * Collins and Soper (1982) J. C. Collins and D. E. Soper, Nucl. Phys. B 194, 445 (1982).
  * Collins and Rogers (2023) J. Collins and T. Rogers (2023), eprint 2309.03346.
  * Pitonyak et al. (2023) D. Pitonyak, C. Cocuzza, A. Metz, A. Prokudin, and N. Sato (2023), eprint 2305.11995.
  * de Florian et al. (2007) D. de Florian, R. Sassot, and M. Stratmann, Phys. Rev. D 76, 074033 (2007), eprint 0707.1506.
  * Hirai et al. (2007) M. Hirai, S. Kumano, T. H. Nagai, and K. Sudoh, Phys. Rev. D 75, 094009 (2007), eprint hep-ph/0702250.
  * Albino et al. (2008) S. Albino, B. A. Kniehl, and G. Kramer, Nucl. Phys. B 803, 42 (2008), eprint 0803.2768.
  * Bertone et al. (2018) V. Bertone, N. P. Hartland, E. R. Nocera, J. Rojo, and L. Rottoli (NNPDF), Eur. Phys. J. C 78, 651 (2018), eprint 1807.03310.
  * Khalek et al. (2021) R. A. Khalek, V. Bertone, and E. R. Nocera (MAP (Multi-dimensional Analyses of Partonic distributions)), Phys. Rev. D 104, 034007 (2021), eprint 2105.08725.
  * Moffat et al. (2021) E. Moffat, W. Melnitchouk, T. C. Rogers, and N. Sato (Jefferson Lab Angular Momentum (JAM)), Phys. Rev. D 104, 016015 (2021), eprint 2101.04664.
  * Stratmann and Vogelsang (1997) M. Stratmann and W. Vogelsang, Nucl. Phys. B 496, 41 (1997), eprint hep-ph/9612250.
  * Salam and Rojo (2009) G. P. Salam and J. Rojo, Comput. Phys. Commun. 180, 120 (2009), eprint 0804.3755.
  * Liu et al. (2023) C. Liu, X. Shen, B. Zhou, and J. Gao, JHEP 09, 108 (2023), eprint 2305.14620.
  * Dulat et al. (2016) S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump, and C. P. Yuan, Phys. Rev. D 93, 033006 (2016), eprint 1506.07443.
  * Zyla et al. (2020) P. A. Zyla et al. (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Kaufmann et al. (2015) T. Kaufmann, A. Mukherjee, and W. Vogelsang, Phys. Rev. D 92, 054015 (2015), [Erratum: Phys.Rev.D 101, 079901 (2020)], eprint 1506.01415.
  * Aaboud et al. (2019) M. Aaboud et al. (ATLAS), Phys. Rev. Lett. 123, 042001 (2019), eprint 1902.10007.
  * Sirunyan et al. (2018) A. M. Sirunyan et al. (CMS), Phys. Rev. Lett. 121, 242301 (2018), eprint 1801.04895.
  * Aad et al. (2021) G. Aad et al. (ATLAS), Phys. Rev. Lett. 126, 072301 (2021), eprint 2008.09811.
  * Sirunyan et al. (2022) A. M. Sirunyan et al. (CMS), Phys. Rev. Lett. 128, 122301 (2022), eprint 2103.04377.
  * Aaboud et al. (2018) M. Aaboud et al. (ATLAS), Phys. Rev. C 98, 024908 (2018), eprint 1805.05424.
  * Aad et al. (2011) G. Aad et al. (ATLAS), Eur. Phys. J. C 71, 1795 (2011), eprint 1109.5816.
  * Aad et al. (2019) G. Aad et al. (ATLAS), Phys. Rev. D 100, 052011 (2019), eprint 1906.09254.
  * Abelev et al. (2014) B. B. Abelev et al. (ALICE), Phys. Lett. B 736, 196 (2014), eprint 1401.1250.
  * Acharya et al. (2020) S. Acharya et al. (ALICE), Phys. Rev. C 101, 044907 (2020), eprint 1910.07678.
  * Acharya et al. (2021) S. Acharya et al. (ALICE), Eur. Phys. J. C 81, 256 (2021), eprint 2005.11120.
  * Agakishiev et al. (2012) G. Agakishiev et al. (STAR), Phys. Rev. Lett. 108, 072302 (2012), eprint 1110.0579.
  * Braunschweig et al. (1989) W. Braunschweig et al. (TASSO), Z. Phys. C 42, 189 (1989).
  * Aihara et al. (1988) H. Aihara et al. (TPC/Two Gamma), Phys. Rev. Lett. 61, 1263 (1988).
  * Akers et al. (1994) R. Akers et al. (OPAL), Z. Phys. C 63, 181 (1994).
  * Buskulic et al. (1995) D. Buskulic et al. (ALEPH), Z. Phys. C 66, 355 (1995).
  * Abreu et al. (1998) P. Abreu et al. (DELPHI), Eur. Phys. J. C 5, 585 (1998).
  * Abe et al. (2004) K. Abe et al. (SLD), Phys. Rev. D 69, 072003 (2004), eprint hep-ex/0310017.
  * Abbiendi et al. (2003) G. Abbiendi et al. (OPAL), Eur. Phys. J. C 27, 467 (2003), eprint hep-ex/0209048.
  * Abreu et al. (2000) P. Abreu et al. (DELPHI), Eur. Phys. J. C 18, 203 (2000), [Erratum: Eur.Phys.J.C 25, 493 (2002)], eprint hep-ex/0103031.
  * Aaron et al. (2007) F. D. Aaron et al. (H1), Phys. Lett. B 654, 148 (2007), eprint 0706.2456.
  * Aaron et al. (2009) F. D. Aaron et al. (H1), Phys. Lett. B 681, 125 (2009), eprint 0907.2666.
  * Abramowicz et al. (2010) H. Abramowicz et al. (ZEUS), JHEP 06, 009 (2010), [Erratum: JHEP 10, 030 (2010)], eprint 1001.4026.
  * Adolph et al. (2017a) C. Adolph et al. (COMPASS), Phys. Lett. B 764, 1 (2017a), eprint 1604.02695.
  * Adolph et al. (2017b) C. Adolph et al. (COMPASS), Phys. Lett. B 767, 133 (2017b), eprint 1608.06760.
  * Pierre (2019) N. Pierre, Ph.D. thesis, IRFU, Saclay (2019).
  * James and Roos (1975) F. James and M. Roos, Comput. Phys. Commun. 10, 343 (1975).
  * Hou et al. (2021) T.-J. Hou et al., Phys. Rev. D 103, 014013 (2021), eprint 1912.10053.
  * Lai et al. (2010) H.-L. Lai, M. Guzzi, J. Huston, Z. Li, P. M. Nadolsky, J. Pumplin, and C. P. Yuan, Phys. Rev. D 82, 074024 (2010), eprint 1007.2241.
  * Borsa et al. (2022) I. Borsa, D. de Florian, R. Sassot, and M. Stratmann, Phys. Rev. D 105, L031502 (2022), eprint 2110.14015.
  * de Florian et al. (2017) D. de Florian, M. Epele, R. J. Hernandez-Pinto, R. Sassot, and M. Stratmann, Phys. Rev. D 95, 094019 (2017), eprint 1702.06353.
  * Bertone et al. (2017) V. Bertone, S. Carrazza, N. P. Hartland, E. R. Nocera, and J. Rojo (NNPDF), Eur. Phys. J. C 77, 516 (2017), eprint 1706.07049.
  * Bierlich et al. (2022) C. Bierlich et al., SciPost Phys. Codeb. 2022, 8 (2022), eprint 2203.11601.