\section{Introduction}
\label{intro}

After a quench from the disordered phase above the critical temperature
$T_c$ to a final temperature $T<T_c$ ferromagnetic materials undergo
phase-ordering~\cite{PuriWad09,Bray94,CRP,Cor15,CorCugYos11}.
The system orders locally inside domains whose typical size $L(t)$ grows in time
until equilibration takes place on a timescale $\tau _{eq}$ which diverges with the system (linear) size ${\cal L}$.
This relaxation process is called coarsening and it is often accompanied by a dynamical scaling symmetry,
which amounts to the physical fact that configurations at different times are statistically similar
upon measuring distances in unit of $L(t)$. The latter usually increases algebraically,
$L(t)\sim t^{1/z}$, where $z$ is a non-equilibrium dynamical exponent that is unrelated to
any equilibrium property. This exponent is also independent of the quench temperature~\cite{Bray94,Bray90_rg}, a property
which is true for any universal quantity, because it can
be shown that temperature is an irrelevant parameter in the sense of the renormalisation group~\cite{Bray90_rg}.
This statement applies for quenches to $0<T<T_c$: indeed when cooling down to $T=T_c$
the process is qualitatively different because the order parameter vanishes in the target equilibrium state, at variance with
what happens when $T<T_c$. Quenches to $T=0$, on the other hand,
may also have peculiar properties because any activated process is forbidden.
To be concrete, let us discuss a system with a scalar non-conserved order parameter.

{\bf Short-range systems --} If interactions are restricted to nearest
neighboring (nn) spins, the Ising model with Glauber single spin-flip kinetics~\cite{Glauber1963}
represents an appropriate description.
Letting space dimension $d>1$ above the lower critical one, in order to have a finite $T_c$,
after quenching to $0<T<T_c$ ordered domains grow at late times with $z=z_{\rm cd}=2$ until the system eventually
attains the equilibrium state in a time $\tau_{eq}$ which
is in most cases $\tau _{eq} \sim {\cal L}^z$. We indicate with $z_{\rm cd}$ this value of $z$ because the
motion of interface in this case is curvature driven~\cite{AllenCahn79,Bray94}.

For quenches to $T=0$, even if the motion of interfaces is no lomger curvature driven but has a diffusive character,
one still observes a growth law with an exponent $z=z_{\rm diff}=z_{\rm cd}=2$.
Another difference between $T=0$ and $T>0$ is that the equilibration in the former quench
may be impeded due to blocking of the system into infinitely lived
metastable states~\cite{BarKraRed09,BlaCorCugPic14}. In $d=2$ these blocked states are stripes
with flat interfaces extending along the lattice directions, as the one denoted with the letter $a$ in Fig.\ref{stab_interf}.
Clearly, the fact that such
flat interfaces are frozen does not only affect the fate of the system in a quench
to $T=0$, but also the preceding dynamics. Indeed it was shown in~\cite{CorLipZan08}
that, although the value of $z$ does not change in going from $T>0$ to $T=0$, since $z_{\rm diff}=z_{\rm cd}$,
some other non-equilibrium exponents related to the geometry of interfaces, do change.

Furthermore,
even if the quench is made to a finite $T$, the different dynamics associated to $T=0$ is observed
in a pre-asymptotic regime which can be rather long if $T$ is small enough.
Similarly, for sufficiently low $T>0$, although the metastable states are eventually escaped, this happens on
timescales that can be huge, greatly delaying the equilibration with respect to what happens at
higher $T$ where, as already mentioned, $\tau_{eq}\sim {\cal L}^z$.
Summarising, the case with $T=0$ displays peculiar features which can strongly affect even the dynamics
of quenches to finite $T$, at least pre-asymptotically.

A separated discussion is deserved by the case $d=1$, because here $T_c=0$, meaning that, strictly speaking $T=0$ is the only possibility to cool the system to a magnetised state. However, even for quenches to a finite $T$, a coarsening stage takes
place because domains keep growing as in a $T=0$ quench until their size $L(t)$ reaches the equilibrium
correlation length $\xi$ or the system size ${\cal L}$, after which equilibration occurs. Also in this case the motion of
interfaces has a diffusive character and one finds $z=z_{\rm diff}=2$~\cite{Glauber1963}
but, at variance with higher dimensions, there are no metastable states due to the trivial $1d$ lattice
geometry.

{\bf Long-range systems --} The presence of long-range interactions changes a lot the above picture.
The equilibrium scenario
with a coupling between spins at distance $r$ decaying as $r^{-(d+\sigma)}$
is well understood. There is an ordered phase
below a finite $T_c$ also in $d=1$ provided, in this case, that $\sigma \le 1$.
Furthermore, the energy is extensive and the system is additive for any $d$ if $\sigma >0$, a regime sometimes
denoted as weak long-range, whereas extensivity and additivity are lost for $\sigma \le 0$, the strong long-range case.
In this paper we will only focus on the weak long-range case $\sigma >0$.

Regarding the non-equilibrium properties, in quenches to a finite $T$ one finds~\cite{BrayRut94,CMJ19,PhysRevLett.125.180601}
$z=z_{\rm adv}=1+\sigma$ for $\sigma \le1$ and  $z=2$ for $\sigma >1$ (with logarithmic corrections right at $\sigma =1$).
This applies down to $d=1$ where, if $\sigma >1$, one has $T_c=0$ and the discussion made above for short-range
interactions regarding equilibration applies.
The different behaviour of the exponent $z$ in crossing $\sigma=1$ can be ascribed to qualitatively different underlying dynamical mechanisms.
When, for $\sigma >1$, $z$ takes the same value $z=2$ as in the case with nn interactions, the motion of interfaces
behaves as for nn interactions, namely it is diffusive for $d=1$ and is governed by the
curvature~\cite{BrayRut94,RutBray94} for $d>1$.
Instead, when $z=z_{\rm adv}$, i.e. for $\sigma \le 1$, the motion
of domains walls is advected
by the drift due to the long-range interactions between far away spins. Notice that
one has $z\to z_{ball}=1$ for $\sigma \to 0^+$. This case corresponds to an advection of interfaces so strong
to produce a completely deterministic motion and hence a ballistic regime. This is related
to the crossover from weak long-range to strong long-range occurring right at $\sigma =0$.

Regarding the differences between $T>0$ and $T=0$ when long-range interactions are present, the
situation is not clear at all. Presently the matter has been well understood only in
$d=1$~\cite{CLP_review,CLP_lambda,PhysRevE.102.020102}.
In this case, when quenching to $T=0$ one finds coarsening with $z=z_{ball}=1$ for any $\sigma$. Let
us recall that, instead, for $T>0$ one has $z\to 1$ only in the limit
$\sigma \to 0^+$. The interpretation is the following: when $T>0$ thermal fluctuation randomise the displacement
of interfaces and, therefore, the motion is not fully deterministic although still advected. In this case, in order
to have full determinism,  one has to go to the strong interaction limit $\sigma \to 0$. Instead, when $T=0$,
even a relatively small drift, present for any $\sigma \le 1$, leads to a deterministic motion with $z=z_{ball}=1$.
Also in this case, the asymptotic ballistic growth which sets in at $T=0$ is observed
as a pre-asymptotic behaviour after quenches to finite $T$ for any $\sigma $.
Finally, let us mention that also with long-range interactions, as in the nn case, metastability is not observed
in $d=1$.

In this paper we take a first step in the direction of understanding the ordering kinetics after quenches to
$T=0$ of systems with long-range interactions in $d>1$. The matter, which is largely unexplored, is relevant because in this case, as we shall
illustrate, the dynamical mechanism is different from the ones discussed above,
thus producing a new value $z={\cal Z}$ of the growth exponent. This is independent of $\sigma$ and
characterises also the
pre-asymptotic evolution in deep quenches to a finite $T$.

The origin of this new growth mechanism can be traced back to metastable configurations.
In order to discuss this point we invite the reader to refer to Fig.~\ref{stab_interf} as a preliminary illustration of the
various shapes of the interfaces, leaving further details contained in this figure, described in the detailed caption,
to the specific discussion that will be conducted in Sec.~\ref{sec.stab}.
As discussed previously, in $d>1$ with nn interactions a flat portion of an interface along a lattice direction
(Fig.~\ref{stab_interf}, letter $a$)
is stable at $T=0$.
However, interfaces are never perfectly flat neither aligned along the lattice directions in the coarsening stage, and hence they can be deformed
by flipping spins on the edges (Fig.~\ref{stab_interf}, letters $b,c,d$). This process typically occurs without any drift and
hence has a diffusive character~\cite{CorLipZan08}, a fact that leads to an exponent $z=z_{\rm diff}=2$.
The addition of long-range interactions changes the situation in two fundamental respects. The first modification
is an increase
of (meta)stable configurations: any globally straight interface (e.g. the one of Fig.~\ref{stab_interf}, letter $b$, where
all steps have the same size)
is now locally stable, even if its direction does not fit the orientation
of the underlying lattice.
This means that also spins on edges can be blocked, at variance with the nn case.
Secondly, if we perturb the constant slope interface (e.g. see Fig.~\ref{stab_interf}, letters $c,d$, where steps
have different lengths),
spins which are now free to flip are subjected to a deterministic drift, similarly to what happens in $d=1$
in the ballistic regime. These two new ingredients have, so to say, contrasting effects. On the one hand the presence
of the drift tends to speed the kinetics with respect to a diffusive case, leading to $z<2$;
on the other hand the drift is contrasted by the abundance of blocked states, causing $z>1$.
The result is a non trivial value, $z={\cal Z}$, with $1<{\cal Z}<2$. Even if this regime is observed unambiguously
in numerical simulation, a precise determination of ${\cal Z}$ is difficult with this technique also because
a pre-asymptotic stage is present. A direct analytical attack of the problem
seems to be difficult as well. This is due to the presence of long-range correlations which mix up with
lattice effects making calculations elusive.
Indeed, as we will discuss further, an approach on the continuum, i.e. off-lattice, which accurately
captures the asymptotic values $z_{\rm adv}, z_{\rm diff}$ in quenches to finite $T$, provides
a wrong value of ${\cal Z}$.
Considering instead a simplified system with a single domain, we are able to obtain a full analytic understanding
of the early stage where ${\cal Z}=3/2$ and a quite clean numerical determination of its asymptotic one with ${\cal Z}=4/3$. While the former value can be traced back to a relatively simple circular geometry of the growing domains, the latter originates from a more peculiar one caused by the
long-range interactions extending along the interfaces.

This paper is organised as follows: In Sec.~\ref{model} we introduce the kinetic model we study, discuss the
implementation of the numerical simulations and define the basic observable quantities we compute.
Sec.~\ref{numer_res} presents the results of numerical simulations of the model and the determination
of the exponent ${\cal Z}$. In Sec.~\ref{simplified} we simplify the problem by considering the evolution
of a single domain. This is studied both numerically and analytically, allowing us to
determine ${\cal Z}$. Finally, in Sec.~\ref{concl}, we conclude the paper by summarising the main results and
discuss some open points. In Appendix~\ref{Ewald}, we give some details of the numerical technique, namely Ewald summation, to incorporate long-range interactions in the system.

\section{The model and its numerical simulation}
\label{model}

We consider an Ising model with the Hamiltonian
\begin{equation}
{\cal H}(\{s_i\})=-\frac{1}{2} \sum _i \sum_j{}^{'} J(r) s_i s_j,
\label{ham}
\end{equation}
where $s_i=\pm 1$ are boolean spin variables on the sites $i$ of a two-dimensional square lattice
of liner size ${\cal L}$, whose spacing we assume to be unitary. The sum over $j$ is primed to indicate that the terms with $i = j$ are excluded, and
\begin{equation}
J(r)=r^{-(2+\sigma)}
\end{equation}
is a ferromagnetic coupling constant with $r=\vert \vec r_i - \vec r_j \vert$; $\vec r_{i}$ being the coordinate of site $i$ of the lattice. The usual nn Ising model arises on setting $J(r)=\delta _{r,1}$.

At equilibrium the long-range model has a para-ferromagnetic phase transition at a finite critical temperature $T_c(\sigma)$.
Critical exponents
match~\cite{PhysRevB.1.2265,Fisher1972critical,PhysRevB.8.281,PhysRevB.56.8945,PhysRevLett.89.025703,PhysRevE.95.012143} those of the corresponding nn model for sufficiently large values of $\sigma$,
i.e. $\sigma \ge \sigma_{nn}$, where $\sigma _{nn}$ has been estimated~\cite{PhysRevB.8.281,PhysRevE.95.012143} to be $\sigma_{nn}=7/4$,
whereas they match the mean field values for $\sigma <1$. In between, for $1< \sigma < \sigma_{nn}$, critical
exponents depend continuously on $\sigma$.

A kinetics is introduced by flipping single spins with Metropolis transition rates
\begin{equation}
w(s_i\to -s_i)={\cal L}^{-2}\min \left\{1,e^{-\frac{\Delta E}{T}}\right \},
\label{metrop}
\end{equation}
where $\Delta E$ is the change in energy due to the spin flip to be attempted, and we have set to unity
the Boltzmann constant. Time is measured in Monte Carlo steps (MCS), each of which corresponds to
${\cal L}^2$ elementary spin-flip attempts.
We consider a quenching protocol where the system is initialised in a configuration sorted from an infinite
temperature equilibrium ensemble, i.e. spins are randomly and independently set to $+1$ or $-1$.
This initial state is evolved at the quench temperature $T$ by means of the transition probabilities~(\ref{metrop})
with $T<T_c(\sigma)$.

Since the spin-spin interaction is long-ranged, the calculation of $\Delta E$ at every spin-flip attempt via Metropolis probability~(\ref{metrop}) is computationally expensive. To speed up the computer simulation, we store the local field for each spin at the beginning of the simulation. In this way we need to update it only once the spin-flip is accepted. This simple trick significantly speeds up the numerical calculations at any given quench temperature \cite{hmu1995,CMJ19}.

The main issue during the simulations of the system with Hamiltonian~(\ref{ham}) is the strong finite-size effect arising due to the long-range character of $J(r)$. One obvious way to diminish these effects is to use periodic boundary conditions via minimum-image convention \cite{fs2002}. In this approach, the square lattice is mapped onto a torus where each spin in the system interacts with other spins up to a certain cut-off distance ($\le {\cal L}/2$). When the interaction decays slowly, this natural cut-off limit on the interaction-range causes artefacts in the simulation results. Therefore, we need a more sophisticated approach to implement periodic boundary conditions in such systems. For this, we envision an infinite $2d$-lattice partitioned into infinite imaginary copies of the original simulation lattice. The central cell of this infinite lattice is the simulation lattice itself, and the imaginary copies, called images, lie across its periodic boundaries in both $x$ and $y$-direction. Implementing periodic boundary conditions with infinite images remove cut-off errors in the simulation results of long-range interacting systems. The effective interaction between two spins inside the simulation lattice can now be expressed as an infinite summation over all images
\begin{equation}
\label{jsum}
J(s_i,s_j) = \sum_{n_x} \sum_{n_y} \frac{1}{\vert \vec n + \vec r_{i} - \vec r_{j} \vert ^{2+\sigma}},
\end{equation}
where displacement vector $\vec n = (n_x,n_y)$ with $n_x, n_y = 0, \pm {\cal L}, \pm 2{\cal L}, ...$, representing the coordinates of image systems and the simulation lattice is located at $\vec n = (0,0)$. The infinite summation involved in (\ref{jsum}) has slow convergence in coordinate space; therefore, it is difficult to handle it directly during the numerical simulations. We have adapted the Ewald summation technique \cite{fs2002,ew1921,PhysRevE.95.012143}, which uses a clever trick to split it into two independent rapidly convergent summations, one in coordinate space and another in reciprocal space (see Appendix~\ref{Ewald}).

The main observable we consider in this paper is the characteristic size of the growing domains,
which we extract from the time dependent spin configurations by means of the equal time
correlation function
\begin{equation}
C(r,t)=\langle s_i(t)s_{i+\vec r}(t)\rangle,
\label{defcor}
\end{equation}
where $\langle \dots \rangle$ is an off-equilibrium average, namely taken over different initial conditions
and thermal histories. If dynamical scaling holds, in quenches below $T_c$ this quantity depends
on a single variable as~\cite{PuriWad09,Bray94,CRP,Cor15,CorCugYos11}
\begin{equation}
C(r,t)=f\left (\frac{r}{L(t)}\right ),
\label{cscal}
\end{equation}
where $f$ is a scaling function and $L(t)$ has the meaning of the size of the growing domains at time $t$.
This quantity can be extracted from the correlation function itself as \cite{PhysRevE.75.011113,PhysRevE.74.041106,CMJ19}
\begin{equation}
C\left (r=L(t),t\right )= \frac{C(0,t)}{2} \equiv \frac{1}{2} .
\end{equation}

The growth law after a quench to $T>0$ was predicted in~\cite{BrayRut94} by Bray and Rutenberg to be characterised by $z=z_{\rm cd}=2$ for $\sigma >1$ and by $z=z_{\rm adv}=1+\sigma$
for $\sigma \le 1$, with logarithmic corrections at $\sigma =1$. Such a prediction was obtained by using a
continuum model, based on a Ginzburg-Landau free energy, assuming a dynamical scaling symmetry and
resorting to a an energy scaling argument. These results have  been recently confirmed by Christiansen et al. \cite{CMJ19} by means of numerical simulations of the model considered in the present paper.

In order to study the growth law $L(t)$ it is useful to consider the effective exponent $z_{\rm eff}(t)$
defined by $1/z_{\rm eff}(t)=d\ln L(t)/d\ln t$.
Since, in order to speed our simulations, observable quantities are computed
only at discrete times $t_i$ (equally spaced in $\ln t$), the effective exponent is computed as follows
\begin{equation}
1/z_{\rm eff}(t_i) = \frac{\ln L(t_i)-\ln L (t_{i-1})}{\ln t_i-\ln t_{i-1}}.
\label{eqeffexp}
\end{equation}
The effective exponent computed by means of the definition~(\ref{eqeffexp}) may have
a noisy character, therefore we will also use the determination based on the symmetric derivative, which
amounts to
\begin{equation}
1/z_{\rm eff}^{(s)}(t_i)=\frac{ \ln L(t_{i+1})- \ln L(t_{i-1}) }{ \ln(t_{i+1}) - \ln(t_{i-1}) }.
\label{eqeffexps}
\end{equation}

\section{Numerical results}
\label{numer_res}

Let us now discuss the outcomes of our numerical simulations. First,  as a benchmark, we show
in Fig.~\ref{asymptoticRB} the results for quenches to relatively high $T$, where we expect
to recover similar results to those found by Christiansen et al.~\cite{CMJ19}, thus confirming the
Bray-Rutenberg asymptotic exponents. In this figure we see that an algebraic growth of $L(t)$
sets in at times $t\gtrsim 10$, for all $\sigma$ values. The value of $z$ matches with
the prediction of Bray and Rutenberg, as it can be appreciated in the main figure and in the inset,
where the effective exponent $1/z_{\rm eff}$ is plotted against $L(t)$.

Let us now move to the core results of this article, namely the effect of quenching to $T=0$.
In Fig.~\ref{preasympt} we summarise the behaviour of $L(t)$ not only in zero temperature quenches,
but also in deep quenches to small but finite $T$, in order to appreciate the existence of pre-asymptotic
effects. Let us start to discuss the case with $\sigma =3$ (lower right panel).
Since this value of $\sigma $ is rather large one could
naively expect to see a behaviour akin to that of the nn model. As we know, such a naive argument is correct
when applied to the shallow quenches discussed above, since for any $\sigma >1$ one recovers the
exponent $z_{\rm cd}=2$ of the nn case. Instead here one observe that for $T=0$ (black curve), for
$t\gtrsim 10$ one has an algebraic increase of $L(t)$ with an exponent $1/z\simeq 0.7$, definitely
larger than the one, $1/z_{\rm diff}=0.5$, of the corresponding quench to $T=0$ in the model with nn interactions
(and of the one $1/z_{\rm cd}=0.5$ predicted by Bray-Rutenberg for quenches to $T>0$ with $\sigma >1$).
Although a precise determination of $z$ is not possible from these simulations, this is surely
enough to establish the existence of a new exponent, that we denote with ${\cal Z}$, associated to the
zero temperature quenches with long  range interactions.
Notice the decrease of such exponent
at large $L(t)$ (meaning very large times) must be attributed to finite size effects start to be
appreciated.

The same panel of Fig.~\ref{preasympt} shows that a quench to
$T=10^{-4}$ behaves very similarly to the case with $T=0$, whereas a quench to $T=10^{-2}$
does so only up to times smaller than $t\simeq 10^2$, after which $L(t)$ slows gradually down until being
compatible with $1/z_{\rm eff}=1/z_{\rm cd}=1/2$. This pattern of behaviours can be interpreted as a crossover
occurring at $L_{\rm cross}(T,\sigma)$ between an early regime with $z={\cal Z}$ and a late stage with the
Bray-Rutenberg exponent. Since $L_{\rm cross}(T,\sigma)$ is a decreasing function of $T$ the crossover
cannot be observed in the quench to $T=10^{-4}$ because it occurs after the longest simulated times.

Moving now to the other values of $\sigma $ considered in Fig.~\ref{preasympt} we observe a similar
pattern of $L(t)$, with a regime characterised by a value
of $1/z$ of order $0.7$. This value, with the precision of our numerical simulations, appears to be
roughly independent of $\sigma$, a fact that suggests that ${\cal Z}$ is a universal exponent.
This can be better appreciated in the top panel of Fig.~\ref{summary_exp}, where we compare $1/z_{\rm eff}$
and $1/z_{\rm eff}^{(s)}$ in zero
temperature quenches with various $\sigma$. This figure shows that $1/z_{\rm eff}$ sets to a value of order
$1/z_{\rm eff}\simeq 2/3$ for $L(t)$ around $L(t)=10$, and then slowly increases towards a value that
exceeds $0.7$, before bending down due to finite size effects. A similar behaviour,
somewhat less noisy, is displayed by $1/z_{\rm eff}^{(s)}$. This observation, together
with the study of single domain models that we will discuss in Section~\ref{simplified}, leads us to the
conjecture that the exponent ${\cal Z}$ toggles between a pre-asymptotic value ${\cal Z}=3/2$ and an
asymptotic one ${\cal Z}=4/3$. These two values are indicated by straight lines in Figs.~\ref{preasympt}
and \ref{summary_exp}.

Comparing the various panels of Fig.~\ref{preasympt}
one can also be convinced that the crossover length $L_{\rm cross}(T,\sigma)$ is a (rather strongly) decreasing
function of $\sigma$. Indeed, for instance, with $\sigma =0.9$ (lower left panel) one has to rise $T$
by a factor $10^2$ (i.e. to set $T=1$) in order to see a crossover pattern similar to the one observed
at $\sigma =3$ (at $T=10^{-2}$).
Moreover, one observes that finite-size effects are more important for the smaller values of $\sigma $
because the longest range of the interactions makes the system feels the periodic boundary conditions  earlier.

Finally, let us discuss the issue of dynamical scaling. This was shown to be obeyed~\cite{CMJ19} in the
Bray-Rutenberg regime occurring in quenches to a finite $T$, by checking the form~(\ref{cscal}) of
the correlation function. In the following we study the same matter in the zero-temperature
quenches. Our results for $C(r,t)$ are shown in Fig.~\ref{corr_log}. The data for different times show
a very good collapse when plotted against the rescaled variable $r/L(t)$, as expected after Eq.~(\ref{cscal}) in
the presence of dynamical symmetry. This is true for all the values of $\sigma $ considered (we show only
a couple of them in Fig.~\ref{corr_log} but similar results are found for the other values).
For $\sigma =0.8$ (left panel), the curve at the longest time $t=500$ starts departing from the
collapse, but this is due to the onset of finite-size effects which can also be detected at these times
from the behaviour of $L(t)$ shown in Fig.~\ref{preasympt}. Hence, we conclude that the dynamical scaling
symmetry is at work also in zero temperature quenches, a fact that will be used in the following. Let us also
mention that metastable states, which are very important with nn interactions, here are absent (or at least greatly
suppressed), as found also in~\cite{Christiansen_arxiv}.

\section{Evolution of a single domain}
\label{simplified}

In the previous section we have shown that, with long-range interactions, one observes a fast growth
regime at $T=0$ regulated by a new exponent $1/{\cal Z}\simeq 0.7$. However, in the absence of
analytic approaches or of a detailed comprehension of the mechanism at work a precise determination
of such exponent on the basis of the sole simulations could not be obtained.

In this section we study the simpler process with a single domain. These approach has been shown
to provide a correct determination of the growth exponent $z$ in systems with nn interactions,
for any $d$~\cite{CorLipZan08,CLP_epl,libro}, and also with long-range interactions in $d=1$~\cite{CLP_review}.
We consider the same model of Sec.~\ref{model} with the only difference that
the initial state is not sorted from an infinite temperature ensemble but is built by hand
with a single domain of size $R$. This configuration is evolved with the
flip probabilities~(\ref{metrop}), which lead to a shrinkage and eventual disappearance
of the bubble. We start with a circular shape of radius $R$ (more precisely,
the lattice approximation to a circle, as in Fig.~\ref{pict_domain}), because we observed that,
starting from different shapes, a roughly circular one tends to be formed during the evolution
after a transient (however, the true asymptotic shape is not perfectly circular, see discussion in Sec.~\ref{Z_shape}).
This is shown in Fig.~\ref{bubble_geom} for an initially square bubble. Notice that the late
geometry of the bubble almost isotropic, at variance with what observed with nn interaction
during coarsening~\cite{PhysRevE.54.R2181} or in metastable states~\cite{GUNTHER1994194}.
The idea behind this single domain approach is based on the dynamical scaling hypothesis which,
in simple words, means that at any time $t$ in a coarsening system there is only a relevant
length $L(t)$. According to this assumption
the average time $t(L)$ needed to close a domain of initial size $L$ in such system scales as $L^z$.
Here we adopt the stratagem is to compute $t(L)$
with a single domain. Since $R$ is the initial size of such domain we indicate
with $t(R)$ this quantity and use $t(R)\sim R^z$ as a proxy to $t(L)$.
This, of course, implies another assumption, namely that the shrinkage of a domain
is not significantly influenced by the presence of the others, so that a single domain configuration suffices.
Previous studies~\cite{CorLipZan08,CLP_epl,CLP_review,libro} proved that such an approach works well
because as long as $\sigma>0$ the interaction is integrable, i.e. $\sum_{\vec r} J(r)$ is finite.
This means that the interaction between spins at distance larger than $L$, in the asymptotic stage
when $L$ has grown large, is negligible.

\subsection{Numerical results}

To start with, we have implemented numerical simulations of the bubble shrinkage process and we have
computed $t(R)$ by averaging over many dynamical trajectories. The results are contained in Fig.~\ref{bubble},
where in each panel we plot $t(R)$ (on the $x$ axis) vs $R$ (on the $y$ axis).
This choice allows us a direct comparison with Fig.~\ref{preasympt} where we plotted $L(t)$ versus $t$ for an
extended multi-domain system. For the same reason we also plot
the effective exponent, defined as previously in Eq.~(\ref{eqeffexp}) (with the replacement $L\to R$), in the inset.

Let us start with the case $\sigma =3$ (lower right panel). Working at $T=0$ one sees that a very clean value
$1/z_{\rm eff}=0.75$ sets in for large $R$, providing ${\cal Z}=4/3$. Instead, for very small sizes the effective
exponent has a zig-zag behaviour. This is perhaps due to the fact that for such small sizes some fine geometrical
details of the initial state become relevant. We noticed indeed that the most pronounced peaks (local maxima)
of $1/z_{\rm eff}$ correspond to values of $R=n^2$, where $n$ is an integer.
As we will further discuss below, as shown with Eq.~(\ref{topterr}), when the bubble diameter is a perfect
square number the largest terrace $\ell _L$ of the domain (the horizontal segment
denoted by 1 in Fig.~\ref{pict_domain})
is naturally an integer, therefore determining a discontinuity in the dynamical process.
This effect clearly reduces and tends to disappear with increasing $R$.
If we neglect these {\it special} maxima, for small $R$ the effective
exponent hits twice the value $1/z_{\rm eff}=2/3$. This might
suggest that ${\cal Z}$ toggles between a pre-asymptotic value ${\cal Z}=3/2$ and an asymptotic one
${\cal Z}=4/3$. This conjecture is going to be further supported.
The asymptotic value $4/3$ is assumed for $R\gtrsim 10^2$. Comparing with
the coarsening data of Fig.~\ref{preasympt} we can conclude that domains of such sizes are only reached
at the very end of the coarsening simulation when finite size effects start to set in. This would explain why in the
coarsening simulations one mostly observe an exponent in between ${\cal Z}=3/2$ and ${\cal Z}=4/3$
and only a hint of the convergence to the late one ${\cal Z}=4/3$ can be observed. Looking
now still in Fig.~\ref{bubble} for $\sigma =3$, but focusing on the data at the finite temperature
$T=0.5$, one sees that the effective exponent attains the value of Bray and Rutenberg, $1/z=0.5$.
This result attests the correctness of the single domain configuration to extract the exponent $z$.

We can now pass to the other panels of Fig.~\ref{bubble}, corresponding to smaller values of $\sigma$.
At finite $T$ one concludes that there is a
crossover phenomenon at a value $R_{\rm cross}(T,\sigma)$, where the exponent $z$
changes from ${\cal Z}$ to the Bray and Rutenberg value (for $\sigma =0.6$, $R_{\rm cross}(T,\sigma)$ is
probably too large to observe the crossover).
At $T=0$
the same value ${\cal Z}=4/3$ is neatly observed also for $\sigma =1.5$ and $2$. For values of $\sigma $ smaller than one
the determination of such exponent turns out to be less precise and there is the tendency to observe
larger values of $1/z_{\rm eff}$ upon decreasing $\sigma$.  This could be due to a stability effect
that we will discuss below. The data for the effective exponent at $T=0$ are summarised in the lower panel
of Fig.~\ref{summary_exp}.

\subsection{Stability of the bubble}
\label{sec.stab}

After having presented the results of the simulations of the bubble shrinkage, which provide a determination
of the growth law in a semi-quantitative agreement with what observed in the full coarsening model, we turn now to a study
of the microscopic kinetic mechanisms producing the closure process, in order to gain some better understanding.
We start by discussing the stability properties of an interface between two regions with differently aligned spins.
The zero-$T$ stability of a (positive, let's say) spin depends on the total field $h$ acting on it.
If $h>0$ the spin is stable (meaning that it cannot flip), if $h<0$ the spin is unstable, if $h=0$ the spin can also flip
but this usually corresponds to a neutral equilibrium, typical of the nn model (see below) but
not of the long-range one, because
a perfect compensation of extended interactions leading to $h=0$ is almost impossible.

In Fig.~\ref{stab_interf} we consider different types of discrete interfaces and
we refer the reader to its detailed caption.
A straight interface parallel to a lattice axis, see $(a)$, is clearly stable because all spins parallel
to a given interface spin, denoted with a $\bullet$, block its reversal while all other spins, above and below such line,
perfectly compensate. This is true for long-range coupling but also for nn coupling:
in the former case all full grey circles block $\bullet$, in the latter case only its two neighbouring ones.
When the interface contains kinks, either keeping a constant slope (i.e. all the steps have the same length) as in $(b)$ or acquiring a finite
curvature as in $(c,d)$, things are different and more complicated.
In the nn case kinks do not interact and each of them can move in both directions even at $T=0$.
With long-range interactions the picture is completely different and we must distinguish between the two
sides of the kink. Let us always focus on the $\bullet$ spins. If the interface has a constant slope it is stable,
this is shown in $(b)$ (see discussion in the caption). Notice that this is true irrespectively of the
value of $\sigma$.
If the interface has a finite curvature as in $(c)$, $\bullet$ spin is unstable while the spin
at its right is stable, meaning that the kink can move uphill (to the left) only. This also occurs
for any $\sigma$.
In $(d)$ the interface has a constant slope locally around $\bullet$, i.e. terraces locally have the same length,
but further uphill and downhill terraces are longer and shorter respectively.
In this case there are short-distance stabilising interactions (due to parallel uncompensated spins)
and longer-distance destabilising interactions (due to anti-parallel uncompensated spins):
the global effect depends on the details of the interface and on the value of $\sigma$:
it is likely that stability prevails for large $\sigma$ and instability
prevails for small $\sigma$. Finally, in $(e)$ it is shown that also the edge spins of a top terrace are unstable.

Things are even more complicated when considering interface spins not at a kink.
As a matter of fact in the limit of small $\sigma$ even bulk spins can flip.
In particular, considering a circular domain of size $R$, the central spin is unstable up to a critical size
$R=R_c\sim 2^{1/\sigma}$. Indeed the spin can't flip when the stabilising interaction
$\sim \int _{r<R}d\vec r \,r^{-(2+\sigma)}$ prevails over the destabilising one
$\sim \int _{r>R}d\vec r \,r^{-(2+\sigma)}$ produced by the anti-aligned spins outside the domain,
which gives the above mentioned result. Notice that in our bubble shrinkage simulations, in order to
have a reliable asymptotic (i.e. large $R$) determination of $t(R)$ one has to consider $R\gg R_c$,
otherwise there is a correction lowering $t(R)$ (because the domain shrinks to zero almost immediately
as soon as $R$ crosses $R_c$). Since for $\sigma \to 0$ it is $R_c\to \infty$, this is the possible explanation
of the observation made, regarding Fig.~\ref{bubble}, that decreasing $\sigma $ one observes values of
$1/{\cal Z}$ larger than $3/4$. In the next Section~\ref{Z_shape} we will clearly show the role of
bulk spin flips.

We conclude this part by stressing that our results are strongly related to the discreteness of the lattice.
In a continuum medium any surface spin of a compact bubble would feel a field favouring the closure of the bubble itself
see Fig.~\ref{continuum}, because the field produced by spins within the bubble is always
compensated by the field due to the mirror domain. Uncompensated spins are anti-parallel to bubble spins
therefore inducing its closure. At zero temperature this uniform driving force would give
a trivial dynamical exponent, $z=1$.

\subsection{Simplified dynamics and determination of ${\cal Z}$}
\label{Z_shape}

In Sec.~\ref{sec.stab} we have discussed at a semi-quantitative level
the role of bulk spin flips which produce finite-size effects
whose importance increases with decreasing $\sigma$. The behaviour of interfacial spins with
three aligned neighbours and the fourth pointing in the opposite direction (namely spins on a locally
flat interface) is more difficult to be investigated. However we can imagine that, at least at a qualitative level,
they should behave similarly to the bulk ones, in that they are stable for large $\sigma $ and bubble sizes $R$, and
become unstable upon decreasing $R$, particularly for small $\sigma$. According to this reasoning, these spins
together with the bulk ones introduce finite bubble size effects that are more severe for small values of $\sigma$.
As already mentioned, we believe these finite size effects to be responsible for the effective exponent
$z_{\rm eff}$ somewhat larger than $4/3$ observed for small $\sigma $ in Fig.~\ref{bubble} and in the lower panel of
Fig.~\ref{summary_exp}.
In order to avoid them, here we offer simulations where
flips in the bulk and with three aligned neighbours are forbidden.

Results of simulations conducted with this simplified dynamics at $T=0$ are shown in Fig.~\ref{bubble_simple}.
Our data clearly prove that the dynamical
exponent with the simplified dynamics is ${\cal Z}=4/3$ for all $\sigma$.
This suggests that this exponent is universal (i.e. independent of $\sigma$) and that $4/3$ is its truly asymptotic value.
However, according to the previous considerations, in order to see this value with a comparable evidence for small
$\sigma$ in full simulations one should access very large bubble sizes which is not feasible with our current numerical resources.

In order to clarify the physical mechanisms leading to ${\cal Z}=4/3$ and to the pre-asymptotic value ${\cal Z}=3/2$, it is necessary to have a more detailed study of the shape of the shrinking bubble.
If we focus on a quadrant and suppose that only edge spins can flip, the dynamics
of the droplet can be reduced to that of the terraces pictorially depicted in Fig.~\ref{pict_domain}.
We indicate with $\ell^{(j)}_r$ the length of the $j$-th terrace exactly at the time step when the top terrace
shrinks to zero, in the configuration when the bubble radius is (at that time) $r$.

Taking into account that the kink on the upper terrace always moves ballistically, we find that the time necessary for the droplet to disappear is
\begin{equation}
t(R) = \sum _{r=2}^{R} \ell_r^{(2)}+\ell_R,
\label{tl}
\end{equation}
where $\ell_R$ is the length of the top terrace in the initial configuration.
The evaluation of $\ell_r^{(2)}$ is complicated because of the ``long-range'' interaction between terraces. An analytical expression for $\ell_r^{(2)}$ can be obtained from the observation (see Fig.~\ref{bubble_geom}) that,
independently of the initial configuration, the droplet asymptotically assumes a quasi circular shape. Assuming that the droplet is {\it exactly} circular at all times,
the quantity $\ell_r^{(2)}$ is obtained from the equation
of the boundary in the continuum, $x^2+y^2=r^2$,
as the value of $x$ corresponding to $y=r-1$,
\begin{equation}
\ell_r^{(2)}=\sqrt{2r-1} \simeq \sqrt{2r}.
\label{topterr}
\end{equation}
Inserting this expression in Eq.~(\ref{tl}) we get $t(R)\sim R^{3/2}$, hence ${\cal Z}$=3/2.
This value, as we know, is observed for small $R$, whereas for large $R$ it crosses over
to ${\cal Z}=4/3$ indicating a deviation from the circular profile of the upper terraces.
This deviation must be attributed to the interaction between different terraces, in particular because of the drive exerted by the ballistic upper terrace transporting to some extent the lower ones. In order to be more quantitative about this,
in Fig.~\ref{terrace} we plot results for $\ell_r^{(j)}$ as a function of $r$ for $j\le 10$.
Let us look at the region $10^2<r<10^3$. In this sector one sees that the size of the
innermost of the considered terraces, namely $j=10$ (lowest curve, magenta) grows as $r^{1/2}$. Recalling Eq.~(\ref{topterr}), and observing that also $\ell ^{(j)}_r\propto \sqrt {r}$, for $r \gg j$,
this is what one would expect for a circular domain indicating that, sufficiently
far from the top of the domain, its shape is circular. However, upon decreasing $j$, that is considering
higher terraces, this slope gradually decreases until, for $j=2$, the top terrace behaves as
$\ell_r^{(2)}\sim r^{1/3}$. Using this result in Eq.~(\ref{tl})
one gets from the behaviour of this  terrace the correct asymptotic value ${\cal Z}=4/3$.
This shows that the origin of this number must be traced back to the non trivial geometry of the domains induced by
non local interactions between terraces. The actual shape of the domain during its shrinkage is shown in Fig.~\ref{profile},
where one can see deviations from circularity in the top and rightmost part of the bubble. Notice that such
deviations are hardly visible to the naked eye, however their effect on ${\cal Z}$ is important enough to change
its value from $3/2$ to $4/3$.

To complete the discussion of Fig.~\ref{terrace}, let us notice that there is a steep increase at $r\gtrsim 10^3$
and also that the algebraic behaviour of the curves is spoiled for small values of $r$. The former effect is due to the fact that the simulation is initialised with a domain of size $\gtrsim 1000$ and it takes some time for the dynamical process to build the domain's shape and to enter the scaling stage. The latter effect occurs because, in the very late stage of the process when the domain is very small, the scaling properties are lost because there are not enough interacting terraces to produce a many-body phenomenon.

\section{Conclusions}
\label{concl}

In this paper, we have shown that during the zero-temperature relaxation dynamics of the two-dimensional long-range Ising model a new dynamical regime appears.  This is characterised by an exponent ${\cal Z}$ that takes a pre-asymptotic value
${\cal Z}=3/2$ and an asymptotic one ${\cal Z}=4/3$. These values are strongly related to the interplay between the discrete nature of the lattice and the long-range nature of the interaction. In fact, if interactions are short-range (nn model)
we would have $z_{\rm cd}=2$ both in a discrete lattice and in a continuum picture (see Sec.~\ref{intro}).
On the other hand, a continuum picture with long-range interactions would trivially give $z=1$ (see Fig.~\ref{continuum} and related discussion at the end of Sec.~\ref{sec.stab}).

The pre-asymptotic value ${\cal Z}=3/2$ [see Eqs.~(\ref{tl}-\ref{topterr})] is simply understandable. For this, we use a discrete approximation of a circular bubble and assume that its shrinkage dynamics is limited by the closure of the largest terraces. However, Fig.~\ref{terrace} shows that the largest terraces scale as $r^{1/3}$ rather than as $r^{1/2}$, as expected for a circular bubble. This results in the true asymptotic value ${\cal Z}=4/3$. The reason for such a behaviour, see Sec.~\ref{Z_shape}, can be associated with the long-range character of the interaction and the necessity to locally preserve a convex shape of the bubble in order for the dynamics to proceed. Apart from this phenomenological understanding of the mechanism producing the asymptotic value ${\cal Z}=4/3$, a rigorous analytical derivation has not yet been possible.

Given the paucity of analytic approaches, further analysis is confined to better numerical simulations of the multi-domain system. To be useful, these must go beyond the time-scales and the statistical accuracy of those presented in this paper.  As we estimate below, this seems to be possible, but entails a huge computational effort.
Fig.~\ref{summary_exp} shows that the effective exponent would reasonably approach the value $4/3$ for sizes
of order $L(t) \simeq 200-300$. Thus, one should check for its further behaviour up to, say, $L(t) \simeq 500-600$.
Fig.~\ref{preasympt} informs us that, in order to reach this point, we should go to times of order $t\simeq 10^4$. Also, from the same figure one understands that finite-size effects start to be important around $L(t) \simeq 100 \simeq {\cal L}/20$. Hence, in order to reach $L(t) \simeq 500$ without severe finite-size effects, the system should have at least a size ${\cal L} \simeq 10^4$. Putting the above facts together, we estimate a required computational effort at least 100 times larger than the one of this study, which is already quite massive. This estimate applies for the larger values of $\sigma$, for smaller values the situation is much worse!

The study of this paper is restricted to a two-dimensional square lattice. However, since we have emphasised that
lattice effects are relevant in determining the $T=0$ growth exponent, an interesting question to be still addressed is the universality of the exponent ${\cal Z}$ with respect to the lattice geometry. In this direction, numerical simulations of $d=2$ systems on different lattices (e.g., triangular or hexagonal lattices) would be useful. Similarly, the dependence on dimensionality should also be considered, although the numerical effort required would be further increased. Finally, the dynamics with extremely long-range interactions (i.e., with $\sigma \le 0$) remains a completely unexplored subject in any dimension.

Few days after submitting this article for publication, a preprint~\cite{Christiansen_arxiv} has been uploaded
where the authors also find the exponent $3/4$ independently of $\sigma$.


\section{Operator renormalization integrals} \label{sec:ori}

The operator convolution integration is shown in Eq.~\eqref{eq:operator} with the splitting functions shown in Eq.~\eqref{eq:splitting}.
The splitting functions can be divided into three parts, while the LO FF have two regions. Thus there are $6$ types of integrals,
\begin{align}
	& \int_z^1 \frac{\mathrm{d} y}{y} \delta (1-y) (1-\frac{z}{y})^{-\epsilon} f_1\left(\frac{z}{y} \right)
		= \lambda^{-\epsilon} f_1(1-\lambda) \,,\\
	& \int_z^1 \frac{\mathrm{d} y}{y} \delta (1-y) (1-\frac{z}{y})^{-1-2\epsilon} f_2\left(\frac{z}{y} \right)
		= \lambda^{-1-2\epsilon} f_2(1-\lambda) \,,\\
	\begin{split}
	& \int_z^1 \frac{\mathrm{d} y}{y} g(y) (1-\frac{z}{y})^{-\epsilon} f_1\left(\frac{z}{y} \right) \\
	& = \lambda^{-\epsilon} \int_0^1 \mathrm{d} y
		h(y,\epsilon)
		\frac{\lambda}{1-\lambda y}
		g\left( \frac{1-\lambda}{1-\lambda y} \right) f_1(1-\lambda y) \,,
	\end{split}\\
	\begin{split}
	& \int_z^1 \frac{\mathrm{d} y}{y} g(y) (1-\frac{z}{y})^{-1-2\epsilon} f_2\left(\frac{z}{y} \right) \\
	&	= \lambda^{-1-2\epsilon}
		\Bigg[ -\frac{1}{2\epsilon} \lambda g(1-\lambda) f_2(1) \\
	& \phantom{=} + \int_0^1 \frac{\mathrm{d} y }{y}
		h(y,2\epsilon)
		\left(
		\frac{\lambda}{1-\lambda y} g\left( \frac{1-\lambda}{1-\lambda y} \right) f_2(1-\lambda y)
		- \lambda g(1-\lambda) f_2(1)
		 \right)
		 \Bigg] \,,
	\end{split}\\
	\begin{split}
	& \int_z^1 \frac{\mathrm{d} y}{y} \frac{1}{[1-y]_+} (1-\frac{z}{y})^{-\epsilon} f_1\left(\frac{z}{y} \right) \\
	& = \lambda^{-\epsilon} \ln \lambda \, f_1(1-\lambda)
	+ \lambda^{-\epsilon} \int_0^1 \frac{\mathrm{d} y}{1-y}
		\left( h(y,\epsilon) f_1(1-\lambda y)
		- \frac{1-\lambda}{1-\lambda y} f_1(1-\lambda ) \right)  \,,
	\end{split}\\
	\begin{split}
	& \int_z^1 \frac{\mathrm{d} y}{y} \frac{1}{[1-y]_+} (1-\frac{z}{y})^{-1-2\epsilon} f_2\left(\frac{z}{y} \right) \\
	& = \lambda^{-1-2\epsilon} \ln \lambda \, f_2(1-\lambda) \\
	& \phantom{=} + \lambda^{-1-2\epsilon}
		\Bigg[\left(-\frac{1}{2\epsilon} - H(-2\epsilon) -\ln (1-\lambda) \right) f_2(1) \\
	& \phantom{=} + \int_0^1 \frac{\mathrm{d} y}{1-y}
		\left(
		h(y,2\epsilon) \frac{ f_2(1-\lambda y) - f_2(1) }{y}
		- \frac{1-\lambda}{1-\lambda y} \left( f_2(1-\lambda ) - f_2(1) \right)
		\right) \Bigg]\,,
	\end{split}
\end{align}
where $\lambda =1-z$, $g(y) = (1-2y)/y + y(1-y) $ comes from the splitting function, $H(-2\epsilon)$ is the Harmonic number which can be expanded in the limit of $\epsilon \to 0$ and
\begin{equation}
	h(y,\epsilon) = y^{-\epsilon} = \sum_i \frac{(-\epsilon)^i}{i!} \ln^i (y)
\end{equation}
is the Taylor series at $\epsilon \to 0$.
Based on these results, we find that  two more behaviors $(1-z)^{-\epsilon}\ln (1-z) $ and $(1-z)^{-1-2\epsilon}\ln (1-z) $ are added to the region at $z\to 1$.

\section{Introduction}

Heavy quarkonium provides an ideal physical system to explore the strong interaction, as the heavy quarkonium production contains both perturbative and non-perturbative effects in QCD. Non-relativistic QCD (NRQCD) factorization \cite{Bodwin:1994jh} is the most widely used theory to explain quarkonium production so far.
When transverse momentum $p_T$ of the produced quarkonium is large, fragmentation mechanism dominates and factorization is easier to hold, because  long-distance interactions between quarkonium and initial-state particles are suppressed then.

Inclusive production differential cross section of a specific hadron $H$ at high $p_T$ can be calculated in collinear factorization \cite{Collins:1989gx},
\begin{equation}
	\mathrm{d} \sigma _{A+B\rightarrow H(p_T)+X} = \sum_i \mathrm{d} \hat{\sigma} _{A+B \rightarrow i(p_T /z)+X'} \otimes D_{i\to H}(z,\mu) + {\cal O} (1/p_T^2) \, ,
\end{equation}
where $i$ sums over all quarks and gluons, $z$ is the light-cone momentum fraction carried by $H$ with respect to the parent parton $i$, and $A$ and $B$ are colliding particles whose effect should be further factorized into partons if they are hadrons.
$d \hat{\sigma} _{A+B \rightarrow i(p_T /z)+X}$ are perturbatively calculable hard parts, while $D_{i\to H}(z,\mu)$ are non-perturbative but universal fragmentation functions (FFs) describing the probability of parton hadronizing to $H$ with momentum fraction $z$. For heavy quarkonium production, an important ${\cal O} (1/p_T^2)$ contributions are double parton FFs which can also be factorized \cite{Kang:2011mg,Kang:2014tta,Kang:2014pya,Fleming:2012wy,Fleming:2013qu}.
In FFs, there is a collinear factorization scale $\mu$ dependence, which cancels with the similar dependence in hard parts perturbatively order by order, leaving physical differential cross section independent of the scale. The evolution of single parton FFs with respect to $\mu$ are controlled by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation \cite{Gribov:1972ri,Altarelli:1977zs,Dokshitzer:1977sg}, and similar evolution equations for double parton FFs are calculated in \cite{Kang:2014tta}. With these evolution equations, the only unknown information for FFs are their values at any chosen factorization scale $\mu=\mu_f$.

When $\mu_f$ is close to the quarkonium mass $m_H$, it is natural to calculate FFs via NRQCD factorization. For single parton FFs that will be considered in this paper, we have
\begin{equation} \label{eq:ff}
	D_{i\to H}(z,\mu_f) = \sum_n d_{i\to Q\bar{Q}(n)} (z,\mu_f) \langle \bar{\mo}^{H}_n\rangle \, ,
\end{equation}
where $d_{i\to Q\bar{Q}(n)}$ represent the perturbative calculable short-distance coefficients (SDCs) to produce a heavy quark-antiquark pair $Q\bar Q$ with quantum number $n$, and $\langle \bar{\mo}^{H}_n\rangle$ are normalized long-distance matrix elements (LDMEs) \footnote{ $\langle \bar{\mo}^{H}_n\rangle$ can be related to the original definition of NRQCD LDME $\mohn$ \cite{Bodwin:1994jh} by the following rules. They are the same if $n$ is color-octet, and $\langle \bar{\mo}^{H}_n\rangle=\mohn/(2N_c)$ if $n$ is color-singlet.}. The quantum number is usually expressed in spectroscopic notation $n=\state{{2S+1}}{L}{J}{c}$, with $c=1, 8$ respectively for color-singlet state and color-octet state. According to velocity scaling rule \cite{Bodwin:1994jh},  $\langle \bar{\mo}^{H}_n\rangle$ is usually suppressed if $L$ is too large. Therefore, the most important states for phenomenological purpose are $S$-wave and $P$-wave states. Because LDMEs are supposed to be process independent, they can be determined by fitting experimental data, while SDCs can be calculated perturbatively through the matching procedure.

All the SDCs for single parton FFs to both $S$-wave and $P$-wave states up to $\mathcal{O}(\alpha_s^2)$ are available \cite{ Braaten:1993mp,Braaten:1993rw,Cho:1994gb,Braaten:1994kd,Beneke:1995yb,Ma:1995vi,Braaten:1996rp,Braaten:2000pc,Hao:2009fa,Jia:2012qx,Bodwin:2014bia} (see  \cite{Ma:2013yla,Ma:2015yka} for a summary and comparison). At $\alpha_s^3$ order, the leading order (LO) SDCs of $g \rightarrow Q\bar{Q}(\CScSa)+X$ and $g \rightarrow Q\bar{Q}(\CSaPa)+X$ are calculated separately in Refs. \cite{Zhang:2017xoj,Braaten:1993rw,Braaten:1995cj,Bodwin:2003wh,Bodwin:2012xc} and Ref. \cite{Sun:2018yam}, and the next-to-leading order (NLO) SDCs of $g\to Q\bar{Q}(\state{1}{S}{0}{1,8})+X$ have also been calculated in Refs. \cite{Zhang:2018mlo,Artoisenet:2018dbs,Feng:2018ulg} recently. To perform phenomenological analysis at a consistent precision, such as that for ${J/\psi}$ or $\chi_{cJ}$ production, at this order, one also needs the NLO SDCs of $g \rightarrow Q\bar{Q}(\state{3}{P}{J}{1,8})+X$ and the next-to-next-to-leading order (NNLO) SDCs of $g \rightarrow Q\bar{Q}(\COcSa)+X$. In this paper, we will focus on the former.

Another motivation to calculate the NLO FFs for processes $g \rightarrow Q\bar{Q}(\state{3}{P}{J}{1,8})+X$ is to test the NRQCD factorization conjecture at two-loop level with a physical process. The test at two-loop  level exists in literature \cite{Nayak:2005rw,Nayak:2005rt,Nayak:2006fm,Bodwin:2019bpf} but within the framework of  eikonal approximation  for infrared (IR) divergences. The eikonal approximation introduces artificial ultraviolet (UV) divergences which makes the extraction of IR divergences  beyond one-loop level highly nontrivial. Because we do not use eikonal approximation, the extraction of IR divergences in our case is straightforward.

In this paper, we aim to calculate NLO SDCs of $g \rightarrow Q\bar{Q}(\state{3}{P}{J}{1,8})+X$ to high precision using similar methods in our previous paper \cite{Zhang:2018mlo}. The rest of the paper is organized as following. In Sec. \ref{sec:full}, we calculate the process to LO and NLO in full-QCD. The calculation method is very similar to that in Ref. \cite{Zhang:2018mlo}. A small improvement is the computation of operator renormalization shown in Sec. \ref{sec:nlo} and Appendix. \ref{sec:ori}.
We argue that the NRQCD factorization at two-loop level without eikonal approximation is verified in Sec. \ref{sec:test}. The $\mathcal{O}(\alpha_s^2)$ perturbative NRQCD matrix elements (PMEs) $\LDME{\CScPj}{\COcSa}_{\mathrm{NLO}}$ and $\LDME{\COcPj}{\COcSa}_{\mathrm{NLO}}$ are matched, and the corresponding IR-safe SDCs are obtained simultaneously in Sec. \ref{sec:match}. Discussion will be presented in Sec. \ref{sec:sdcdis}.

\section{Full-QCD Calculation} \label{sec:full}

\subsection{Definition}

FF of a gluon to a hadron (quarkonium) is defined by Collins and Soper \cite{Collins:1981uw},
\begin{align}\label{eq:defFF}
\begin{split}
    D_{g \rightarrow H}(z,\mu_0)=
    & \frac{-g_{\mu\nu}z^{D-3}}{2 \pi P_c^{+}(N_{c}^{2}-1)(D-2)} \int_{-\infty}^{+\infty}\mathrm{d}x^{-} e^{-i P_c^{+} x^{-}} \\
    & \times \langle 0 | G_{c}^{+\mu}(0) \mathcal{E}^{\dag}(0,0,\boldsymbol{0}_{\perp})_{cb} \mathcal{P}_{H(P)} \mathcal{E}(0,x^{-},\boldsymbol{0}_{\perp})_{ba} G_{a}^{+\nu}(0,x^{-},\boldsymbol{0}_{\perp}) | 0 \rangle \, ,
\end{split}
\end{align}
where $G^{\mu\nu}$ is the gluon field-strength operator, $P$ and $P_c$ are respectively the momenta of the produced hadron $H$ and the initial-state fragmenting gluon $g$, and $z=P^+/P_c^+$ is the ratio of momenta along the ``+'' direction.
It is convenient to choose the frame in which the hadron has zero transverse momentum, $P = (z P_c^+, M^2/(2 z P_c^+),\boldsymbol{0}_{\perp})$, with $P^2=2P^+P^-=M^2$.
The projection operator $\mathcal{P}_{H(P)}$ is defined by
\begin{equation}\label{eq:projectH}
\mathcal{P}_{H(P)} = \sum_X |H(P)+X \rangle \langle H(P)+X|\,,
\end{equation}
where $X$ sums over all unobserved particles.
The gauge link $\mathcal{E}(x^{-})$ is an eikonal operator that involves a path-ordered exponential of gluon field operators along a light-like path,
\begin{equation}
  \mathcal{E}(0,x^{-},\boldsymbol{0}_{\perp})_{ba}= \mathrm{P} \, \text{exp} \left[+i g_s \int_{x^{-}}^{\infty}\mathrm{d}z^{-} A^{+}(0,z^{-},\boldsymbol{0}_{\perp}) \right]_{ba} \, ,
\end{equation}
where $g_s=\sqrt{4 \pi \alpha_s}$ is the QCD coupling constant and $A^{\mu}(x)$ is the matrix-valued gluon field in the adjoint representation: $[A^{\mu}(x)]_{ac} = i f^{abc} A^{\mu}_{b}(x)$.

Since the SDCs $d$'s in the NRQCD factorization formula (\ref{eq:ff}) are independent of the final state $H$, we can replace $H$ by on-shell $Q\bar{Q}$ states to extract $d_{i\to Q\bar{Q}[n]}$ through the matching procedure [See below (\ref{eq:matchlo}) and (\ref{eq:matchnlo})].

Choosing the factorization scale $\mu_0 \sim 2m_Q$ with $m_Q$ being the heavy quark mass, we can calculate FF for $g\to Q\bar{Q}[n]$ in the full-QCD with the Feynman gauge.
Most of the Feynman rules are the same as those of QCD, while the others related to the eikonal line can be found in Ref.~\cite{Zhang:2018mlo}. Feynman amplitudes are denoted as $\mathcal{M}_{\lambda_{{Q}}\lambda_{\bar{Q}}\lambda_i}(q)$, where $\lambda_{{Q}}$ and $\lambda_{\bar{Q}}$ are respectively spins of produced on-shell heavy quark and heavy antiquark, $\lambda_i$ ($i=1, 2, \dots$) are spins of the initial-state virtual gluon or final-state unobserved light particles, and $q$ is the heavy quark momentum in the rest frame of heavy quark pair.
To project the free $Q\bar{Q}$ pair to the state with spin-triplet and specific color, one can multiply amplitudes by projection operators and obtain
\begin{equation}
	\mathcal{M}^\alpha_{\lambda_i}(q)=\text{Tr}\left[ \Gamma_c \Gamma_\lambda^\alpha \overline{\mathcal{M}}_{\lambda_{{Q}}\lambda_{\bar{Q}}\lambda_i}(q)  \right]\, ,
\end{equation}
where $\overline{\mathcal{M}}_{\lambda_{{Q}}\lambda_{\bar{Q}}\lambda_i}$ denotes amplitude with external heavy-(anti)quark spinors removed and the projection operators are defined as
\begin{align}
\begin{split}
	&\Gamma_{c=1}
	 = \frac{1} {\sqrt{N_{c}}} \,, \\
	&\Gamma_{c=8}
	 = \frac{\sqrt{2} T^{a}} {\sqrt{N_{c}^{2}-1}} \,, \\
	&\Gamma_\lambda^\alpha
	 = \frac{1} {\sqrt{2 E} (E + m_Q)}
		(\slashed{\overline{p}} - m_Q)
		\frac{2 E - \slashed{P}}{4 E}
		\gamma^\alpha
		\frac{2 E + \slashed{P}}{4 E}
		(\slashed{p} - m_Q)	\, ,
\end{split}
\end{align}
where $p=P/2+q$ and $\overline{p}=P/2-q$ are momenta respectively for $Q$ and $\bar{Q}$,  $m_Q$ is the mass of heavy quark, and $E$ is the energy of $Q$ or $\bar{Q}$ in the rest frame of $Q\bar{Q}$. With the on-shell conditions of $Q$ and $\bar{Q}$, we have
\begin{align}
	P\cdot q =0,  \quad \quad P^2=4E^2=4(m_Q^2-q^2).
\end{align}
To project the amplitude into P-wave, we also need to take the derivation over $q^\beta$,
\begin{equation}
	\mathcal{M}^{\alpha \beta}_{\lambda_i} =
	 \frac{\partial}{\partial q^\beta}
	\mathcal{M}^\alpha_{\lambda_i}(q) |_{q\to 0} \,.
\end{equation}
And by summing over spin and color of initial-state and final-state particles, we get the squared amplitude
\begin{equation}
	\mathcal{A}_J  = \frac{N_{\mathrm{CS}}}{N_J} \overline{\sum}
	P_J^{\alpha \alpha' \beta \beta'} \,
	\mathrm{Re} \left[ \mathcal{M}^{\alpha \beta}_{\lambda_i} \mathcal{M}^{*\alpha' \beta'}_{\lambda_i} \right]\,,
\end{equation}
where $N_{\mathrm{CS}}=\frac{z^{D-2}}{(N_{c}^{2}-1)(D-2)} $ with $D=4-2\epsilon$ is the space-time dimension.
For different $J$, $N_J$ and $P_J^{\alpha \alpha' \beta \beta'}$ are respectively
\begin{equation}
	N_J = \left\{
	\begin{aligned}
 		& 1 \,, & \text{for } J=0 \\
 		& \frac{1}{2}(D-1)(D-2) \,, & \text{for }    J=1 \\
 		& \frac{1}{2}(D+1)(D-2)\,, & \text{for } J=2 \\
 		& (D-1)^2 \,, & \text{for } \sum_J
 	\end{aligned}
	\right. \,,
\end{equation}
and
\begin{equation} \label{eq:project}
	P_J^{\alpha \alpha' \beta \beta'} = \left\{
	\begin{aligned}
 		& \frac{1}{D-1} I^{\alpha \beta} I^{\alpha' \beta'} \,, & \text{for } J=0 \\
 		& \frac{1}{2}\left( I^{\alpha \alpha'} I^{\beta \beta'}-I^{\alpha \beta'} I^{\alpha' \beta} \right) \,, & \text{for }    J=1 \\
 		& \frac{1}{2}\left( I^{\alpha \alpha'} I^{\beta \beta'}-I^{\alpha \beta'} I^{\alpha' \beta} \right) - \frac{1}{D-1} I^{\alpha \beta} I^{\alpha' \beta'} \,, & \text{for } J=2 \\
 		& I^{\alpha \alpha'} I^{\beta \beta'} \,, & \text{for } \sum_J
 	\end{aligned}
	\right. \,,
\end{equation}
where
\begin{equation}
	I^{\alpha\beta} = -g^{\alpha \beta}
			+ \frac{P^{\alpha} P^{\beta}} {P^{2}} \, .
\end{equation}
Then, fragmentation function in the full-QCD can be denoted as
\begin{equation} \label{eq:FFint}
	D[g\to Q\bar{Q} (\state{3}{P}{J}{1/8}])] = \int \mathrm{d} \Phi \mathcal{A}_J \, .
\end{equation}
Here, $\mathrm{d} \Phi$ is the final-state phase space,
\begin{equation}
    \mathrm{d} \Phi =
  	\frac {P^{+}}{z^{2} S}
  	\delta \left( \frac{1-z}{z} P^{+} - \sum_i k_{i}^{+} \right)
    \prod_{i} \frac{\mathrm{d}^{D} k_i}{(2\pi)^{D-1}}
        \delta _+ (k_i^{2})
         \,,
\end{equation}
where $S$ is the symmetry factor for final-state particles and $k_i$'s are momenta of final-state light particles.
In the rest of the paper, we will take the case of $g\to Q\bar{Q} (\CScPz)+X$ as an explicit example to explain our calculations.

\subsection{LO FFs} \label{sec:lo}
The Feynman diagrams of gluon fragmenting into $Q\bar{Q}[\state{3}{P}{0}{1}]$ at the LO in $\alpha_s$ are shown in Figure  \ref{fig:lofd}.
\begin{figure}[htb!]
 \begin{center}
 \includegraphics[width=120pt]{loFD-eps-converted-to.pdf}
  \caption{One of the two Feynman diagrams of gluon fragmenting into $Q\bar{Q}[\state{3}{P}{0}{1}]$ at LO in $\alpha_s$. Another diagram can be obtained by permuting the heavy quark and anti-quark.  \label{fig:lofd}}
 \end{center}
\end{figure}
These diagrams can be easily calculated, and the obtained FF can be expressed as
\begin{equation}\label{eq:lo}
	D_{\mathrm{LO}}[g\to Q\bar{Q} (\CScPz)] = \frac{N_{\mathrm{LO}}}{1-\epsilon}
 		\left[ (1-z)^{-\epsilon} f_1(z) + (1-z)^{-1-2\epsilon} f_2(z) \right] \, ,
\end{equation}
where
\begin{equation}\label{eq:factorlo}
	N_{\mathrm{LO}} = \frac{ \, \alpha_s^2 }{(D-1) N_c m_Q^5 }
 		\left( \frac{\pi \mu_r^2}{m_Q^2} \right)^{\epsilon} \Gamma(1+\epsilon)
\end{equation}
and
\begin{align}
\begin{split}
	f_1(z) & = \frac{3(5-3z)}{\epsilon }+\frac{3}{2} (3 z-4)-\frac{9 z  }{2}\epsilon \,,\\
	f_2(z) & = \frac{3(1-z) (3z-5)}{\epsilon }+\frac{13 z^3-42 z^2+15 z+18}{3}  +2 \left(2 z^2-2 z-1\right) z \epsilon +\frac{2 z^3 \epsilon ^2}{3}

  \,.
\end{split}
\end{align}
With the following decomposition
\begin{equation}
	(1-z)^{-1-2\epsilon} = -\frac{\delta(1-z)}{2 \epsilon }+\frac{1}{ [1-z]_+}-2 \epsilon  \left[ \frac{\ln(1-z)  }{1-z}\right]_+ + \mathcal{O} (\epsilon^2) \,,
\end{equation}
the FF can be expressed in the limit of $\epsilon \to 0$ as
\begin{align}\label{eq:loexp}
\begin{split}
	D_{\mathrm{LO}}[g\to Q\bar{Q} (\CScPz)] = & N_{\mathrm{LO}}
	\Big( -\frac{2}{3 \epsilon} \delta(1-z)+\frac{1}{3} \delta(1-z)+\frac{4}{3} \frac{z}{[1-z]_+} \\
	& \phantom{ N_{\mathrm{LO}}\Big( }
	+\frac{1}{6} (85-26 z) z+3 (5-3 z) \ln(1-z) \Big) \,.
\end{split}
\end{align}
It is clear that there is an IR divergence in the full-QCD result.
This divergence can be absorbed by the color-octet PME, which will be explained in Sec. \ref{sec:match}.

\subsection{NLO FFs} \label{sec:nlo}

The NLO full-QCD calculations of this process resemble closely to the NLO calculations of $g\to Q\bar{Q}(\state{1}{S}{0}{1,8})+X$ in Ref. \cite{Zhang:2018mlo}.
The relevant Feynman diagrams are shown in Figure \ref{fig:nlofd} and Figure \ref{fig:nloekfd}. In our calculation, we need to cut these diagrams at some specific positions. In Figure \ref{fig:cutfd}, we give an example to show how to obtain cut diagrams for the real and  the virtual corrections.
\begin{figure}[htb!]
 \begin{center}
 \includegraphics[width=\textwidth]{nloFD-eps-converted-to.pdf}
  \caption{Typical Feynman diagrams without gluon eikonal vertex for $g\to Q\bar{Q}(\state{3}{P}{0}{1})+X$. The other diagrams can be obtained by permuting the heavy quark and anti-quark. The real and virtual diagrams can be obtained by cutting at different positions.  \label{fig:nlofd}}
 \end{center}
\end{figure}

\begin{figure}[htb!]
 \begin{center}
 \includegraphics[width=0.8\textwidth]{nloekFD-eps-converted-to.pdf}
  \caption{Typical Feynman diagrams with gluon eikonal vertex for $g\to Q\bar{Q}(\state{3}{P}{0}{1})+X$. The other diagrams can be obtained by permuting the heavy quark and anti-quark. The real and virtual diagrams can be obtained by cutting at different positions.  \label{fig:nloekfd}}
 \end{center}
\end{figure}

\begin{figure}[htb!]
 \begin{center}
 \includegraphics[width=0.5\textwidth]{cutFD-eps-converted-to.pdf}
  \caption{ Real and virtual cut diagrams by cutting at different positions for the last diagrams in Figure \ref{fig:nloekfd}.  \label{fig:cutfd}}
 \end{center}
\end{figure}

The calculation is almost the same as that in Ref. \cite{Zhang:2018mlo}. Firstly we generate the amplitude from the Feynman diagrams for real or virtual corrections. Then reduce them into linear combinations of master integrals (MIs). There are 95 MIs for real corrections and 66 MIs for virtual corrections.
Fortunately, these MIs are the same as those in Ref. \cite{Zhang:2018mlo}, which can be calculated by using differential equations (DEs) method \cite{Kotikov:1990kg,Bern:1992em,Remiddi:1997ny,Gehrmann:1999as,Henn:2013pwa, Lee:2014ioa,Adams:2017tga,Caffo:2008aw,Czakon:2008zk,Mueller:2015lrx,Lee:2017qql,Liu:2017jxz,Liu:2020kpc}.
We can estimate values of these MIs in regions $0\sim1/4$, $1/4\sim3/4$ and $3/4\sim1$ respectively by the asymptotic expansions of them at $z=0,1/2$ and $1$.
Finally, the real and virtual corrections can be expressed by asymptotic expansions in different regions as
\begin{equation}\label{eq:rvexp}
	\sum_{s} \sum_{i=0}^{n_{s}} (z-z_0)^{s} \ln^{i}(z-z_0) \sum_{j=0}^{\infty} I^{ s \, i \, j}(\epsilon) (z-z_0)^j \, ,
\end{equation}
where $s$ are linear functions of $\epsilon$, $n_s$ are integers determined by $s$, and coefficients $I^{ s \, i \,j }(\epsilon)$ are functions of $\epsilon$, which can be computed to sufficient high order with high precision.

Renormalization has some differences from that in Ref. \cite{Zhang:2018mlo}, because $1/\epsilon$ pole already presents at LO FF.
Firstly, the on-shell renormalization constants of heavy quark field and mass ($\delta_2$ and $\delta_m$) should be expanded at least to $\mathcal{O}(\epsilon)$.
To simplify the procedure, we use expression with exact $\epsilon$ dependence,
\begin{equation}
	\delta_i = \frac{\alpha_s}{\pi} \Gamma(1+\epsilon)
				\left( \frac{\pi \mu_r^2}{m_Q^2} \right)^{\epsilon}
				\hat{\delta} _i \, ,
\end{equation}
with
\begin{align}\label{eq:zren}
\begin{split}
	\hat{\delta}_2 & = -\frac{C_F}{4^{1-\epsilon}}
						\left( \frac{1}{\epsilon_{\mathrm{UV}}} \frac{1+2 \epsilon}{1-2 \epsilon } +\frac{2}{\epsilon_{\mathrm{IR}}} \right) \,, \\
	\hat{\delta}_m & = -\frac{C_F}{4^{1-\epsilon}}
						\frac{1}{\epsilon_{\mathrm{UV}}} \frac{3-2 \epsilon }{1-2 \epsilon } \,,
\end{split}
\end{align}
where $C_F=(N_c^2-1)/(2N_c) $.

Apart from the calculations of the QCD counter terms,  some tricks are needed in the calculation of operator renormalization.
In the $\overline{\mathrm{MS}}$ scheme, the operator counter term is given by
\begin{equation}\label{eq:operator}
	D_{\mathrm{NLO} }^{\mathrm{Operator}} (z) =
		-\frac{\alpha_s}{2\pi}
		\frac{\Gamma(1+\epsilon)}{\epsilon}
		\left( \frac{4 \pi \mu_r^2}{\mu_0^2} \right)^{\epsilon}
		\int_z^1 \frac{\mathrm{d} y}{y} P_{gg} (y)
			D_{\mathrm{LO}} \left(\frac{z}{y} \right) \, ,
\end{equation}
where $\mu_0$ is the collinear factorization scale, and the Altarelli-Parisi splitting function $P_{gg} (z)$ is
\begin{equation}\label{eq:splitting}
	P_{gg} (z) = b_0 \, \delta (1-z) +
				 2 N_c \left( \frac{z}{(1-z)_+} + \frac{1-z}{z} + z(1-z) \right) \, ,
\end{equation}
where $b_0 = (11N_c-2n_f)/6$.
If one takes the Eq.~\eqref{eq:loexp} as LO FF $D_{\mathrm{LO}}(z)$, there will be two plus functions convolved in the integral in (\ref{eq:operator}), which is difficult to define and calculate.
However, we can take Eq.~\eqref{eq:lo} as the input form of the $D_{\mathrm{LO}}(z)$ to overcome this difficulty, with corresponding integrated results given in Appendix \ref{sec:ori}.
In this way, operator renormalization results can be easily obtained in each region, which are convenient to be added to the real and virtual corrections.

\subsection{Verification of the NRQCD factorization at two-loop level}
\label{sec:test}

In the full-QCD calculations, the FF can be expressed by the phase space integral of amplitude squared as shown in Eq.~\eqref{eq:FFint}. Before taking the derivative over $q$, the amplitude squared at ${\cal{O}}(q^2)$ can be expressed as
\begin{equation}\label{eq:ampqn}
	\mathcal{A} (q,q') =  (q\cdot q')\mathcal{A}_1 + \frac{(q\cdot n) (q'\cdot n)}{(P\cdot n) ^2} \mathcal{A}_2 \,,
\end{equation}
where $n$ is the light-like vector that defines the direction of the eikonal line of fragmentation function. Because NRQCD LDMEs should be independent of the direction of the eikonal line of fragmentation function, and all remained IR divergences in full-QCD calculation must be absorbed into the LDMEs if NRQCD factorization is valid, there should be no IR divergences in $\mathcal{A}_2$. On the other hand, if one can show that there is no IR divergence in $\mathcal{A}_2$, then all IR divergences in full-QCD FF can be absorbed into the NRQCD LDMEs in principle, and the NRQCD factorization can be verified at two-loop level.

Two-loop verification of NRQCD factorization was also done in Refs. \cite{Nayak:2005rw,Nayak:2005rt,Nayak:2006fm,Bodwin:2019bpf}, in which eikonal approximation was used to simplify the calculation. The eikonal approximation however introduced artificial UV divergences that makes the extraction of IR divergences highly nontrivial.

In this subsection, we verify the NRQCD factorization without applying eikonal approximation. We take the process $g\to Q\bar{Q} (\CScPj)+X$ with $J$ summed as an example.
To extract the part $\mathcal{A}_2$ only, we use the following projection operator
\begin{equation}
	P_{qn}^{\alpha \alpha' \beta \beta'} =  I^{\alpha \alpha'}
	\left( (D-2)I^{\beta \beta'} -(D-1)d^{\beta \beta'} \right)  \,,
\end{equation}
where
\begin{equation}
	d^{\beta \beta'} = -g^{\beta \beta'} + \frac{P^{\beta} n^{ \beta'} + P^{ \beta'} n^{\beta}} {P\cdot n}
			- \frac{P^2 n^{\beta} n^{ \beta'}} {(P\cdot n)^{2}} \,,
\end{equation}
and we get
\begin{equation}
	P_{qn}^{\alpha \alpha' \beta \beta'} \, \frac{\partial}{\partial q^\beta \, \partial q'^{\beta'}} \mathcal{A} (q,q')
	= I^{\alpha \alpha'} (D-2)\frac{1}{P^2} \mathcal{A}_2  \,.
\end{equation}
Performing the calculations similar to those introduced in Sec.~\ref{sec:full}, we find that all divergences are canceled out after adding real corrections, virtual corrections as well as UV renormalization contributions. As a result, IR divergences in $\mathcal{A} (q,q')$ are free of the gauge line vector $n$. Thus the NRQCD factorization is verified to hold at two-loop level, which is consistent with the observation in Refs. \cite{Nayak:2005rw,Nayak:2005rt,Nayak:2006fm,Bodwin:2019bpf}.

\section{Matching procedure and the SDCs}

\label{sec:match}

As  SDCs are independent of specifics of the hadronic states,  we can match the $\mathcal{O}(\alpha_s^2)$ PMEs from the divergences in full-QCD FFs and obtain the IR-safe SDCs simultaneously, thanks to the fact that NRQCD factorization holds at this order.

For $g\to Q\bar{Q}(\state{3}{P}{J}{1/8} )$ process, the LO FF can be expressed by
\begin{align}\label{eq:matchlo}
\begin{split}
	D_{\mathrm{LO}}[g\to Q\bar{Q} (\state{3}{P}{J}{1/8})] =
	& \phantom{+} d_{\mathrm{LO}}[g\to Q\bar{Q} (\state{3}{P}{J}{1/8})] \LDME{\state{3}{P}{J}{1/8}}{\state{3}{P}{J}{1/8}}_{\mathrm{LO}} \\
	& + d_{\mathrm{LO}}[g\to Q\bar{Q} (\COcSa)] \LDME{\state{3}{P}{J}{1/8}}{\COcSa}_{\mathrm{LO}} \,,
\end{split}
\end{align}
where $\LDME{\CScPj}{\CScPj}_{\mathrm{LO}}=\LDME{\COcPj}{\COcPj}_{\mathrm{LO}} =1$ is normalized. The LO SDC of gluon to $\COcSa$ in (\ref{eq:matchlo}) has been calculated in Ref. \cite{Braaten:1996rp,Braaten:2000pc}, which is given in $D-$dimension by
\begin{equation}
	d_{\mathrm{LO}}[g\to Q\bar{Q} (\COcSa)] =
	\frac{\pi  \alpha_s \delta (1-z)}{ (D-1) (N_c^2-1) m_Q^3} \,,
\end{equation}
and the $\mathcal{O}(\alpha_s)$ PMEs in the $\overline{\mathrm{MS}}$ subtraction scheme can be expressed as
\begin{align}
	\LDME{\CScPj}{\COcSa}_{\mathrm{LO}} &=
	- C_F \frac{4\alpha_s}{3 \pi m_Q^2}
	\frac{\Gamma(1+\epsilon)}{\epsilon_\mathrm{IR}}
	\left( \frac{4 \pi \mu_r^2}{\mu_\Lambda^2} \right)^{\epsilon} \,,\\
	\LDME{\COcPj}{\COcSa}_{\mathrm{LO}} &=
	- B_F \frac{4\alpha_s}{3 \pi m_Q^2}
	\frac{\Gamma(1+\epsilon)}{\epsilon_\mathrm{IR}}
	\left( \frac{4 \pi \mu_r^2}{\mu_\Lambda^2} \right)^{\epsilon} \,,
\end{align}
where $B_F=(N_c^2-4)/(4N_c) $ and $\mu_\Lambda$ is the NRQCD factorization scale.
With full-QCD FF in Eq.~\eqref{eq:loexp}, the $D-$dimension SDC $d_{\mathrm{LO}}[g\to Q\bar{Q} (\CScPj)]$ can be matched as
\begin{align}\label{eq:sdclomu}
\begin{split}
	d_{\mathrm{LO}}[g\to Q\bar{Q} (\CScPz)] =&\frac{ \, \alpha_s^2 \Gamma(1+\epsilon)}{(D-1) N_c m_Q^5 }
 		\left( \frac{\pi \mu_r^2}{m_Q^2} \right)^{\epsilon} \times \Bigg(  \left(\frac{1}{3} - \frac{2}{3} \ln \left(\frac{\mu_\Lambda ^2}{4 m_Q^2} \right) \right) \delta(1-z)  \\
	&+\frac{4}{3} \frac{z}{[1-z]_+}
	 +\frac{1}{6} (85-26 z) z+3 (5-3 z) \ln(1-z) \Bigg) \,,
\end{split}
\end{align}
which equals to the finite part of the full-QCD FF in Eq.~\eqref{eq:loexp} if factorization scale is chosen as $\mu_\Lambda = 2m_Q$.
We find that the IR divergence in full-QCD FF $D_{\mathrm{LO}}[g\to Q\bar{Q} (\CScPz)]$ is absorbed into the $\mathcal{O}(\alpha_s)$ PME $\LDME{\CScPj}{\COcSa}_{\mathrm{LO}}$, and the SDCs are IR-safe up to LO.
This result is consistent with that in Ref. \cite{Ma:2013yla}.

Similarly, $D_{\mathrm{NLO}}[g\to Q\bar{Q} (\state{3}{P}{J}{1/8})]$ can be expressed in the NRQCD factorization as
\begin{align}\label{eq:matchnlo}
\begin{split}
	D_{\mathrm{NLO}}[g\to Q\bar{Q} (\state{3}{P}{J}{1/8})] = \phantom{+}
	& d_{\mathrm{NLO}}[g\to Q\bar{Q} (\state{3}{P}{J}{1/8})] \LDME{\state{3}{P}{J}{1/8}}{\state{3}{P}{J}{1/8}}_{\mathrm{LO}} \\
	+ & d_{\mathrm{LO}}[g\to Q\bar{Q} (\state{3}{P}{J}{1/8}))] \LDME{\state{3}{P}{J}{1/8}}{\state{3}{P}{J}{1/8}}_{\mathrm{NLO}} \\
	+ & d_{\mathrm{NLO}}[g\to Q\bar{Q} (\COcSa)] \LDME{\state{3}{P}{J}{1/8}}{\COcSa}_{\mathrm{LO}} \\
	+ & d_{\mathrm{LO}}[g\to Q\bar{Q} (\COcSa)] \LDME{\state{3}{P}{J}{1/8}}{\COcSa}_{\mathrm{NLO}} \,,
\end{split}
\end{align}
where all LO terms are known. For NLO terms, $D_{\mathrm{NLO}}[g\to Q\bar{Q} (\state{3}{P}{J}{1/8})]$ are calculated in this work, $\LDME{\state{3}{P}{J}{1/8}}{\state{3}{P}{J}{1/8}}_{\mathrm{NLO}}$ are well-known, which have only Coulomb divergences  and vanish in dimensional regularization with $\overline{\mathrm{MS}}$ subtraction scheme, and $d_{\mathrm{NLO}}[g\to Q\bar{Q} (\COcSa)]$ has been calculated in Refs. \cite{Braaten:2000pc,Ma:2013yla}, whose $D-$dimension form is given by
\begin{align}
\begin{split}
	d_{\mathrm{NLO}}[g\to Q\bar{Q} (\COcSa)] = &
	\frac{ \alpha_s^2 \Gamma(1+\epsilon)}{ 4(D-1) C_F m_Q^3}
	\left( \frac{\pi \mu_r^2}{m_Q^2} \right)^{\epsilon}
	\Bigg[ A(\mu_r) \delta(1-z) + \frac{1}{N_c} P_{gg} (z) \times  \\
	& \left( \ln \left(\frac{\mu_0 ^2}{4 m_Q^2} \right) -1 \right)+\frac{2(1-z)}{z} - \frac{4(1-z+z^2)^2}{2} \left[\frac{\ln(1-z)}{1-z} \right]_+
	 \Bigg] \,,
\end{split}
\end{align}
with
\begin{equation}
	A(\mu_r) = \frac{b_0}{N_c} \left( \ln \left(\frac{\mu_r ^2}{4 m_Q^2} \right) +\frac{13}{3} \right)
	+ \frac{4}{N_c^2} -\frac{\pi^2}{3}+ \frac{16}{3}\ln 2 \,.
\end{equation}
Therefore, SDCs $d_{\mathrm{NLO}}[g\to Q\bar{Q} (\state{3}{P}{J}{1/8})]$ and PMEs $\LDME{\state{3}{P}{J}{1/8}}{\COcSa}_{\mathrm{NLO}}$ can be extracted from Eq.~\eqref{eq:matchnlo} thanks to the fact that the former one is finite and the later one is IR-divergent defined in $\overline{\mathrm{MS}}$ subtraction scheme.

The $\mathcal{O}(\alpha_s^2)$ PMEs $\LDME{\state{3}{P}{J}{1/8}}{\COcSa}_{\mathrm{NLO}}$ are obtained numerically with high precision.
Using PSLQ algorithm \cite{ferguson1999analysis,Bailey:1999nv}, analytical expressions are obtained as
\begin{align}
	\LDME{\CScPj}{\COcSa}_{\mathrm{NLO}} &=
	- C_F \frac{\alpha_s^2}{27 \pi^2 m_Q^2}
	\left( \frac{4 \pi \mu_r^2}{\mu_\Lambda^2} \right)^{2\epsilon}
	\Gamma(1+\epsilon)^2
	\left( \frac{9b_0}{\epsilon^2}
	+\frac{5n_f+(12\pi^2-47)N_c}{\epsilon} \right)   \,,\\
		\LDME{\COcPj}{\COcSa}_{\mathrm{NLO}} &=
	- B_F \frac{\alpha_s^2}{27 \pi^2 m_Q^2}
	\left( \frac{4 \pi \mu_r^2}{\mu_\Lambda^2} \right)^{2\epsilon}
	\Gamma(1+\epsilon)^2
	\left( \frac{9b_0}{\epsilon^2}
	+\frac{5n_f+(3\pi^2-47)N_c}{\epsilon} \right)   \,.
\end{align}

Setting factorization scales $\mu_0=\mu_\Lambda=2m_Q$, the NLO SDC can be expressed as
\begin{align}\label{eq:sdcnlo}
\begin{split}
	d_{\mathrm{NLO}}[g\to Q\bar{Q} (\CScPz)] =
	\frac{ \alpha_s^3 }{3 \pi N_c m_Q^5 }
	\Bigg(
	& p_\delta \delta (1-z)
	+ \sum_{i=0}^2 p_i \left[\frac{\ln^i(1-z)}{1-z}\right]_+
	+ p(z) \\
	& +\ln \left(\frac{\mu_r^2}{4m_Q^2} \right) b_0 \, d_{\mathrm{LO}}[g\to Q\bar{Q} (\CScPz)]
	\Bigg) \,,
\end{split}
\end{align}
where the coefficients $p_\delta,p_i$ can be fitted with the PSLQ algorithm from high-precision numerical values as following
\begin{align}
\begin{split}
	p_\delta & = \frac{95-24 \pi ^2}{162}  n_f
		+\frac{756 \zeta (3)+66 \pi ^2-1141-576 \ln ^2 2-1032 \ln 2}{324} N_c
		- \frac{2(4 \ln 2-1)}{3} \frac{1}{N_c} \,,\\
	p_0 & = -n_f
		+\frac{-12 \pi ^2+119+64 \ln 2 }{18}  N_c
		+ \frac{8}{3} \frac{1}{N_c} \,,\\
	p_1 & = \frac{4}{9}  n_f
		-\frac{64 }{9}  N_c \,,\\
	p_2 & = -\frac{16 }{3}  N_c  \,,
\end{split}
\end{align}
and $p(z)$ can be expressed as a piecewise function
\begin{equation}
  	p(z) =
  	\left\{
 	\begin{aligned}
 		&\frac{63-8 \pi ^2}{36 z} N_c + \sum_{i=0}^{2} \sum_{j=0}^{\infty} \ln^i z \, (2z)^j \left(A_{ij}^{f} \, n_f + A_{ij} \, N_c +  \frac{A_{ij}^{N}}{N_c} \right) \,, & \text{for } 0<z<\frac{1}{4} \\
 		&\sum_{j=0}^{\infty} (2z-1)^j \left(B_{j}^{f} \, n_f + B_{j}\, N_c + \frac{B_{j}^{N}}{N_c} \right) \,, & \text{for }    \frac{1}{4} \le z \le \frac{3}{4} \\
 		&\sum_{i=0}^{3} \sum_{j=0}^{\infty} \ln^i (1-z) \, (2-2z)^j \, \left(C_{ij}^{f} \, n_f + C_{ij} \, N_c + \frac{C_{ij}^{N}}{N_c} \right) \,. & \text{for } \frac{3}{4}<z<1
 	\end{aligned}
 	\right. \,.
\end{equation}
The coefficients $A_{ij}^{k},B_{j}^{k},C_{ij}^{k}$ can be evaluated numerically with high precision. In the attached ancillary file, we give these coefficients calculated up to $j=500$ to obtain the results with $140$-digit precision at any value of $z$.

The SDCs of $g\to Q\bar{Q} (\CScPj)+X$ with $J=1,2$ and $g\to Q\bar{Q} (\COcPj)+X$ with sum $J$ are calculated similarly, with results given as attachments in the ancillary file.

\section{Numerical Results and Discussion} \label{sec:sdcdis}

To see the effects of the NLO corrections, we choose parameters with $m_b=4.75\mathrm{~GeV}$, $N_c=3$, $n_f=4$, and $\alpha_s (\mu_r=2m_b)=0.181 $.
In Figure \ref{fig:res3pj18ren}, we plot the curves of LO SDCs and LO+NLO SDCs for the $P$-wave FFs with $\mu_0=\mu_\Lambda = 2m_b$.
\begin{figure}
\begin{minipage}{0.5\textwidth}
  \centerline{$g\to b\bar{b} (\CScPz)+X $}
  \centerline{\includegraphics[width=1\textwidth]{res3p01ren-eps-converted-to.pdf}}
\end{minipage}
\hfill
\begin{minipage}{0.5\textwidth}
  \centerline{$g\to b\bar{b} (\CScPa)+X $}
  \centerline{\includegraphics[width=1\textwidth]{res3p11ren-eps-converted-to.pdf}}
  \end{minipage}
\vfill
\vspace{5pt}
\begin{minipage}{0.5\textwidth}
  \centerline{$g\to b\bar{b} (\CScPb)+X $}
  \centerline{\includegraphics[width=1\textwidth]{res3p21ren-eps-converted-to.pdf}}
  \end{minipage}
\hfill
\begin{minipage}{0.5\textwidth}
  \centerline{$g\to b\bar{b} (\COcPj)+X $}
  \centerline{\includegraphics[width=1\textwidth]{res3pj8ren-eps-converted-to.pdf}}
\end{minipage}
\caption{ SDCs of the fragmentation functions of $g\to b\bar{b} (\state{3}{P}{J}{1/8}))+X $ at LO and NLO. The dashed line is for $d_{\mathrm{LO}}(z)\times c^{[1/8]}$, and the solid line is for $(d_{\mathrm{LO}}(z)+ d_{\mathrm{NLO}}(z)) \times c^{[1/8]}$, with scale choices $\mu_r=\mu_0=\mu_\Lambda=2m_b$, $c^{[1]}=9m_b^5$ for color-singlet and $c^{[8]}=144m_b^5/5$ for color-octet. The bands are obtained by varying the renormalization scale $\mu_r$ by a factor of 2. }
\label{fig:res3pj18ren}
\end{figure}
The overall factor $9m_b^5$ for color-singlet  and $144 m_b^5/5$ for color-octet comes from normalization factors at NLO.

The sensitivities of LO and LO+NLO FFs with respective to the renormalization scale $\mu_r$ are illustrated Fig.~\ref{fig:res3pj18ren} with bands corresponding to vary $\mu_r$ from $m_b$ to $4m_b$. We find that there are still large theoretical uncertainties after the NLO corrections.

With $\mu_r = \mu_0=\mu_\Lambda = 2m_b$, we also provide K-factors (the ratio of LO+NLO over LO) for some special values of $z$ in Table \ref{table:Kfactor} , where we find that the K-factors are negative for most values of $z$.
\begin{table}[htb!]
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}*{1}{*{5}{c}}}
\hline \hline
$z$ & $K_0^{[1]}$ & $K_1^{[1]}$ & $K_2^{[1]}$ & $K_J^{[8]}$  \\
\hline
$0.10$ & $-5.746823334$ & $-24.60770220$ & $-79.84997431$ & $-38.52164437$ \\
$0.15$ & $-2.285338484$ & $-10.32801893$ & $-36.18197374$ & $-17.35843659$ \\
$0.20$ & $-2.644346438$ & $-5.288432534$ & $-19.87280192$ & $-9.481286949$ \\
$0.25$ & $-4.379318130$ & $-2.914573185$ & $-11.89225918$ & $-5.572752664$ \\
$0.30$ & $-7.198593861$ & $-1.593420613$ & $-7.390964744$ & $-3.288936921$ \\
$0.35$ & $-11.09174086$ & $-0.7748378276$ & $-4.639679071$ & $-1.813295182$ \\
$0.40$ & $-15.40618698$ & $-0.2304903644$ & $-2.878388543$ & $-0.8000547312$ \\
$0.45$ & $-17.98353020$ & $0.1472252728$ & $-1.722389454$ & $-0.08193653748$ \\
$0.50$ & $-16.79362533$ & $0.4131713893$ & $-0.9577973668$ & $0.4303450620$ \\
$0.55$ & $-13.06191610$ & $0.5961117069$ & $-0.4588481917$ & $0.7878655480$ \\
$0.60$ & $-9.270458907$ & $0.7101131977$ & $-0.1500652088$ & $1.020569545$ \\
$0.65$ & $-6.453740962$ & $0.7593385828$ & $0.01246796303$ & $1.144600900$ \\
$0.70$ & $-4.605946556$ & $0.7391584729$ & $0.05021547782$ & $1.164601499$ \\
$0.75$ & $-3.488264428$ & $0.6339314219$ & $-0.03511661296$ & $1.072131730$ \\
$0.80$ & $-2.912955012$ & $0.4091668085$ & $-0.2661028908$ & $0.8383434112$ \\
$0.85$ & $-2.799683132$ & $-0.01116815401$ & $-0.7094089942$ & $0.3920893805$ \\
$0.90$ & $-3.226904801$ & $-0.8227021717$ & $-1.548987722$ & $-0.4556362309$ \\
$0.95$ & $-4.779112355$ & $-2.756468771$ & $-3.481118189$ & $-2.413493587$ \\
\hline \hline
\end{tabular*}
\caption{K-factors at different values of $z$. $K_j^{[c]}$ denotes the K-factor of $g\to Q\bar{Q} (\state{3}{P}{j}{c})+X$ and the symbol $J$ means the sum of $j=0,1,2$. \label{table:Kfactor}}
\end{table}

As shown in Figure \ref{fig:res3pj18ren}, LO FFs are divergent at $z=1$.
These divergences come from the $\delta(1-z)$ and $1/[1-z]_+$ terms in Eq.~\eqref{eq:sdclomu}.
Besides, the NLO FFs are negative and divergent at both $z=0$ and $z=1$. The leading behavior at $z\to 0$ is $1/z$, which means that the total fragmenting probability obtained by integrating the NLO FF over $z$ from $0$ to $1$ is divergent. But fortunately, physical cross sections are not sensitive to the behavior at $z\to 0$ after convolving FFs with partonic hard parts, which behave as $z^n$ in the small $z$ region with $n$ usually larger than $4$.
As shown in Eq.~\eqref{eq:sdcnlo}, the divergences at $z \to 1$ behavior as $\delta(1-z)$ and plus functions with leading power $[\ln^2(1-z)/(1-z) ]_+$, which comes from the expansion of $(1-z)^{-1-2\epsilon}\epsilon^{-2} $.

Obviously, the distribution of FF will not give us direct information about physical cross sections. To see the impact of the NLO results to the physical cross sections, we evaluate the $n$-th moments of the SDCs,
\begin{equation}\label{eq:moments}
	\int_{z_0}^1 \mathrm{d} z \, z^n (d_{\mathrm{LO}}(z)+ d_{\mathrm{NLO}}(z))\times c^{[1/8]} \, ,
\end{equation}
where $c^{[1]}=9m_b^5$ and $c^{[8]}=144m_b^5/5$.
We plot the numerical results of integration K-factors with different $z_0$ and $n$, which are shown in Figure \ref{fig:intz0n}.
\begin{figure}
	\begin{minipage}{0.5\textwidth}
	  \centerline{$g\to b\bar{b} (\CScPz)+X $}
	  \centerline{\includegraphics[width=1\textwidth]{intz0n3p01-eps-converted-to.pdf}}
	\end{minipage}
	\hfill
	\begin{minipage}{0.5\textwidth}
	  \centerline{$g\to b\bar{b} (\CScPa)+X $}
	  \centerline{\includegraphics[width=1\textwidth]{intz0n3p11-eps-converted-to.pdf}}
	\end{minipage}
	\vfill
	\vspace{5pt}
	\begin{minipage}{0.5\textwidth}
	  \centerline{$g\to b\bar{b} (\CScPb)+X $}
	  \centerline{\includegraphics[width=1\textwidth]{intz0n3p21-eps-converted-to.pdf}}
	\end{minipage}
	\hfill
	\begin{minipage}{0.5\textwidth}
	  \centerline{$g\to b\bar{b} (\COcPj)+X $}
	  \centerline{\includegraphics[width=1\textwidth]{intz0n3pj8-eps-converted-to.pdf}}
	\end{minipage}
	\caption{ The K-factors calculated by the integration results with Eq.~\eqref{eq:moments} for different $z_0$ and $n$. }
	\label{fig:intz0n}
\end{figure}
We find that for the same $n$, the K-factors are quite steady for different $z_0$, meaning that cross sections are dominated by FFs at large $z$.
Besides, the K-factors decrease with the increase of $n$. Especially, the color-singlet K-factors become negative for large $n$. And there may be a cancellation between the LO and the NLO FFs at about $7$-th moment for color-singlet quarkonium productions. Further more, it is clear that though the K-factors and values of FFs at each position $z=z_0$ are quite different for $j=0,1,2$ as shown in Table~\ref{table:Kfactor} and in Figure \ref{fig:res3pj18ren}, they are almost the same in the sense of integration as shown in Figure \ref{fig:intz0n}.

For the special case $z_0=0$, we can estimate the impact of NLO calculations on physical cross sections. The numerical results are shown in Table \ref{table:probability}.
\begin{table}[htb!]
	\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}} c c c c c c c}
		\hline \hline
		state          			  & SDCs$*c^{[1/8]}$ 		  & $z^2$ & $z^4$ & $z^6$ & $z^8$ & $z^{10}$ \\
		\hline
		\multirow{3}{*}{$\CScPz$} & LO & $-0.128268$ & $-0.135763$ & $-0.141953$ & $-0.147253$ & $-0.151898$ \\
			& LO+NLO & $-0.0982734$ & $-0.0544311$ & $-0.0143470$ & $0.0226121$ & $0.0570568$ \\
			& K-factor & $0.766154$ & $0.400927$ & $0.101068$ & $-0.153559$ & $-0.375626$ \\
		\hline
		\multirow{3}{*}{$\CScPa$} & LO & $-0.0860886$ & $-0.102339$ & $-0.114018$ & $-0.123176$ & $-0.130722$ \\
			& LO+NLO & $-0.0735848$ & $-0.0436727$ & $-0.0119884$ & $0.0195102$ & $0.0501917$ \\
			& K-factor & $0.854756$ & $0.426745$ & $0.105145$ & $-0.158393$ & $-0.383957$ \\
		\hline
		\multirow{3}{*}{$\CScPb$} & LO & $-0.106291$ & $-0.118021$ & $-0.126989$ & $-0.134287$ & $-0.140456$ \\
			& LO+NLO & $-0.0872125$ & $-0.0507874$ & $-0.0152633$ & $0.0187392$ & $0.0511634$ \\
			& K-factor & $0.820505$ & $0.430326$ & $0.120194$ & $-0.139546$ & $-0.364266$ \\
		\hline
		\multirow{3}{*}{$\COcPj$} & LO & $-0.101999$ & $-0.114765$ & $-0.124328$ & $-0.132024$ & $-0.138483$ \\
			& LO+NLO & $-0.211825$ & $-0.196364$ & $-0.173108$ & $-0.147720$ & $-0.121978$ \\
			& K-factor & $2.07674$ & $1.71101$ & $1.39235$ & $1.11888$ & $0.880816$ \\
		\hline \hline
	\end{tabular*}
	\caption{Moments and K-factors of SDCs.  \label{table:probability}}
\end{table}
For large $n$, the influence of the behavior at $z\to 1$ becomes more and more prominent and the NLO correction becomes more and more significant.
In addition, the K-factors of color-octet process are much larger than those of color-singlet process.

The above discussion implies that resummation at large $z$ region may be crucial. This can be done within the framework of soft gluon factorization \cite{Ma:2017xno}, which we leave for future study.

\input{appendix.tex}

\acknowledgments

We thank Xiao Liu and Hao-Yu Liu for useful discussions.
The work is supported in part by the National Natural Science Foundation of China (Grants No. 11475005, No. 11075002, and No. 11875071), the National Key Basic Research Program of China (No. 2015CB856700), and High-performance Computing Platform of Peking University.

\input{paper.bbl}

\end{document}